Craniofacial anomalies (CFAs) are a common group of congenital anomalies (CAs) caused by the abnormal development of skull and/or facial bones. The most common CFAs, cleft lip with or without cleft palate (CL/P) and cleft palate alone, make up one third of all CAs in the United States. They occur in ~ 16 in 10,00 births within the United States and are a major contributor to infant mortality.¹ Likewise, craniosynostosis, another common CFA, is estimated to affect 1 in 2100–2500 births and can result in abnormal brain growth causing neurological dysfunction.² The frequency, contribution to infant mortality, and life-long associated health problems associated with CFAs makes understanding of genetic underpinnings of the risk and presentations of CFAs essential.

Despite a great deal of research into the genetic etiology of CFAs, clinical genetic testing for CFAs has a relatively low diagnostic yield. There is little understanding of what drives the variable expression and incomplete penetrance of CFA syndromes. Combined diagnostic approaches of karyotype, chromosomal microarray, and exome sequencing can detect the genetic cause of 22.5% of individuals with orofacial clefts.^3,4 Similarly, comprehensive clinical diagnostic genetic testing of individuals with craniosynostosis can detect a genetic cause in about a quarter (25%) of affected individuals.⁵ These results suggest that, at best, for every four patients with a CFA who access comprehensive genetic testing, one will receive a diagnostic results for the genetic cause of their condition. Even when a diagnostic test identifies a causal variant in a Mendelian gene, there is often incomplete penetrance and/or variable expressivity, suggesting the presence of complex genomic and environmental interactions that impact phenotypic manifestations. As sequencing technologies have advanced and sample sizes have increased, several studies have examined how polygenic variation impacts the penetrance of rare variants associated with CFAs such as cleft lip and palate.^6,7

To better understand how polygenicity impacts CFAs, large sample sizes are essential. Resources such as Vanderbilt University Medical Center’s BioVU and the Electronic Medical Records and Genomics (eMERGE) network cohort, which link genotype data to electronic health record (EHR) data, could contribute to these efforts.^8,9 We therefore undertook the first use of these resources to evaluate the genetically predicted changes in gene expression of CFA-associated Mendelian genes as well as the genetically predicted gene expression of individuals with CFAs.

Vanderbilt University Medical Center’s EHR-system provides a unique resource with a wealth of de-identified health information in the synthetic derivative (SD). The SD contains ~ 3.5 million individuals, with ~ 300k individuals having genetic information in BioVU, VUMC’s deidentified EHR-linked biobank. 6,131 Individuals with CAs and CFAs were identified from Vanderbilt University Medical Center (VUMC)’s de-identified EHR-linked DNA biobank, BioVU.^8,12 This work is deemed non-human subjects work by the VUMC IRB and received all necessary approvals. We restricted our cases and controls used in analyses to a medical home population (n = 1,275,576), a high medical use population defined as individuals with an ICD-9 or 10 billing code collected from at least three unique visit dates over three or more years. This medical home definition ensures a study population with substantial phenotypic information within their EHR. All genetic data for included participants in BioVU were genotyped using the Illumina MEGAEX Array and the included population was restricted to participants of European genetic ancestry based on clustering in principal component analysis (PCA) using genetic data from the 1000 Genomes Project as reference populations, as previously described.¹³ This restriction was done to minimize confounding results and maximize sample size. Future work to include individuals of non-European genetic ancestry is ongoing.

We identified individuals with CFAs from the Electronic Medical Records and Genomics (eMERGE) Network, which combines DNA biorepositories with EHR data.^9,14 The eMERGE GWAS cohort includes 64,536 individuals from five institutions (VUMC, Columbia University Irving Medical Center, Northwestern Medical Center, Mass General Brigham, and Cincinnati Children’s Hospital Medical Center (CCHMC)). Analyses were restricted to participants of European genetic ancestry from CCHMC in the eMERGE GWAS cohort. Because 75% of CFA cases in eMERGE were from CCHMC we restricted our secondary analysis to this population. ICD billing codes for each individual are provided as part of available data to eMERGE consortium researchers. We mapped the CM phecodes from the ICD code data using the PheWAS R package as previously described.¹¹

We curated a list of genes associated with CFAs, using previously published genome-wide association study (GWAS) results and the Concert Genetics’ registry of clinical genetic tests.^15,16 This testing registry aggregates all clinically available genetic panels, lists all the genes included in each panel, and allows for comparison of included genes between different companies’ panels. These panels reflect a largely comprehensive list of genes with a known monogenic CFA association. A certified genetic counselor searched Concert Genetics’ Test Registry using the search term “craniofacial” and compiled a list of all the unique genes that are offered on the resulting Craniofacial Panel Tests. These panels included both syndromic and non-syndromic CFA genes. Supplemental Table 1 lists curated known CFA-associated genes.

To identify individuals in BioVU and eMERGE that have a CFA, we used phecodes from the “Tongue Mouth and Pharynx” and the “Skull Face and Jaw” phecodeX Chaps. 1⁰ We identified individuals that had at least two instances of a CFA phecode of the same parent code and/or family head code Supplemental Table 2. Figure 1B demonstrates the breakdown of the parent code, family head code, and specific code portions of a phecode. We defined controls as individuals who did not have a single congenital anomaly phecode in their record.

Using machine learning models such as PrediXcan, UTMOST, and Joint Tissue Imputation (JTI), we calculated genetically predicted gene expression (GPGE) using GTEx (version 8) as a reference population in both BioVU and eMERGE.^17,18,19,20 Using the single best performing model out of these three models for each gene tissue pair (r² > 0.01), we conducted a transcriptome-wide association study (TWAS). We tested the association of CFA status with GPGE in a logistic regression model adjusting for age, sex, number of visits, and the first ten principal components of ancestry.

To determine how variation in GPGE of known CFA genes may increase risk for other CAs, we conducted gene-based phenotype-wide association studies (PheWAS) of the known CFA genes.^21,22,23 Of the 391 known CFA genes, we had quality prediction for 341 (r² > 0.1). We tested whether a diagnosis of a CA is associated with GPGE of a known CFA gene adjusting for age, sex, the first ten principal components of ancestry, and number of visits in a logistic regression mode.

Because most previous genetic studies of CFAs focus on rare variant/ monogenic causes of disease, this work investigated the GPGE of individuals with CFAs and how the GPGE of known CFA genes relate to CAs. We perform a TWAS for individuals with a CFA and a CA-wide PheWAS for the GPGE of genes with a known CFA association. Overall, the results of these analyses suggest that in addition to rare variants, polygenic variation impacting gene expression may contribute to many CAs and may play a role in the penetrance and expressivity of CA syndromes.

As sample sizes for genetic studies have increased, so has our understanding of the complexity of the genetic architecture driving phenotypes such as CFAs through our expanded understanding of how common variation contributes to the diversity of facial morphology and shape.^15,26,27 Studies of CFAs and polygenic architecture have identified that there are shared genetic features between cleft lip and palate and the size of facial features as well as common variation that affects the penetrance of known cleft clip and palate variants in PDGFRA (*173490).^6,7 Additionally, researchers have proposed that Mendelian CFAs are extreme phenotypes on a continuum of phenotypic variation in facial morphology and that integrating common variation into the study of these phenotypes is essential to understanding their genetic drivers.^27,28 All of these studies leverage GWAS and polygenic scores to examine common variations contributions to CFAs. However, using an approach that utilizes GPGE allows us to capture the predicted effects of common variant gene expression in a TWAS to identify genes whose altered GPGE is associated with CFA status. By conducting a gene-based analysis of CFAs, we can obtain more biologically interpretable results for our CFA associations.

Out of the 391 genes known to cause CFAs in a monogenic fashion, over 90% are not significantly associated with CFAs in a transcriptome-wide fashion (p < 0.01), highlighting the gap in knowledge when studying CFAs that present in a Mendelian or syndromic fashion. The two genes whose GPGE showed a significant association (p < 0.01) with CFAs in both cohorts, VAV1 (*164875) and CYP3A7 (*605340), have not yet been associated with any human disease, but do have known roles in fetal development. VAV1 (*164875) is a proto-oncogene that is involved with hematopoiesis and T and B cell signaling. The well-established relationship between disrupted oncoprotein signaling, cancer, and congenital anomalies suggests that differential expression of the oncogene VAV1 (*164875) could be driving the development of some CFAs.^{29, 30,24} Additional research on the differential expression of proto-oncogenes in individuals with congenital anomalies could give additional insight into why those with congenital anomalies are more likely to develop cancer. CYP3A7 (*605340) encodes a cytochrome P450 “super family” enzyme, CYP3A7 (*605340), which metabolizes a diverse array of endogenous and exogenous substances, including prescription medications that many pregnant individuals need to maintain their own health such as carbamazepine, diltiazem, caffeine, and nifedipine. CYP3A7 (*605340) is primarily expressed in fetal liver tissue, being detected as early as 50 days of gestation and decreasing in expression until 24 months of age postnatally.^31,32 Additionally, CYP3A7 (*605340) plays a key role in the production of a critical pregnancy hormone, estriol, which has been shown to be an important epigenetic modifier in mice fetuses. Variable expression of CYP3A7 (*605340) could have dramatic effects on fetal development, and further research can assess the complex interactions of environmental risks and genetic predispositions to CFAs and other CAs.

The observations that only 10% of the previously identified CFA associated genes had a significant GPGE association with CFAs and the strong overlap (21.1%) of individuals with CFAs having other CAs drove us to analyze what CAs were significantly associated with the GPGE of the curated CFA genes. These known CFA genes show more significant GPGE associations with other organ system CAs than the CFAs themselves, suggesting that there are shared genetic and environmental susceptibilities across CAs. The many genetic syndromes that contain multiple CAs support the idea of shared risk factors for many types of CAs. Taken together these results suggest that polygenic variation in CFA-associated genes may relate to developmental changes more broadly and are not necessarily restricted to craniofacial development.

Many genetic syndromes that cause multiple types of CAs demonstrate variable expression and incomplete penetrance, yet the factors causing variable expression/ incomplete penetrance is not well understood. For example, ~ 65% of individuals with 22q11.2 deletion syndrome have a congenital heart defect (CHD) and ~ 67% have a palate abnormality.^33,34 As another example, we noted that the predicted expression of GLI2’s (*165230) is significantly (p < 0.05) associated with many congenital heart defects (CHDs), such as congenital pulmonary valve stenosis (p = 0.0002), congenital malformations of heart valves (p = 0.0002), congenital insufficiency of the aortic valve (p = 0.0005), and ten others (Table 3). GLI2 (*165230) is classified as a known CFA gene due to its association with two congenital malformation syndromes, Culler-Jones syndrome and Holoprosencephaly 9.³⁵ Both syndromes can present with several congenital anomalies, such as cleft lip/palate, microcephaly, polydactyly, but to date CHDs are not associated with either syndrome.^36,37 Yet GLI2 (*165230) does have a well-established role in cardiomyogenesis and there is a group of individuals with CHDs that have GLI2 (*165230) missense variants shown to dysregulate sonic hedgehog signaling, which is crucial for fetal development.^,38,39,40 The expression of GLI2 (*165230) in the developing heart suggests CHDs could possibly be a phenotypic expansion for the two GLI2-related congenital anomaly syndromes and suggests examining its role in cardiac development. Differences in gene expression could be one factor causing variable expressivity and reduced penetrance that is characteristic of many congenital anomaly syndromes.

One of the main limitations of studying CFAs at biobank scale is that they have a relatively low prevalence and are caused by large-effect rare variants. Both issues affect our statistical power to detect phenotype associations. One way that our analysis attempted to address this limitation is by conducting gene-based analyses and defining our CFA phenotype across multiple phecodes that describe different types of CFAs. Despite using a less stringent p-value threshold, conducting gene-based analyses such as TWAS provides a more interpretable biological unit than single variant analyses. While the analyses in this paper are underpowered, they still provide biologically and clinically meaningful results.^22,23

Throughout this study we leveraged a set of known CFA genes that were compiled from clinical diagnostic testing for CFAs as well as a GWAS of face shape that was curated and reviewed by a certified genetic counselor. While we tried to make this gene set as comprehensive as possible, we are limited by which genes are currently have a well-established associated with CFAs. The main goal of the work in this study is to try to better understand the genetic drivers of CFAs.

Overall, our results support that both rare and common genetic variants in CFA Mendelian genes may contribute to a variety of CAs and highlights the complexities of the CA phenotypes, suggesting there are shared underlying genetic and environmental risk factors. Further research of CAs through GPGE could help better explain variable presentations and penetrance of CA syndromes.

TWAS of individuals with CFAs in the BioVU and eMERGE cohorts identified relatively few previously identified CFA-associated genes. The two genes whose GPGE had the strongest association in both cohorts have potential roles in the complex genetic drive of CAs. The predicted expression of genes that have a known Mendelian-association with CFAs are more often significantly associated with other types of CAs in the BioVU and eMERGE cohorts. For example, the predicted expression of GLI2, which is associated with a syndrome that can include CFAs, is significantly associated with several CHDs. The results of both analyses suggest that there are overlapping polygenic causes of many types of CAs and that with further research may help explain the variability in how CA syndromes can present.