PRICKLE1 × FOCAD Interaction Revealed by Genome-Wide vQTL Analysis of Human Facial Traits

The human face is a highly complex and variable structure resulting from the intricate coordination of numerous genetic and non-genetic factors. Hundreds of genomic loci impacting quantitative facial features have been identified. While these associations have been shown to influence morphology by altering the mean size and shape of facial measures, their effect on trait variance remains unclear. We conducted a genome-wide association analysis for the variance of 20 quantitative facial measurements in 2,447 European individuals and identified several suggestive variance quantitative trait loci (vQTLs). These vQTLs guided us to conduct an efficient search for gene-by-gene (G × G) interactions, which uncovered an interaction between PRICKLE1 and FOCAD affecting cranial base width. We replicated this G × G interaction signal at the locus level in an additional 5,128 Korean individuals. We used the hypomorphic Prickle1^Beetlejuice (Prickle1^Bj) mouse line to directly test the function of Prickle1 on the cranial base and observed wider cranial bases in Prickle1^Bj/Bj. Importantly, we observed that the Prickle1 and Focadhesin proteins co-localize in murine cranial base chondrocytes, and this co-localization is abnormal in the Prickle1^Bj/Bj mutants. Taken together, our findings uncovered a novel G × G interaction effect in humans with strong support from both epidemiological and molecular studies. These results highlight the potential of studying measures of phenotypic variability in gene mapping studies of facial morphology.

Introduction

Human genetic studies have been remarkably successful at detecting variants with an impact on mean trait values. On the other hand, the genetic basis of phenotypic variability remains largely uncharacterized. Variance quantitative trait loci (vQTLs) are genetic variants exhibiting inter-individual intra-genotypic variability, where one of the alleles is associated with a larger phenotypic variance compared to the other (Rönnegård and Valdar, 2012). vQTLs are seldom examined in conventional genome-wide association studies (GWASes), which usually assume variance homogeneity across genotype groups and aim to detect differences in group means. Despite this, knowledge of vQTLs may lead to a deeper understanding of genotype–phenotype relationships by providing insights into the genetic control of phenotypic variance and revealing possible interactions among genetic variants.

The human face is a highly complex structure resulting from the intricate coordination of multiple genetic and epigenetic factors. Beyond the individual effect of common single nucleotide polymorphisms (SNPs) identified through GWAS (Liu et al., 2012; Paternoster et al., 2012; Adhikari et al., 2016; Cole et al., 2016; Shaffer et al., 2016; Lee et al., 2017; Cha et al., 2018; Claes et al., 2018; Crouch et al., 2018; White et al., 2020), gene-by-gene (G × G) interactions are also expected to play an essential role. Statistical interaction is typically defined as the deviation from additivity in a linear model (Moore and Williams, 2009). The biological interpretation of a statistical interaction is often a challenging task and requires careful experiments to demonstrate a real crosstalk between the genetic factors in the relevant context (Moore, 2005; Moore and Williams, 2005). Although recognized as one of the critical potential sources of missing heritability (Zuk et al., 2012), G × G interactions remain largely uncharacterized in most human traits due to low statistical power, a high computational burden, and difficulties in uncovering the biological relevance for statistical interactions. A recent study of facial morphology reported some evidence of the potentially coordinated actions of variants on facial variation (White et al., 2020). We now know that facial morphology is influenced by variants at over 200 loci (White et al., 2020), which, in combination, give rise to a vast and usually unknowable interaction map. As a genome-wide search for G × G interactions would be intractable, narrowing the focus to a subset of candidate SNPs is necessary, and effective strategies are needed to guide this pre-selection process.

One sensible strategy is variance prioritization. It has been shown, both theoretically and empirically, that a heterogeneous phenotypic variance across genotype groups can arise from genetic interaction effects (Struchalin et al., 2010; Forsberg and Carlborg, 2017). Although these two non-additive genetic inheritance patterns are not always related, their relationship suggests an enrichment of the interaction effects among variance-controlling SNPs at vQTLs, indicating that these SNPs can be prioritized for interaction tests. Indeed, this variance prioritization strategy has been shown to be useful in studying other quantitative human phenotypes (Paré et al., 2010; Hulse and Cai, 2013; Rask-Andersen et al., 2017; Sarkar et al., 2019; Wang et al., 2019; Young et al., 2018; Johnson et al., 2020; Marderstein et al., 2021).

The aim of this study was twofold: to characterize the genetic control of facial trait variability and to investigate G × G interaction effects by leveraging the identified vQTLs. We conducted genome-wide vQTL analyses for 20 facial measurements in 2,447 individuals. Among the suggestive vQTLs, one was located within PRICKLE1, a gene known to contribute to the control of craniofacial development in mouse mutants (Yang et al., 2013; Liu et al., 2014; Gibbs et al., 2016; Wan et al., 2018). We followed up on this vQTL by testing whether it interacts with other loci across the genome. We discovered a novel G × G interaction effect between PRICKLE1 and FOCAD impacting cranial base width, which was statistically replicated in an independent cohort and experimentally verified by showing the co-expression of both genes during the critical stages of mouse craniofacial development. Our findings highlight the importance of studying facial variability in addition to facial mean differences for gene discovery and mechanistic exploration.

Materials and Methods

Discovery Cohort

The cohort, facial phenotyping, and SNP genotyping were described in a previous GWAS (Shaffer et al., 2016). In brief, three-dimensional (3D) facial images were collected using digital stereophotogrammetry on a cohort consisting of 2,447 unrelated white individuals of European ancestry from the United States. The cohort ranged in age from 3 to 49 years, had a female proportion of 62%, and was free of any conditions known to affect the face or head. Facial landmarks were manually registered on the surface images, from which we calculated 20 linear distances (Supplementary Figure 1 and Supplementary Table 1). The cohort was genotyped on the Illumina OmniExpress + Exome v1.2 array (Illumina, San Diego, CA, United States) and was then fully imputed to the 1,000 Genomes Project phase 3 reference panel. Ancestry principal components (PCs) based on linkage disequilibrium (LD)-pruned autosomal SNPs were constructed and used to control for population stratification.

SNPs and Significance Threshold in the Discovery Analysis

Detecting differences in group variances demands larger sample sizes than when detecting differences in group means. Given that our cohort is smaller than the typical sample size of previous vQTL studies, we increased the minimum minor allele frequency (MAF) cutoff for the sake of maintaining statistical power. SNPs with MAF <0.2 were excluded, giving a total of 3,104,639 qualified autosomal SNPs (genotyped and imputed). Because SNPs with MAF >0.2 represent much less of the genome-wide variation than typically interrogated in a GWAS, to reduce burden of multiple comparisons, we used a relaxed threshold of 5 × 10^–7 to identify suggestive hits for the sake of not missing potential signals. We calculated the empirical type I error rate at the 5 × 10^–7 level after randomly shuffling the phenotypic data (see section “Results”), and as we expect some false associations, to prioritize SNPs for the G × G interaction stage, we further considered the biological relevance of the loci where the vQTLs are located rather than relying on their p-value alone.

Statistical Methods

We used the Levene’s test with median implemented in the OSCA software package¹ to identify SNPs associated with the variance of 20 quantitative facial measurements. The Levene’s test is equivalent to a test of differences in mean deviation across subgroups, with the deviation calculated as the absolute difference between an individual’s phenotypic value and their group-specific mean/median. We used the median as the central measure because it provides more robustness to non-normality and outliers than does the mean. To control for covariates, age, age², sex, height, weight, facial size, and four genetic ancestry PCs were first regressed on facial traits in a linear regression model, and the residuals were then used as the phenotype in statistical tests.

The G × G interaction test was done in PLINK (Purcell et al., 2007) with the command –epistasis, which uses a linear regression to model phenotypic residuals on the main effects of each inspected SNP and an interaction term for the pair. Since interaction tests at the vQTL SNP would inevitably lead to a violation of the assumption of variance homogeneity in the generalized linear model, we examined Q–Q plots of the interaction p-values, and when there was evidence of inflation (quantified by the genomic inflation factor, λ), we used the robust standard error to obtain unbiased standard errors of the coefficients under variance heteroscedasticity (R package sandwich). This technique was able to correct the inflation and generate well-behaved Q–Q plots. In the discovery analysis, all G × G interaction analyses were performed between a single SNP and all other SNPs across the genome, and therefore a p-value below the conventional threshold of 5 × 10^–8 was considered significant. Co-localization plots for the vQTL analysis and the G × G interaction summary statistics were generated by R package LocusCompareR (Liu et al., 2019).

Although the current study had a primary focus on G × G interactions as the underlying genetic process for variance heterogeneity, we also considered possible gene–environment (G × E) interactions by using a linear regression model in PLINK. Both the 1 degree-of-freedom (df) test for the interaction effect only and the 2 df test for the main and the interaction effects jointly were used. Only a limited set of environmental factors were collected in our cohort, and we chose to examine the interactions of SNPs with age and sex.

Replication

The replication cohort consisted of 5,128 unrelated individuals of Asian ancestry collected for the Korean Genome and Epidemiology Study (KoGES) in South Korea from 2009 to 2012 (Ansan–Ansung cohort), for which standardized two-dimensional (2D) frontal photographs were available. The cohort, phenotyping approach, and genotype data have been described in detail previously (Cha et al., 2018). The maximum upper width of the face was represented by the zygion-to-zygion distance (zyR–zyL), which, in a 2D representation, is approximately equivalent to the 3D cranial base width measurement in the discovery sample. ZyR–zyL was pre-adjusted for age, sex, height, and weight and was analyzed using the same statistical approaches applied in the discovery stage. We considered first a SNP-level replication where the exact same SNPs identified in the discovery analysis or, when unavailable, their proxies identified based on LD were tested in the replication data. Next, we sought for a locus-level replication where the SNPs located within 500 kb on either side of the discovery signals were considered. The significance threshold was determined by the Bonferroni method correcting for the effective number of tests (M_eff), which was computed according to the eigenvalue-based procedure by Li and Ji (2005), informed by the 1,000 Genomes Project East Asian LD structure. For the locus-level replication of the G × G interaction signal, a p-value below 0.05 divided by the product of the M_eff of the two regions under test was considered significant.

Prickle1^Bj Mouse Husbandry and Sample Preparation

The Prickle1^Bj mouse line is maintained as heterozygotes outbred to wild-type C57/Bl^6J individuals. All animal work is approved by the University of Pittsburgh Institutional Animal Care and Use Committee. To collect embryos, we performed timed mating, where the presence of the plug was designated as embryonic day 0.5 (E0.5). On the appropriate day of pregnancy, we euthanized the pregnant dam and collected the embryos via caesarian section. We confirmed the embryonic development stage based on morphology. Genotyping was performed using custom SNP Taqman assays (AH7041R, Invitrogen, Waltham, MA, United States).

Embryos were fixed overnight in 4% paraformaldehyde, followed by dehydration and storage in 70% ethanol until further processing. For the tissue sections, the embryos were embedded in paraffin wax, sectioned at 8-μm thickness, and placed on (3-aminopropyl) triethoxysiloxane (TESPA)-treated slides that were stored at 4°C until use.

Mouse Morphometrics

The superior view images of 7 Prickle1^+/+, 10 Prickle1^Bj/+, and 13 Prickle1^Bj/Bj heads at E14.5 were captured on a Leica M165FC dissecting microscope with a DFC 450 camera at 1.6× magnification (Leica, Wetzlar, Germany). Care was taken to avoid excess canting of the heads during imaging. Six landmarks were placed around the cranium using 3D Slicer (Fedorov et al., 2012): two bilaterally at the eyes, two bilaterally at the ears, and one each at the midline on the anterior and posterior borders of the cranium. Using MorphoJ (Klingenberg, 2011), the landmark configurations were subjected to generalized Procrustes superimposition in order to translate, rotate, and scale them to unit centroid size and to facilitate relative comparison of cranial base width by genotype. A linear distance capturing cranial base width was calculated between the left and right landmarks placed at the junction of the posterior end of the eye globes where they meet the cranium (Figure 2A); this distance approximated our measure of cranial base width in humans. Linear regression and Levene’s test were then performed to test for differences in the mean and variance, respectively, of cranial base width across the three genotype groups under the additive genetic model while adjusting for litter and centroid size (as a global measure of head size) as covariates. Analysis was performed in the R statistical environment. Differences were considered statistically significant at p < 0.05.

Mouse Immunofluorescence Staining and Imaging

For double immunofluorescence labeling, we used the anti-Focadhesin (HPA055015, Sigma, St. Louis, MO, United States) and the anti-Prickle1 (sc-393034, Santa Cruz, Santa Cruz, CA, United States) antibodies and performed citrate buffer antigen retrieval. After antigen retrieval, the slides were briefly blocked and then incubated overnight in the fridge with the primary antibodies (diluted 1:100). The next day, the slides were washed followed by incubation with secondary antibodies (Alexa Flour, 1:200). The sections were mounted with Prolong Gold with DAPI (Invitrogen). Fluorescent images were captured on a Nikon TE2000 inverted fluorescent microscope (Nikon, Tokyo, Japan) and are represented as a single slice of an image stack. Prickle1^+/+ and Prickle1^Bj/Bj images were captured at the same exposure and magnification.

Ethics Statement

Institutional ethics (IRB) approval was obtained at each recruitment site (University of Pittsburgh Institutional Review Board #PRO09060553 and #RB0405013; UT Health Committee for the Protection of Human Subjects #HSC-DB-09-0508; Seattle Children’s Institutional Review Board #12107; University of Iowa Human Subjects Office/Institutional Review Board #200912764 and #200710721; and the Korea Institute of Oriental Medicine #I-2007/006-002). All adult subjects gave written informed consent prior to participation; for children, written consent was obtained from a parent or legal guardian. All procedures performed in this study were conducted in accordance with the guidelines of the Declaration of Helsinki. All experimental protocols using mice were approved by the University of Pittsburgh Institutional Animal Care and Use Committee (#17050839) and carried out in accordance with institutional animal care protocols.

Results

vQTL Search for 20 Facial Measurements Identified PRICKLE1

In a cohort of 2,447 unrelated individuals with 3D facial surface data, we conducted a genome-wide search of vQTLs for 20 facial measurements (Supplementary Figure 1 and Supplementary Table 1) using the Levene’s test of medians. No genome-wide significant vQTLs were identified. Ten suggestive loci with p-values below 5 × 10^–7 were observed (Table 1). None of these loci have been implicated in prior facial GWAS or candidate gene association studies. Most of the lead SNPs are intergenic, and the nearby genes are not known to have roles in craniofacial morphogenesis. However, 12q12 showed evidence of association with variance in cranial base width, and its lead SNP rs1796391 was located in an intron of PRICKLE1 (Figure 1). The Levene’s test statistics for cranial base width were well-behaved and there was no sign of genomic inflation or deflation (Figure 1A). PRICKLE1, when disrupted, has been linked to craniofacial malformations, such as orofacial clefting (OFC), in animal models and in humans (Yang et al., 2013, 2014; Wan et al., 2018; Ahsan et al., 2019). Interestingly, excessive cranial base/upper facial width has been reported as a feature in families with a history of OFC (Fraser and Pashayan, 1970; Nakasima and Ichinose, 1983; Blanco et al., 1992; Mossey et al., 1998; Suzuki et al., 1999; Chatzistavrou et al., 2004; McIntyre and Mossey, 2004; Yoon et al., 2004; Weinberg et al., 2006, 2008; Roosenboom et al., 2017).

TABLE 1

Table 1. Suggestive variance quantitative trait loci (vQTLs) of facial features (p < 5 × 10^–7).

FIGURE 1

Figure 1. Results of the genome-wide variance quantitative trait locus (vQTL) analysis for cranial base width. (A) Q–Q plot and Manhattan plot. There was no sign of inflation. The arrow points to the PRICKLE1 locus on chr12. The blue horizontal line indicates suggestive threshold 5 × 10^{– 7}. The bottom two plots show the association –log₁₀(p-value) at the PRICKLE1 locus with the variance (B) and the mean (C) of cranial base width.

In light of the potential function of PRICKLE1 in facial development, we decided to directly test its role using a hypomorphic Prickle1 mouse allele. The Prickle1^Bj mouse allele has a missense mutation in Prickle1 at c:G482T (p:C161F) (Gibbs et al., 2016; Wan et al., 2018). This mutation causes the substitution of cysteine with a phenylalanine in a cysteine knot that is required for protein–protein interactions. The Prickle1^Bj allele is a hypomorphic allele and bypasses the pre-gastrulation lethality of the conventional knockout allele (Tao et al., 2009). Prickle1^Bj/Bj survivors die shortly after birth due to cardiovascular and craniofacial anomalies (Gibbs et al., 2016; Wan et al., 2018).

We explored to what extent the Prickle1^Bj allele affected the cranial base width in a dosage-dependent fashion using data from 30 mice. At E14.5, Prickle1^Bj/Bj mice exhibited significantly wider cranial bases (p = 0.01) compared to Prickle1^+/+ and Prickle1^Bj/+ littermates (Figure 2). The variance of this phenotype in the three genotype groups was not significantly different (p = 0.72), which could be due to a lack of power given the limited number of mice.

FIGURE 2

Figure 2. The Prickle1^Bj mouse allele dosage was associated with a wider cranial base. (A–C) Superior views of littermates at embryonic day 14.5 showing craniofacial shape changes in Prickle1^Bj/Bj compared to Prickle1^+/+. Cranial base width measurement is indicated by the black line in (A). Scale bar, 1 mm (and applies to A–C). (D) Box plot of the cranial base width. Thick horizontal lines represent group medians. The p-value was obtained from the test of group mean differences assuming an additive genetic model, adjusting for head size and litter.

We observed 11 SNPs with a p < 5 × 10^–7 after randomly shuffling the cranial base width among our samples (Supplementary Figure 2), indicating that some of the 10 identified suggestive vQTLs are likely spurious associations. We therefore selected PRICKLE1 rs1796391, the only signal which has been previously studied in craniofacial disorder, as the prioritized locus to be further examined in the G × G interaction stage and did not pursue an analysis of the other nine.

G × G Interaction of PRICKLE1 and Cranial Base Width

We followed up on the human PRICKLE1 vQTL to explore the source of the phenotypic variance heterogeneity, with a particular focus on its G × G interaction effect. We utilized the robust standard error (see section “Materials and Methods”) to ensure well-behaved statistics for the genome-wide search of SNPs interacting with the lead PRICKLE1 SNP rs1796391 (Supplementary Figure 3A). While no SNP pairs reached the genome-wide significance threshold (Supplementary Table 2 and Supplementary Figure 3B), the top interaction signal was seen at the gene FOCAD on chromosome 9, with its intronic SNP rs10511683 interacting with the PRICKLE1 vQTL rs1796391 (p = 5.82 × 10^–7) (Supplementary Figure 3C). Figure 3 illustrates the form of this G × G interaction effect and its relationship with the vQTL signal at the PRICKLE1 locus. The minor allele homozygote group (Figure 3A, rightmost) of PRICKLE1 rs1796391 displayed the largest phenotypic variance, which is a result of the differential phenotypic means when stratified by the genotype at the interacting FOCAD SNP rs10511683 (Figure 3B). The right table in Figure 3B displays the association effect sizes and the p-values of rs10511683 (FOCAD) for cranial base width residuals in the unstratified and the stratified data. The minor allele of rs10511683 (FOCAD) was associated with a wider cranial base in people with the AA genotype at rs1796391 (PRICKLE1), whereas it showed an opposite effect of similar magnitude in the PRICKLE1 GG group. FOCAD SNP rs10511683 had small p-values for its association with the phenotypic mean in the stratified analyses of the AA and GG groups of rs1796391 (PRICKLE1), but not in the combined sample, where the opposite effects offset each other (p = 0.13) (also see Supplementary Figure 3D).

FIGURE 3

Figure 3. The cranial base width variance heterogeneity at rs1796391 (PRICKLE1) was induced by its gene-by-gene (G × G) interaction effect with rs10511683 (FOCAD) on phenotypic mean. (A) Box plots show the residual of cranial base width [after regressing out age, age², sex, height, weight, facial size, and four genetic principal components (PCs)] in the three genotype groups of the variance quantitative trait locus (vQTL) rs1796391 (PRICKLE1), colored by the genotype at the interacting SNP rs10511683 (FOCAD). Standard deviations are shown above boxes. (B) Box plots show the distribution of the phenotypic residual in all nine combinations of the genotypes at rs1796391 (PRICKLE1) and rs10511683 (FOCAD). The x-axis, y-axis, and the color scheme are the same as those in (A). Medians are represented by diamonds and connected to form segments, the slopes of which indicate the group-specific effects of rs10511683 (FOCAD). The non-parallel pattern is the hallmark of statistical interaction. Table to the right shows the association (with the mean of cranial base width) beta coefficients, p-values, and sample sizes for rs10511683 (FOCAD) in subgroups defined by the genotype at rs1796391 (PRICKLE1) and the combined sample.

We queried the human cell line experimental data from the Roadmap Epigenomics Project and the ENCODE Project for potential regulatory roles of the identified non-coding SNPs (Supplementary Figure 4). The genomic sequence enclosing rs10511683 (FOCAD) overlaps with weak signals of histone modifications, signifying enhancer activity in chondrocyte cells and osteoblast primary cells. SNP rs1796391 (PRICKLE1) also showed a potential enhancer signature in osteoblasts. These cell types are important players during craniofacial morphogenesis, and studying cell type-specific roles of PRICKLE1 and FOCAD in future work could potentially provide more insights.

The Relationship Between Variance Heterogeneity and G × G Interaction at the PRICKLE1 Locus

We provided two pieces of evidence supporting that this G × G interaction effect can at least partially account for the observed variance heterogeneity. Firstly, a scatter plot with the local LD pattern taken into account demonstrated a co-localization of the vQTL and G × G interaction test summary statistics at the FOCAD locus (Figure 4), in agreement with the hypothesis that the interaction effect underlies the observed variance heterogeneity. Secondly, a stratified vQTL analysis was performed for PRICKLE1 SNPs in each of the three genotype groups defined by rs10511683 (FOCAD) separately. If the variance heterogeneity at rs1796391 (PRICKLE1) was induced by its interaction with rs10511683 (FOCAD), we would expect to see weakened signals in each group. In accordance with this expectation, the group-specific vQTL signals were either absent or much attenuated compared to the unstratified analysis (Supplementary Figure 5). Specifically, the p-values of rs1796391 were 1.83 × 10^–4, 0.05, and 0.6, respectively (Supplementary Table 3). We further downsampled the two larger strata and repeated the analysis in three equally sized groups, the results of which are shown in Supplementary Table 3. We caution that the stratification inevitably led to a reduced power and a tendency for higher variance in smaller groups, which complicated the interpretation of this result. Nonetheless, the results from these two analyses suggest that the interaction with FOCAD SNPs is a viable explanation for the observed variance heterogeneity at the PRICKLE1 locus with no contradicting evidence.

FIGURE 4

Figure 4. LocusCompare plot at the PRICKLE1 locus for cranial base width. (Left) Joint distribution of the –log₁₀(p-value) from the Levene’s test and the gene-by-gene (G × G) interaction test, colored by linkage disequilibrium (LD) with rs1796391. Co-localization of the signals is supported by the clustering of dots at top right. (Top right) –log10(p-value) from the Levene’s test on the y-axis and genomic coordinate on the x-axis. (Bottom right) -log₁₀(p-value) from the G × G interaction test between rs10511683 (FOCAD) and SNPs at the PRICKLE1 locus on the y-axis.

G × G Interaction of FOCAD for Cranial Base Width

We then turned our focus to FOCAD rs10511683 and asked whether it also interacts with other loci. A genome-wide search did not identify significant interacting loci other than PRICKLE1. However, we did find a second interaction signal at the PRICKLE1 locus in close proximity to (7 kb away), but demonstrating low LD (r² = 0.11) with the original signal at rs1796391 (Supplementary Table 4 and Figure 5). The lead SNP of this second peak, rs10880322 (PRICKLE1), had a significant interaction p-value (1.88 × 10^–10) smaller than that of the original one. SNP rs10880322 had neither a strong mean effect (GWAS p = 0.025) nor a variance effect (Levene’s test p = 0.428) on cranial base width when tested alone; however, when taking into account its interaction with rs10511683 (FOCAD), both the main and interaction effects of rs10880322 became highly significant (Supplementary Table 5 and Supplementary Figure 6). Note that seeing a significant interaction effect with neither SNP being a vQTL, such as the case of rs10880322 × rs10511683, is not a paradox because not all forms of interactions would manifest as variance heterogeneity. We also note that, since the assumption of variance homogeneity was not violated at rs10511683 and there was no sign of inflation in the Q–Q plot (Figure 5A), we need not have used the robust standard error technique in the interaction model, and this gave a smaller p-value (3.36 × 10^–8) for the original rs1796391 (PRICKLE1) × rs10511683 (FOCAD) interaction than that reported in the section G × G interaction of PRICKLE1 for cranial base width (5.82 × 10^–7) above (where the robust standard error technique was used).

FIGURE 5

Figure 5. Results of a genome-wide gene-by-gene (G × G) interaction search for loci interacting with FOCAD rs10511683. (A) Q–Q plot and Manhattan plot. The red horizontal line indicates 5 × 10^{– 8}. (B) Regional association at the PRICKLE1 locus, colored by the linkage disequilibrium (LD) with the suggestive variance quantitative trait locus (vQTL) rs1796391. (C) The same region as in (B), colored by the LD with the other lead SNP rs10880322.

Other Possible Mechanisms for the Variance Heterogeneity at the PRICKLE1 Locus

Interaction is only one of the possible explanations for variance heterogeneity; other genetic processes can also give rise to genotype-specific variances. One group of such processes reflects genuine genetic effects, such as the presence of linked causal variants with a mean effect, nearby rare variants with a mean effect, and G × E interactions (Ek et al., 2018). Another group includes statistical artifacts that may confound the analysis and should be managed properly, such as the presence of phenotypic outliers. In our analysis of cranial base width, we noticed a potential outlying measurement for one participant and thus conducted a sensitivity analysis to evaluate the robustness of the signals at rs1796391 (PRICKLE1). The outlier is visible in Figure 3A, where one individual in the minor allele homozygous group had the maximum phenotypic residual. After removing this individual, rs1796391 (PRICKLE1) yielded a p = 9.03 × 10^–7 in the Levene’s test, and its interaction with rs10511683 (FOCAD) had a p = 1.07 × 10^–6 (using robust standard error). Supplementary Figure 7 is a version of Figure 3 with the outlying observation excluded. As expected, both signals were only slightly weakened compared to those in the original analyses, in line with the robustness of the Levene’s test with median to non-normality and outliers.

We further examined alternative explanations involving three different kinds of genetic processes. Firstly, SNP rs1796391 (PRICKLE1) did not show a mean effect on cranial base width (p = 0.18), nor did any common SNP (MAF > 1%) located within ±500 kb have a small p-value (Figure 1C). Secondly, to examine whether there were nearby low-frequency SNPs (0.02% < MAF < 1%) associated with the phenotype, we conducted gene-based tests for PRICKLE1 coding variants using the sequence kernel association test (SKAT) (Lee et al., 2012) and the combined multivariate and collapsing (CMC) method (Li and Leal, 2008). There was no evidence of association (SKAT p = 0.83, CMC p = 0.32), although variants that are rarer or undetected could still potentially have an effect. Thirdly, the results from the G × E tests did not support an interaction effect of rs1796391 with either sex or age, which were the only two non-genetic factors available in this cohort (Supplementary Table 6). These results indicated that the genotype-specific variances at the PRICKLE1 locus were unlikely to be induced by the surrounding mean-controlling SNPs, nor could they be possibly explained by interaction effects with the two non-genetic factors. It should be noted that, theoretically, one cannot exhaust all factors that can lead to heterogeneous variance by assorted mechanisms. For example, just like common and rare SNPs, the genetic mean effect of copy number variants could also confound a variance analysis at nearby loci. Moreover, there are likely plenty of non-genetic factors influencing facial morphology by interacting with individuals’ genetic composition, yet current knowledge is limited and the data necessary for examining how they played a role in our analysis are not available. Therefore, our goal was not to exhaustively scrutinize all possible explanations for the PRICKLE1 vQTL, but rather to identify the most likely ones as far as the data allow and further verify them in replication and functional experiments.

Replication Strategy

The replication cohort consisted of 5,128 unrelated individuals of Asian ancestry collected in South Korea for which a suitable measure of maximum upper facial width was available. We conducted both a SNP-level and a locus-level replication of the PRICKLE1 vQTL signal and the interaction between PRICKLE1 and FOCAD. The SNP-level replication considered the exact same SNPs detected in the discovery stage or, when unavailable, proxy SNPs based on LD. The locus-level replication expanded the genetic markers under consideration using a window size of 500 kb on either side of the genes of interest. We used a minimum MAF of 0.1 for SNP inclusion, which is lower than that of the discovery analysis as the replication cohort is roughly twice as large. For simplicity and clarity, below, we follow a SNP rs identifier with either “PRICKLE1” or “FOCAD” in parentheses to indicate its locus. Note that this notation does not imply that the SNP falls inside the gene itself, but reflects our locus-level strategy considering a surrounding region (approximately 1 Mb) encompassing the gene. A total of 906 SNPs in and nearby PRICKLE1 and 796 SNPs in and nearby FOCAD were available for replication analysis. The numbers of independent SNPs were estimated to be 138 and 162, respectively, according to the eigenvalue-based procedure by Li and Ji (2005). The significant p-value thresholds were therefore set to be 0.05/138 = 3.62 × 10^–4 for replicating the PRICKLE1 vQTLs and 0.05/138/162 = 2.24 × 10^–6 for replicating the PRCKLE1 × FOCAD effect.

Replication of PRICKLE1 vQTLs

No PRICKLE1 locus SNPs yielded a significant replication p-value in the Levene’s test (Supplementary Figure 8 and Supplementary Table 7). The discovery lead SNP, rs1796391, was not available in this replication cohort. We identified four proxy SNPs in perfect LD with rs1796391 (r² = 1) (Supplementary Table 8), yet none of them showed evidence of influencing phenotypic variance (p > 0.1 for all) (Supplementary Table 9).

SNP-Level Replication of G × G Interactions

Similarly, the lead SNP rs10511683 (FOCAD) interacting with PRICKLE1 in the discovery cohort was not available in the replication cohort. The best proxy SNPs had only moderate LD with rs10511683 (r² = 0.63–0.77) (Supplementary Table 8). Supplementary Table 10 displays the G × G interaction test p-values between each pair of the proxy SNPs in the two genes. There was no evidence of replication when only these 16 SNP pairs were taken into consideration.

Locus-Level Replication of G × G Interactions

The distribution of the 721K valid p-values from the G × G interaction test was consistent with the null expectation, except for those at the tail (Supplementary Figure 9), indicating the presence of significant interactions between the loci at some fine-scale locations. Fifteen pairs of SNPs passed our significance threshold of 2.24 × 10^–6 (Supplementary Table 11), and these interactions were mainly seen between one of two FOCAD locus SNPs (rs10964862 and rs10123324) and the different PRICKLE1 locus SNPs (Table 2). SNP rs10964862 (FOCAD) showed significant interaction with 10 PRICKLE1 SNPs, and the top two interacting SNP pairs surpassed the conventional genome-wide threshold of 5 × 10^–8: rs10964862 (FOCAD) × rs11181736 (PRICKLE1) with p = 4.53 × 10^–9 and rs10964862 (FOCAD) × rs11181735 (PRICKLE1) with p = 9.33 × 10^–9. These interaction signals were several hundred kilobase pairs away from the original SNP pair identified in the discovery samples (rs1796391 × rs10511683), and the signal at the PRICKLE1 locus did not overlap with the known coding sequence (Figure 6A). The G × G replication lead SNP rs10964862 (FOCAD) is 9.4 kb toward the 5′-UTR of the IFNW1 gene and overlaps with a sequence showing enhancer signatures in mesenchymal stem cell-derived chondrocyte cells (Supplementary Figure 4). The other SNP at the FOCAD locus showing several interactions was rs10123324, an intronic SNP of the MLLT3 gene. It had the lowest interaction p-value with a long non-coding RNA intronic SNP nearby PRICKLE1. Signals involving this SNP (rs10123324) were also a few hundred kilobase pairs distant to the target SNP pair (Figure 6B).

TABLE 2

Table 2. Locus-level replication analysis of FOCAD × PRICKLE1 detected significantly interacting SNP pairs between one of the two SNPs (rs10964862 and rs10123324) at the FOCAD locus and different PRICKLE1 SNPs.

FIGURE 6

Figure 6. Significant loci in the locus-level replication analysis of FOCAD × PRICKLE1. (A) –log10(p-value) of the interaction between rs10964862 (FOCAD locus) and single nucleotide polymorphisms (SNPs) at the PRICKLE1 locus in the replication cohort. (B) –log10(p-value) of the interaction between MLLT3 intronic rs10123324 and SNPs at the PRICKLE1 locus in the replication cohort.

Replication in the United States Cohort of Interactions With rs10964862 Identified in the South Korean Cohort

We tried to replicate the interacting SNP pairs identified in the South Korean cohort in our United States sample. One of the G × G replication SNPs, rs10123324 (FOCAD), failed the quality control procedure in the United States cohort. Since most of the SNPs that the other G × G replication, rs10964862 (FOCAD), interacted with (Supplementary Table 11) were unavailable and had no suitable proxy, we conducted a PRICKLE1 locus-wide test for rs10964862. Out of the 1,405 SNPs tested, the most significant interaction was found for rs10736004 (PRICKLE1) with p = 6.0 × 10^–4, which did not pass a Bonferroni threshold correcting for the effective number of tests [120; by the method of Li and Ji (2005)]. Ten SNPs yielding the lowest interaction p-values can be found in Supplementary Table 12.

Co-localization of Prickle1 and Focadhesin in the Mouse Cranial Base

To locate Prickle1 and Focadhesin (the protein encoded by Focad) in the cranial base, we performed immunofluorescence staining with antibodies to Prickle1 and Focadhesin at E12.5 and E15.5 cranial base. Focad was observed ubiquitously throughout the Prickle1^+/+ cranial base cells at E12.5 (Figure 7A). In contrast, it was enriched in the puncta in the Prickle1^Bj/Bj cranial base (Figure 7C). In the more mature chondrocytes at E15.5, both Prickle1 and Focadhesin proteins were enriched near the cleavage furrow of the Prickle1^+/+ proliferating chondrocytes (Figure 7B). Consistent with the E12.5 results, in the E15.5 Prickle1^Bj/Bj cells, Focadhesin was enriched in the puncta and was not consistently found in the cleavage furrow (Figure 7D). These data suggest that Prickle1 and Focadhesin are co-localized during normal chondrogenesis and that Prickle1 protein contributes to the control of the Focadhesin protein localization during cranial base development.

FIGURE 7

Figure 7. Prickle1 and Focadhesin proteins co-localize in the mouse cranial base mesenchyme. Dual immunofluorescence images of Prickle1 (green) and Focadhesin (red) proteins, with the nuclei stained with DAPI (blue). The areas within the white boxes are presented as single channels beside the combined image. (A,C) Embryonic day 12.5 cranial bases: in the Prickle1^+/+ mesenchyme, the Prickle1 and Focadhesin proteins are co-localized (yellow, arrows), while in the Prickle1^Bj/Bj mesenchyme, the signals are co-localized in bright foci, but in the cell periphery. (B,D) In Prickle1^+/+ chondrocytes, the signals co-localize (yellow, arrows) at the cleavage furrow. In Prickle1^Bj/Bj chondrocytes, Focadhesin is found in numerous puncta and in some cleavage furrows.

Discussion

By focusing on the variance of human facial features, this study uncovered a novel G × G interaction effect between PRICKLE1 and FOCAD in shaping the cranial base/upper facial width and provided strong biological support for this statistical interaction by showing the co-expression of the two genes in relevant craniofacial tissues during early development. Our findings suggest that future genetic studies of human facial morphology should expand their perspective from a sole focus on phenotypic mean to include measures of phenotypic variance and to move beyond the marginal effects of genetic variants in order to investigate their joint and interaction effects.

PRICKLE1 serves as a signaling factor in the non-canonical Wnt pathway, the disruption of which is known to cause cleft palate and to stunt limb growth (He and Chen, 2012; Mostowska et al., 2012; Yang et al., 2013; Liu et al., 2014; Gibbs et al., 2016; Wan et al., 2018). A recent study characterized the functions of Prickle1a and Prickle1b in zebrafish cranial neural crest cell development during epithelial-to-mesenchymal transition and migration. Studies in mice also support its essential role in craniofacial development. Prickle1 missense allele mutant mice were microcephalic and displayed several craniofacial defects including a cleft lip, incompletely penetrant cleft palate, and a shorter proximal–distal axis of the head (Wan et al., 2018). These phenotypes were a result of the abnormal migration and differentiation of osteoblast precursors in the frontal bone in the absence of a functioning Prickle1 protein. Sequencing studies of patients with craniofacial syndrome implicated both rare and common PRICKLE1 variants (Yang et al., 2014). In accordance with the facial measurements significantly associated in the current study, an enlarged upper facial width was frequently reported in patients with OFC and their unaffected relatives (Fraser and Pashayan, 1970; Nakasima and Ichinose, 1983; Blanco et al., 1992; Mossey et al., 1998; Suzuki et al., 1999; Chatzistavrou et al., 2004; McIntyre and Mossey, 2004; Yoon et al., 2004; Weinberg et al., 2006, 2008; Roosenboom et al., 2017). Our human data results suggested that PRICKLE1 may not operate by directly altering the width of the face, but instead via its control over how variable the phenotype can be. In the mouse cranial base, we observed the opposite, i.e., the Prickle1^Bj allele had a mean effect, but not a variance effect. This discrepancy may be explained by a combination of developmental timing, inherent species differences, and low power in detecting the variance effect in a small number of mice.

FOCAD encodes a focal adhesion complex protein (Focadhesin) (Brockschmidt et al., 2012). There is, so far, little known about its biological function, except that it is a potential tumor suppressor gene highly expressed in brain tissues (Weren et al., 2015). Intriguingly, two recent GWASes by our group in samples of different ancestries both reported genome-wide significant SNPs associated with facial morphology at 9p21.3, which encompasses FOCAD. One of the GWASes was performed in an African population, and the lead SNP at 9p21.3 was ∼357 kb upstream of FOCAD (Liu et al., 2021). The association signals were seen for the shape variations of the zygoma, nose bridge, and area the surrounding eyes, and the multidimensional representation of these facial regions likely captured some information about the width of the cranial base. The other GWAS, conducted among Europeans, found a significant association at 9p21.3 with the shape of the eye area (Hoskens et al., 2021). In addition to its association with facial variations in healthy populations, FOCAD was also suggested to be potentially involved in craniofacial malformation in a study where a parent-of-origin interaction effect between FOCAD and maternal smoking was reported for cleft lip with or without cleft palate (Haaland et al., 2019). Our findings about FOCAD and the existing research results together provided multiple lines of evidence converging on it being a potential player in craniofacial-related processes and traits.

Our observational results in humans were echoed, to some extent, in the mouse model. The Prickle1^Bj/Bj mice had significantly wider heads and abnormal localization of the Focadhesin protein. These data, together, strongly support the relevance of Prickle1 and Focad and their potential interaction in the development of the cranial base. Proximal distal growth of the cranial base is accomplished by the expansion of the paired synchondroses. Growth plate development involves a process where chondrocytes in a quiescent resting zone mature into proliferating chondrocytes that line up in columns that are oriented in the proximal distal axis of the cranial base (Kronenberg, 2003). In the cranial base, the proliferating chondrocytes initially divide in the medial–lateral axis and then slide into the proximal–distal columns through a process that requires cell adhesion and Wnt/PCP signaling (Romereim et al., 2014). Our observation of the co-localization of Prickle1 and Focadhesin in the cleavage furrow of normal proliferating chondrocytes supports the hypothesis that Prickle1 and Focadhesin determine the angle of cell division. The Prickle1^Bj allele affects protein function (Gibbs et al., 2016; Wan et al., 2018). Furthermore, our immunofluorescence data suggest that Prickle1 and Focadhesin may interact in this process. Future studies are needed to further elucidate the mechanism by which the two proteins operate together during craniofacial development.

In humans, despite the lack of variance effect of PRICKLE1 in the replication cohort, we achieved a locus-level replication of its interaction with FOCAD. The failure of our SNP-level replication was not surprising. Considering the two cohorts of different ancestries, differences in the LD structure, the non-availability of the discovery lead SNPs, and the differences in the genetic and environmental backgrounds may all have a role in shaping the distribution of the variance and the interaction effects and/or our ability to detect them. These factors may also explain why the interaction signals were detected at different locations relative to PRICKLE1 and FOCAD in the two cohorts—introns of the genes in discovery and outside of the genes in replication. The statistical signals in the two cohorts may implicate the same underlying biology, or they may in fact reflect different mechanisms through which genetic variants influence a facial feature. There was some evidence from existing epigenomic ChIP-Seq data that the non-coding regions with which the lead interacting SNPs overlapped may function as enhancers in relevant cell types (Supplementary Figure 4). Several major links need to be further established, including how the SNPs identified disrupt those enhancers, whether those enhancers target PRICKLE1 and FOCAD, and how PRICKLE1 and FOCAD together affect early craniofacial development. Despite these unknowns, the results from our mouse experiments provided a good reason to hypothesize that, even if different regulatory elements were implicated in the discovery and the replication analysis, these elements likely all targeted PRICKLE1 and FOCAD.

Although the statistical results of our human data that initially led to the experimental interrogation of Prickle1 and Focad were not decisively strong, the fact that we found strong evidence for a genuine biological interaction indicates that the lack of significance was likely due to the limited power. It is well recognized that detecting differences in group variances demands larger sample sizes than when detecting differences in group means. Given that a previous GWAS in the same cohort identified only a limited number of significant loci (Shaffer et al., 2016), we did not expect to see more than a few significant vQTLs. Indeed, the most significant vQTLs in our analysis were only at the suggestive level (5 × 10^–7). It is a limitation that we followed up on a strictly speaking non-significant signal (PRICKLE1), and in the next step of the analysis, its interaction with FOCAD again did not reach the genome-wide significance threshold. However, three things together motivated us to pursue a replication and mouse experiments. Firstly, the suggestive p-values of a true positive might be a mere result of insufficient power. Secondly, the disruption of PRICKLE1 was known to associate with craniofacial phenotypes, but nearby SNPs were never significant in the GWAS, implying a role of a rare mutation or a complicated mode of action, such as interaction. Thirdly, the genome-wide search of loci interacting with FOCAD did identify PRICKLE1 as the most and only significant locus (p = 1.88 × 10^–10).

The strategy of variance prioritization showed a clear benefit in studying G × G interactions in the present study. With a sample size of 2,447 individuals, an exhaustive search for pairwise interactions between genome-wide SNPs would not have sufficient power given the huge burden of multiple comparisons, even with the scope restricted to SNPs with a MAF >0.2. As PRICKLE1 and FOCAD were not known to relate in any way before this study, few other pre-selection strategies would have been able to sort them out and make a discovery. Note that the significant rs10880322 (PRICKLE1) × rs10511683 (FOCAD) interaction exemplified how interaction effects need not cause variance heterogeneity at either locus. Likewise, variance heterogeneity can arise from genetic processes other than interaction, as well as statistical artifacts such as outliers. We carefully scrutinized alternative explanations for the PRICKLE1 vQTL (see section “Results”), ensuring that, although not impossible, they were not likely to explain our finding. Despite the lack of definite relationship between heterogeneous variance and interaction, a detected vQTL does signify a possible presence of unmodeled statistical interaction, which can be further examined relatively easily. Our results highlighted the potential of vQTL analysis in revealing G × G interactions and underscored the need to explore how genetic factors crosstalk in facial genetic studies.

Data Availability Statement

Publicly available datasets were generated for this study. These data can be found here: the genetic data and the 20 facial distances evaluated in the United States cohort are available to the research community through the dbGaP controlled-access repository (https://dbgap.ncbi.nlm.nih.gov/) at accession phs000949.v1.p1. The raw source data for the phenotypes – the 3D facial surfaces – are available for the 3D Facial Norms dataset through the FaceBase Consortium (www.facebase.org). The genotypic and phenotypic data of the Korean sample is available upon request from the Korean Genome and Epidemiology Study (KoGES) of the Center for Disease Control and from the Korea Medicine Data Center of the Korea Institute of Oriental Medicine, respectively. The mouse data is available upon request from HSR.

Ethics Statement

The studies involving human participants were reviewed and approved by University of Pittsburgh Institutional Review Board #PRO09060553 and #RB0405013; UT Health Committee for the Protection of Human Subjects #HSC-DB-09-0508; Seattle Children’s Institutional Review Board #12107; University of Iowa Human Subjects Office/Institutional Review Board #200912764 and #200710721; and the Korea Institute of Oriental Medicine #I-2007/006-002. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin. The animal study was reviewed and approved by University of Pittsburgh Institutional Animal Care and Use Committee (#17050839).

Author Contributions

DL, SW, and JS conceived and designed the study. EF, MM, SW, and JS provided funding for data collection. JH, GW, and LM provided the discovery human data. H-JB and SC provided the replication human data and involved in the replication with Korean subjects. DL spearheaded the statistical analysis and wrote the manuscript. H-JB and ML contributed to the statistical analysis. AE and HS-R conducted the mouse experiments. All authors reviewed and approved the final manuscript.

Funding

This study was supported by the National Institute for Dental and Craniofacial Research (NIDCR; http://www.nidcr.nih.gov/) through the following grants: U01-DE020078 to SW and MM, R01-DE016148 to SW and MM, R01-DE027023 to SW and JS, and X01-HG007821 to MM, SW, and EF. Funding for the initial genomic data cleaning of the United States cohort by the University of Washington was provided by contract #HHSN268201200008I from the NIDCR awarded to the Center for Inherited Disease Research (CIDR; https://www.cidr.jhmi.edu/). The study that initially collected the Korean data was supported by the National Research Foundation of Korea (NRF-2018M3E3A1057298). The mouse work was supported by Momentum Funds from the University of Pittsburgh, P30-NIH-DE020740, and start-up funds from the University of Pittsburgh School of Dental Medicine (HLSR).

Conflict of Interest

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank all the dedicated staff, collaborators, and participants for their contribution to the study. The genotype data of the Korean sample were provided with bioresources from the National Biobank of Korea, the Centers for Disease Control and Prevention, South Korea (4845-301, 4851-302, and 307).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.674642/full#supplementary-material

Footnotes

1. ^ http://cnsgenomics.com/software/osca

References

Adhikari, K., Fuentes-Guajardo, M., Quinto-Sánchez, M., Mendoza-Revilla, J., Camilo Chacón-Duque, J., 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahsan, K., Singh, N., Rocha, M., Huang, C., and Prince, V. E. (2019). Prickle1 is required for EMT and migration of zebrafish cranial neural crest. Dev. Biol. 448, 16–35. doi: 10.1016/j.ydbio.2019.01.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Blanco, R., Cifuentes, L., Maldonado, M. J., Rameau, X., and Muñoz, M. A. (1992). [Cleft lip and cleft palate: cephalometric characteristics of affected individuals, their relatives and a control population]. Rev. Med. Chil. 120, 13–19.

Google Scholar

Brockschmidt, A., Trost, D., Peterziel, H., Zimmermann, K., Ehrler, M., Grassmann, H., et al. (2012). KIAA1797/FOCAD encodes a novel focal adhesion protein with tumour suppressor function in gliomas. Brain 135, 1027–1041. doi: 10.1093/brain/aws045

PubMed Abstract | CrossRef Full Text | Google Scholar

Cha, S., Lim, J. E., Park, A. Y., Do, J.-H., Lee, S. W., Shin, C., et al. (2018). Identification of five novel genetic loci related to facial morphology by genome-wide association studies. BMC Genomics 19:481. doi: 10.1186/s12864-018-4865-4869

CrossRef Full Text | Google Scholar

Chatzistavrou, E., Ross, R. B., Tompson, B. D., and Johnston, M. C. (2004). Predisposing factors to formation of cleft lip and palate: inherited craniofacial skeletal morphology. Cleft. Palate. Craniofac. J. 41, 613–621. doi: 10.1597/03-090.1

CrossRef Full Text | Google Scholar

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-54

CrossRef Full Text | Google Scholar

Cole, J. B., Manyama, M., Kimwaga, E., Mathayo, J., Larson, J. R., Liberton, D. K., et al. (2016). Genomewide association study of African children identifies association of SCHIP1 and PDE8A with facial size and shape. PLoS Genet 12:e1006174. doi: 10.1371/journal.pgen.1006174

PubMed Abstract | CrossRef Full Text | Google Scholar

Crouch, D. J. M., Winney, B., Koppen, W. P., Christmas, W. J., Hutnik, K., Day, T., et al. (2018). Genetics of the human face: identification of large-effect single gene variants. Proc. Natl. Acad. Sci. U S A. 115, E676–E685. doi: 10.1073/pnas.1708207114

PubMed Abstract | CrossRef Full Text | Google Scholar

Ek, W. E., Rask-Andersen, M., Karlsson, T., Enroth, S., Gyllensten, U., and Johansson, Å (2018). Genetic variants influencing phenotypic variance heterogeneity. Hum. Mol. Genet. 27, 799–810. doi: 10.1093/hmg/ddx441

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedorov, A., Beichel, R., Kalpathy-Cramer, J., Finet, J., Fillion-Robin, J.-C., Pujol, S., et al. (2012). 3D slicer as an image computing platform for the quantitative imaging network. Magn. Reson. Imaging 30, 1323–1341. doi: 10.1016/j.mri.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Forsberg, S. K. G., and Carlborg, Ö (2017). On the relationship between epistasis and genetic variance heterogeneity. J. Exp. Bot. 68, 5431–5438. doi: 10.1093/jxb/erx283

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraser, F. C., and Pashayan, H. (1970). Relation of face shape to susceptibility to congenital cleft lip. a preliminary report. J. Med. Genet. 7, 112–117. doi: 10.1136/jmg.7.2.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibbs, B. C., Damerla, R. R., Vladar, E. K., Chatterjee, B., Wan, Y., Liu, X., et al. (2016). Prickle1 mutation causes planar cell polarity and directional cell migration defects associated with cardiac outflow tract anomalies and other structural birth defects. Biol. Open 5, 323–335. doi: 10.1242/bio.015750

PubMed Abstract | CrossRef Full Text | Google Scholar

Haaland, ØA., Romanowska, J., Gjerdevik, M., Lie, R. T., Gjessing, H. K., and Jugessur, A. (2019). A genome-wide scan of cleft lip triads identifies parent-of-origin interaction effects between ANK3 and maternal smoking, and between ARHGEF10 and alcohol consumption. F1000Res 8:960. doi: 10.12688/f1000research.19571.2

PubMed Abstract | CrossRef Full Text | Google Scholar

He, F., and Chen, Y. (2012). Wnt signaling in lip and palate development. Front. Oral. Biol. 16:81–90. doi: 10.1159/000337619

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoskens, H., Liu, D., Naqvi, S., Lee, M. K., Eller, R. J., Indencleef, K., et al. (2021). 3D facial phenotyping by biometric sibling matching used in contemporary genomic methodologies. PLoS Genet. 17:e1009528. doi: 10.1371/journal.pgen.1009528

PubMed Abstract | CrossRef Full Text | Google Scholar

Hulse, A. M., and Cai, J. J. (2013). Genetic variants contribute to gene expression variability in humans. Genetics 193, 95–108. doi: 10.1534/genetics.112.146779

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, R. A., Sotoudeh, R., and Conley, D. (2020). Polygenic scores for plasticity: a new tool for studying gene-environment interplay. bioRxiv [Preprint]. doi: 10.1101/2020.08.30.274530

CrossRef Full Text | Google Scholar

Klingenberg, C. P. (2011). MorphoJ: an integrated software package for geometric morphometrics. Mol. Ecol. Resour. 11, 353–357. doi: 10.1111/j.1755-0998.2010.02924.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kronenberg, H. M. (2003). Developmental regulation of the growth plate. Nature 423, 332–336. doi: 10.1038/nature01657

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M. K., Shaffer, J. R., Leslie, E. J., Orlova, E., Carlson, J. C., Feingold, E., et al. (2017). Genome-wide association study of facial morphology reveals novel associations with FREM1 and PARK2. PLoS One 12:e0176566. doi: 10.1371/journal.pone.0176566

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S., Wu, M. C., and Lin, X. (2012). Optimal tests for rare variant effects in sequencing association studies. Biostatistics 13, 762–775. doi: 10.1093/biostatistics/kxs014

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Leal, S. M. (2008). Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. Am. J. Hum. Genet. 83, 311–321. doi: 10.1016/j.ajhg.2008.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., and Ji, L. (2005). Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. Heredity (Edinb) 95, 221–227. doi: 10.1038/sj.hdy.6800717

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, B., Gloudemans, M. J., Rao, A. S., Ingelsson, E., and Montgomery, S. B. (2019). Abundant associations with gene expression complicate GWAS follow-up. Nat. Genet. 51, 768–769. doi: 10.1038/s41588-019-0404-400

CrossRef Full Text | Google Scholar

Liu, C., Lee, M. K., Naqvi, S., Hoskens, H., Liu, D., White, J. D., et al. (2021). Genome scans of facial features in east africans and cross-population comparisons reveal novel associations. PLoS Genet. (in press).

Google Scholar

Liu, C., Lin, C., Gao, C., May-Simera, H., Swaroop, A., and Li, T. (2014). Null and hypomorph Prickle1 alleles in mice phenocopy human robinow syndrome and disrupt signaling downstream of Wnt5a. Biol. Open 3, 861–870. doi: 10.1242/bio.20148375

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Marderstein, A. R., Davenport, E. R., Kulm, S., Van Hout, C. V., Elemento, O., and Clark, A. G. (2021). Leveraging phenotypic variability to identify genetic interactions in human phenotypes. Am. J. Hum. Genet. 108, 49–67. doi: 10.1016/j.ajhg.2020.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

McIntyre, G. T., and Mossey, P. A. (2004). Parental craniofacial morphology in orofacial clefting. Eur. J. Orthod. 26, 375–384. doi: 10.1093/ejo/26.4.375

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, J. H. (2005). A global view of epistasis. Nat. Pub. Group 37, 13–14. doi: 10.1038/ng0105-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, J. H., and Williams, S. M. (2005). Traversing the conceptual divide between biological and statistical epistasis: systems biology and a more modern synthesis. Bioessays 27, 637–646. doi: 10.1002/bies.20236

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, J. H., and Williams, S. M. (2009). Epistasis and its implications for personal genetics. Am. J. Hum. Genet. 85, 309–320. doi: 10.1016/j.ajhg.2009.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Mossey, P. A., McColl, J., and O’Hara, M. (1998). Cephalometric features in the parents of children with orofacial clefting. Br. J. Oral Maxillofac. Surg. 36, 202–212. doi: 10.1016/s0266-4356(98)90498-90493

CrossRef Full Text | Google Scholar

Mostowska, A., Hozyasz, K. K., Biedziak, B., Wojcicki, P., Lianeri, M., and Jagodzinski, P. P. (2012). Genotype and haplotype analysis of WNT genes in non-syndromic cleft lip with or without cleft palate. Eur. J. Oral Sci. 120, 1–8. doi: 10.1111/j.1600-0722.2011.00938.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakasima, A., and Ichinose, M. (1983). Characteristics of craniofacial structures of parents of children with cleft lip and/or palate. Am. J. Orthod. 84, 140–146. doi: 10.1016/0002-9416(83)90178-90171

CrossRef Full Text | Google Scholar

Paré, G., Cook, N. R., Ridker, P. M., and Chasman, D. I. (2010). On the use of variance per genotype as a tool to identify quantitative trait interaction effects: a report from the Women’s genome health study. PLoS Genet. 6:e1000981. doi: 10.1371/journal.pgen.1000981

PubMed Abstract | CrossRef Full Text | Google Scholar

Paternoster, L., Zhurov, A. I., Toma, A. M., Kemp, J. P., Pourcain, B. S., Timpson, N. J., 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. doi: 10.1016/j.ajhg.2011.12.021

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Rask-Andersen, M., Karlsson, T., Ek, W. E., and Johansson, Å (2017). Gene-environment interaction study for BMI reveals interactions between genetic factors and physical activity, alcohol consumption and socioeconomic status. PLoS Genet. 13:e1006977. doi: 10.1371/journal.pgen.1006977

PubMed Abstract | CrossRef Full Text | Google Scholar

Romereim, S. M., Conoan, N. H., Chen, B., and Dudley, A. T. (2014). A dynamic cell adhesion surface regulates tissue architecture in growth plate cartilage. Development 141, 2085–2095. doi: 10.1242/dev.105452

PubMed Abstract | CrossRef Full Text | Google Scholar

Rönnegård, L., and Valdar, W. (2012). Recent developments in statistical methods for detecting genetic loci affecting phenotypic variability. BMC Genet. 13:63. doi: 10.1186/1471-2156-13-63

PubMed Abstract | CrossRef Full Text | Google Scholar

Roosenboom, J., Indencleef, K., Hens, G., Peeters, H., Christensen, K., Marazita, M. L., et al. (2017). Testing the face shape hypothesis in twins discordant for nonsyndromic orofacial clefting. Am. J. Med. Genet. 173, 2886–2892. doi: 10.1002/ajmg.a.38471

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarkar, A. K., Tung, P.-Y., Blischak, J. D., Burnett, J. E., Li, Y. I., Stephens, M., et al. (2019). Discovery and characterization of variance QTLs in human induced pluripotent stem cells. PLoS Genet. 15:e1008045. doi: 10.1371/journal.pgen.1008045

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Struchalin, M. V., Dehghan, A., Witteman, J. C., van Duijn, C., and Aulchenko, Y. S. (2010). Variance heterogeneity analysis for detection of potentially interacting genetic loci: method and its limitations. BMC Genet. 11:92. doi: 10.1186/1471-2156-11-92

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, A., Takenoshita, Y., Honda, Y., and Matsuura, C. (1999). Dentocraniofacial morphology in parents of children with cleft lip and/or palate. Cleft. Palate. Craniofac. J. 36, 131–138.

Google Scholar

Tao, H., Suzuki, M., Kiyonari, H., Abe, T., Sasaoka, T., and Ueno, N. (2009). Mouse prickle1, the homolog of a PCP gene, is essential for epiblast apical-basal polarity. Proc. Natl. Acad. Sci. U S A. 106, 14426–14431. doi: 10.1073/pnas.0901332106

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, Y., Lantz, B., Cusack, B. J., and Szabo-Rogers, H. L. (2018). Prickle1 regulates differentiation of frontal bone osteoblasts. Sci. Rep. 8:18021. doi: 10.1038/s41598-018-36742-36740

CrossRef Full Text | Google Scholar

Wang, H., Zhang, F., Zeng, J., Wu, Y., Kemper, K. E., Xue, A., et al. (2019). Genotype-by-environment interactions inferred from genetic effects on phenotypic variability in the UK Biobank. Sci. Adv. 5:eaaw3538. doi: 10.1126/sciadv.aaw3538

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinberg, S. M., Maher, B. S., and Marazita, M. L. (2006). Parental craniofacial morphology in cleft lip with or without cleft palate as determined by cephalometry: a meta-analysis. Orthod. Craniofac. Res. 9, 18–30. doi: 10.1111/j.1601-6343.2006.00339.x

PubMed Abstract | CrossRef Full Text | Google Scholar

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. A 146A, 409–420. doi: 10.1002/ajmg.a.32177

PubMed Abstract | CrossRef Full Text | Google Scholar

Weren, R. D. A., Venkatachalam, R., Cazier, J.-B., Farin, H. F., Kets, C. M., de Voer, R. M., et al. (2015). Germline deletions in the tumour suppressor gene FOCAD are associated with polyposis and colorectal cancer development. J. Pathol. 236, 155–164. doi: 10.1002/path.4520

PubMed Abstract | CrossRef Full Text | Google Scholar

White, J. D., Indencleef, K., Naqvi, S., Eller, R. J., Hoskens, H., Roosenboom, J., et al. (2020). Insights into the genetic architecture of the human face. Nat. Genet. 66, 101–109. doi: 10.1038/s41588-020-00741-747

CrossRef Full Text | Google Scholar

Yang, T., Bassuk, A. G., and Fritzsch, B. (2013). Prickle1 stunts limb growth through alteration of cell polarity and gene expression. Dev. Dyn. 242, 1293–1306. doi: 10.1002/dvdy.24025

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, T., Jia, Z., Bryant-Pike, W., Chandrasekhar, A., Murray, J. C., Fritzsch, B., et al. (2014). Analysis of PRICKLE1 in human cleft palate and mouse development demonstrates rare and common variants involved in human malformations. Mol. Genet. Genomic Med. 2, 138–151. doi: 10.1002/mgg3.53

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoon, Y.-J., Perkiomaki, M. R., Tallents, R. H., Barillas, I., Herrera-Guido, R., Fong, C.-T., et al. (2004). Transverse craniofacial features and their genetic predisposition in families with nonsyndromic unilateral cleft lip and palate. Cleft. Palate. Craniofac. J. 41, 256–261. doi: 10.1597/02-134.1

CrossRef Full Text | Google Scholar

Young, A. I., Wauthier, F. L., and Donnelly, P. (2018). Identifying loci affecting trait variability and detecting interactions in genome-wide association studies. Nat. Genet. 50, 1608–1614.

Google Scholar

Zuk, O., Hechter, E., Sunyaev, S. R., and Lander, E. S. (2012). The mystery of missing heritability: genetic interactions create phantom heritability. Proc. Natl. Acad. Sci. U S A. 109, 1193–1198. doi: 10.1073/pnas.1119675109

PubMed Abstract | CrossRef Full Text | Google Scholar