A Genome-Wide Search for Gene-Environment Effects in Isolated Cleft Lip with or without Cleft Palate Triads Points to an Interaction between Maternal Periconceptional Vitamin Use and Variants in ESRRG

Background: It is widely accepted that cleft lip with or without cleft palate (CL/P) results from the complex interplay between multiple genetic and environmental factors. However, a robust investigation of these gene-environment (GxE) interactions at a genome-wide level is still lacking for isolated CL/P. Materials and Methods: We used our R-package Haplin to perform a genome-wide search for GxE effects in isolated CL/P. From a previously published GWAS, genotypes and information on maternal periconceptional cigarette smoking, alcohol intake, and vitamin use were available on 1908 isolated CL/P triads of predominantly European or Asian ancestry. A GxE effect is present if the relative risk estimates for gene-effects in the offspring are different across exposure strata. We tested this using the relative risk ratio (RRR). Besides analyzing all ethnicities combined (“pooled analysis”), separate analyses were conducted on Europeans and Asians to investigate ethnicity-specific effects. To control for multiple testing, q-values were calculated from the p-values. Results: We identified significant GxVitamin interactions with three SNPs in “Estrogen-related receptor gamma” (ESRRG) in the pooled analysis. The RRRs (95% confidence intervals) were 0.56 (0.45–0.69) with rs1339221 (q = 0.011), 0.57 (0.46–0.70) with rs11117745 (q = 0.011), and 0.62 (0.50–0.76) with rs2099557 (q = 0.037). The associations were stronger when these SNPs were analyzed as haplotypes composed of two-SNP and three-SNP combinations. The strongest effect was with the “t-t-t” haplotype of the rs1339221-rs11117745-rs2099557 combination [RRR = 0.50 (0.40–0.64)], suggesting that the effects observed with the other SNP combinations, including those in the single-SNP analyses, were mainly driven by this haplotype. Although there were potential GxVitamin effects with rs17734557 and rs1316471 and GxAlcohol effects with rs9653456 and rs921876 in the European sample, respectively, none of the SNPs was located in or near genes with strong links to orofacial clefts. GxAlcohol and GxSmoke effects were not assessed in the Asian sample because of a lack of observations for these exposures. Discussion/Conclusion: We identified significant interactions between vitamin use and variants in ESRRG in the pooled analysis. These GxE effects are novel and warrant further investigations to elucidate their roles in orofacial clefting. If validated, they could provide prospects for exploring the impact of estrogens and vitamins on clefting, with potential translational applications.


INTRODUCTION
With a prevalence of 3.4-22.9 per 10,000 live births (Mossey and Castilla, 2003), cleft lip with or without cleft palate (CL/P) ranks among the most common birth defects in humans. It is widely accepted that CL/P results from the complex interplay between multiple genetic and environmental factors (GxE), but only recently have data and practical approaches become available to enable an investigation of these effects at a genomewide level. Identifying GxE interactions may not only provide important insights into the etiology of orofacial clefts, but may also be important from a public health perspective because of the potential for interventions on environmental risk factors alone, particularly in genetically more susceptible subgroups of the population. This rationale has long been demonstrated in animal models (Millicovsky and Johnston, 1981a,b;Juriloff, 2002), but the evidence in humans is less conclusive. Just as some murine strains are more susceptible to external teratogens, human fetuses harboring high-risk alleles may also be more sensitive to particular environmental risk factors, and identifying those that interact with genetic risk variants may lead to important inroads in our understanding of the causes of CL/P.
Over the years, a growing list of maternal exposures has been reported to influence the risk of isolated CL/P. In particular, cigarette smoking (Zeiger and Beaty, 2002;Little et al., 2004;Lie et al., 2008), alcohol consumption (DeRoo et al., 2008), folic acid and other B-complex vitamin supplementation (Hayes, 2002;Munger, 2002;Munger et al., 2004;Wilcox et al., 2007), and anti-folate medication (Hernández-Díaz et al., 2000;Holmes et al., 2001) have been among the most widely studied exposures. Associations with other environmental risk factors have also been reported (reviewed in Dixon et al., 2011), but these have been less consistent across studies and no consensus has yet emerged on the harmful effects of these exposures.
A wide variety of study designs have been used to study GxE effects in orofacial clefts, one of which involves treating a caseparent triad as the unit of analysis. This family-based "triad design" handles bias due to population stratification by using non-transmitted parental alleles as controls, to be compared with the alleles transmitted to the case child (Gjessing and Lie, 2006). We and others have used the triad design to investigate GxE effects in orofacial clefts (Jugessur et al., 2003;Shi et al., 2007;Wu et al., 2010). In terms of statistical power, a unit comprising one case and one control provides approximately the same power as a complete triad when studying single-SNP associations (Schaid, 2002). For haplotype reconstruction, however, the triad design offers the additional advantage of allowing haplotypes to be deduced from the family structure (Gjessing and Lie, 2006).
A handful of genome-wide association studies (GWAS) have now been published on orofacial clefts (reviewed in Dixon et al., 2011;Mangold et al., 2011;Rahimov et al., 2012;Beaty et al., 2016), offering unprecedented opportunities for investigating GxE effects at a genome-wide level. Using data on case-parent triads from a previously published GWAS (Beaty et al., 2010), Wu et al. (2014) and Beaty et al. (2011) screened for GxE effects in one category of orofacial clefts: isolated cleft palate only (CPO). Here, we used the same GWAS dataset to search for GxE effects in the larger category of isolated CL/P.

Study Participants
Participants in this study stem from an international cleft collaboration involving seven Asian and six European/US populations. Characteristics of the study populations and details of the GWAS have been provided in Beaty et al. (2010). Briefly, genotyping was performed on an Illumina Human610-Quad R platform and genotypes for 589,945 SNPS were deposited in the Database of Genotypes and Phenotypes (dbGaP; http://www.ncbi.nlm.nih.gov/gap) under study accession ID phs000094.v1.p1. Besides genotypes, information was also available on maternal cigarette smoking, alcohol intake, and vitamin use in the periconceptional period (i.e., 3 months prior to conception through the first trimester of pregnancy). Interviews and questionnaires were used to assess maternal exposures and the data were coded as simple yes/no responses for cigarette smoking, any reported alcohol consumption, and any use of vitamin supplements. For a detailed description of these maternal exposures, see the recent work on isolated CPO by Wu et al. (2014) and the study outline at dbGAP under study accession ID phs000094.v1.p1. Quality control was performed using the approach outlined in Haaland et al. (2017). The number of SNPs before and after the pruning process is provided in Table 1. Figure 1 shows the distribution of triads by ethnicity and environmental exposure after quality control was performed. Note that we did not analyze GxSmoke or GxAlcohol in the Asian sample because very few of the mothers reported smoking cigarettes or drinking alcohol in the periconceptional period. The FIGURE 1 | Distribution of families according to ethnicity and maternal exposure status. The gray box shows the pooled sample consisting of all ethnicities (N = 1,908). White boxes show the number of families according to ethnicity and exposure (yes: exposed; no: unexposed; NA: data not available). Note that we did not consider GxSmoke or GxAlcohol analyses in the Asian sample because of a very low number of observations for these exposures.
final version of the dataset included genotypes for 1,908 nuclear families; 825 were of European ancestry, 1,024 of Asian descent, and 59 of other ethnicities. Table 2 provides an overview of triad completeness and ethnicity. Among the 1,908 nuclear families, 1,594 were complete mother-father-child triads and 314 were parent-child dyads only.

Statistical Analysis
We are continuously extending a comprehensive R-package, Haplin, for analyzing different constellations of family-based data (Gjessing and Lie, 2006). Haplin uses maximum likelihood to estimate the relative risk (RR), p-value, and 95% confidence intervals for a given risk allele or haplotype. To identify GxE effects, a Wald test is used to test whether RR estimates for geneeffects in the offspring are significantly different across strata of exposed and unexposed mothers, as described in our previous works (Skare et al., 2012;Gjerdevik et al., 2017). A GxE effect is assessed as the ratio of relative risks (RRR). That is, RRR = RR(1)/RR(0), where RR(1) and RR(0) are the relative risks for the offspring of exposed and unexposed mothers, respectively. If RRR = 1, the RR among the exposed is different from the RR among the unexposed, which would indicate a GxE effect. See Haaland et al. (2017) for more details on the calculation of RRR. Since its inception, Haplin has incorporated haplotype estimation for triads. By using the expectation-maximization (EM) algorithm (Dempster et al., 1977), Haplin reconstructs haplotypes from the multi-SNP data even though phase is not known from the observed markers alone. The EM-algorithm also enables Haplin to account for missing parental genotypes, thus allowing the 314 parent-child dyads in the current dataset to be included in the analyses.
If an offspring was homozygous for the risk allele at a locus, we assumed a multiplicative dose-response model. In other words, the relative risk with two copies of the risk allele was assumed to be simply RR × RR. For each maternal exposure (smoking, alcohol use, vitamin use), we performed separate analyses on Asians (n = 1,024), Europeans (n = 825), and all participants combined irrespective of their ethnicity ("pooled sample"; N = 1,908). The different sets of analyses are illustrated in Figure 1.
Further, we conducted haplotype analyses for SNPs that had very low q-values. This was done using the haplinSlide function in Haplin, as described in our previous works on candidate genes for clefts (Jugessur et al., 2008(Jugessur et al., , 2009b(Jugessur et al., , 2012a. Next, we performed statistical power calculations using the hapPowerAsymp function in Haplin. It is based on the asymptotic variance-covariance structure of the Wald test, as described in our recent work . The sample sizes used in these calculations reflect those that were available in the current GWAS dataset (Figure 1).
To facilitate the interpretation of the large number of statistical tests being performed simultaneously, we generated Quantile-Quantile (QQ) plots (Wilk and Gnanadesikan, 1968) using the pQQ function in Haplin to screen for significant p-values. This was done separately for the analysis of each exposure. To control the false discovery rate, we used the method described by Storey and Tibshirani (2003) where q-values are computed from the p-values. The function qvalue from the Rpackage "q-value" in the Bioconductor repository [https://www. bioconductor.org/, Gentleman et al., 2004] was used to generate the q-values. A low q-value indicates a low likelihood that a statistically significant association is a false positive. A q-value of 0.1, for instance, corresponds to a false discovery rate of 10%.
A regional plot was generated to map the chromosomal area flanking the most significant SNP ("lead SNP"). This plot provides a visual display of linkage disequilibrium between the lead SNP and other SNPs, along with their positions. It also provides the position of the closest genes and the recombination rates in the region. We generated the plot by modifying an R-script available at http://www.broadinstitute.org/files/shared/ diabetes/scandinavs/assocplot.R (Pruim et al., 2010).
The plethora of publicly available data makes it feasible to explore the connections between different entities, e.g., genes and vitamins. Hetionet integrates different sources of information from actively maintained meta-databases and allows the visualization of possible paths between various diseases, genes and compounds (Himmelstein et al., 2017, https://neo4j. het.io/). For example, to explore possible links between a given disease, compounds and gene nodes, Hetionet uses the following databases: Disease Ontology (Kibbe et al., 2015), Drug Bank (Law et al., 2014), and Entrez Gene (Maglott et al., 2005). Furthermore, Hetionet can also integrate information from the Human Interactome Database (Rolland et al., 2014), CMap (https://clue.io/) or Bgee (Bastian et al., 2008) to identify potential gene-gene interactions, compound-gene interactions and gene expression in specific tissue, respectively. We queried Hetionet to visualize the relationship between genes and maternal exposures for the cleft lip phenotype. The complete query and output codes are provided in the online Supplementary Text 1.

Ethics Approvals
Ethics approvals for the International Cleft Consortium were obtained from the respective institutional review boards of the participating sites. The consortium was formed in 2007 and each participating institution approved research protocols for the recruitment of case-parent triads from 13 individual sites. All participants have granted their written informed consents.

RESULTS
Tables 3-5 provide the relative risk ratios (RRRs) and 95% confidence intervals (CIs) for the top 20 SNPs in each set of analyses depicted in Figure 1. These SNPs were ranked by p-value, and their corresponding q-values are also provided in these tables. RRRs were calculated as previously described by Haaland et al. (2017). For SNPs that are not associated with isolated CL/P, -log 10 (p-values) in the QQ-plots are expected to fall along the straight diagonal line representing the null distribution, within the 95% confidence band (gray area). Conversely, -log 10 (p-values) would fall above this line, outside the confidence band and in the upper right-hand corner of the plot. Several SNPs were significantly associated with isolated CL/P in the GxVitamin analyses of the pooled sample ( Figure 2). In the other analyses, none of the SNPs were outside the confidence bands, suggesting that the distribution of pvalues is as expected if there are no associations. Still, in each of the GxVitamin and GxAlcohol analyses of the European sample (Tables 3, 5), two SNPs had markedly lower p-values than expected (see also Figures 5, 7). However, none of these SNPs was located in or near genes with obvious links to orofacial clefts. In the remainder of this section, therefore, we focus mainly on the results of the GxVitamin analysis of the pooled sample.
All the top 20 SNPs in the GxVitamin analysis of the pooled sample had q-values below 0.54 ( Table 3). The two top SNPs among these, rs1339221 and rs11117745, both had a qvalue of 0.01 and are located in the gene for "Estrogen-related receptor gamma" (ESRRG). These SNPs are in strong linkage disequilibrium (LD) with one another (R 2 = 0.77 and D' = 1.00). A third SNP, rs2099557 (q = 0.037), is also located in ESRRG but is in weaker LD with the above two SNPs (R 2 =0.33 and D' = 0.57 with rs1339221; R 2 = 0.28 and D' = 0.61 with rs11117745). A regional plot centered around the lead SNP in ESRRG (rs1339221) was generated to visualize the strength of the association signals and regional information around that SNP (Figure 8). Next, we performed stratified analyses of the SNPs in ESRRG by first testing for an overall child effect in the unstratified sample, followed by the effects among children who were exposed and unexposed to maternal vitamin use, respectively ( Table 6). All three SNPs exhibited a so-called "qualitative interaction" (Clayton, 2009), in that the effect of the SNP-allele among mothers taking vitamins was  in the opposite direction of that among mothers not taking vitamins. Specifically, RR>1 among non-takers, whereas RR<1 among takers, and none of the 95% CIs included 1 ( Table 6).
There were no statistically significant overall effects of the child's allele alone for any of the SNPs; all the 95% CIs included RR = 1. As mentioned earlier, the total GxE effect was measured as the ratio of the RRs in each stratum of vitamin use (i.e., RR Vitamins /RR No vitamins ). RRRs and 95% CIs were 0.56 (0.45-0.69) with the variant allele at rs1339221 (q = 0.011), 0.57 (0.46-0.70) with the variant allele at rs11117745 (q = 0.011), and 0.62 (0.50-0.76) with the variant allele at rs2099557 (q = 0.037) ( Table 6). Next, we considered haplotypes of the top three SNPs in ESRRG (Table 7) and analyzed the following two-SNP and three-SNP combinations: (i) rs2099557-rs1339221, (ii) rs1339221-rs11117745, and (iii) rs2099557-rs1339221-rs11117745. The RRRs and 95% CIs were 0.52 (0.41-0.66) for SNP combination (i), 0.55 (0.44-0.68) for SNP combination (ii), and 0.50 (0.40-0.64) for SNP combination (iii) ( Table 7). The pattern of RRRs was similar to that in the single-SNP analyses. Provided that vitamin use itself is not harmful, the haplotypes with the lowest p-values appeared to be protective among vitamin-takers, but detrimental among non-takers ( Table 7). The RRs and RRRs were even further away from 1, and the p-values were lower than in the single-SNP analyses. The most prominent effect was with the "tt-t" haplotype in the analysis of SNP combination (iii) above, suggesting that the observed effects with SNP combinations (i) and (ii), as well as those observed in the single-SNP analyses, are likely to be driven by this haplotype.
Our in silico analyses using Hetionet revealed that the cleft lip phenotype is influenced by vitamins A and D through a network of genes connected to ESRRG (Figure 9). Except for "Platelet derived growth factor subunit A" (PDGFA), few of the genes in the network have previously been linked to orofacial clefts. Figure 9 shows that ESRRG is expressed in three major tissues/organs and that the cleft lip phenotype is localized to the embryo, head, and telencephalon (the most highly developed and anterior part of the forebrain). ESRRG regulates or interacts with various genes that are regulated by the levels of vitamin A and D. Note that several genes are up-or down-regulated in the telencephalon (marked by blue arrows in Figure 9) and they also interact directly with ESRRG.

DISCUSSION
The current genome-wide search for GxE effects in isolated CL/P is based on the largest available GWAS dataset on orofacial cleft triads (Beaty et al., 2010). We identified a statistically significant GxE effect between maternal periconceptional vitamin use and genetic variants in ESRRG in the pooled analysis. We also identified potential GxVitamin and GxAlcohol effects in the European sample (Tables 3, 5), but none of the SNPs were outside the 95% confidence bands in the QQ plots (Figures 2, 4), nor were they located in or near genes with obvious connections to orofacial clefts. Because of a lack of observations, we were unable to perform GxSmoke and GxAlcohol analyses in the Asian sample. The low level of alcohol intake and cigarette smoking appears to be a general trend among Asian women, and is likely to be even lower among those who are pregnant or planning pregnancy (Yang et al., 2010;Ng et al., 2014;Haaland et al., 2017). Although GxVitamin analysis of the Asian sample was not hampered by a lack of observations, vitamin intake was nevertheless markedly lower in this group compared to European women (14 vs. 56%) (Figure 1).
The current genome-wide scan for GxE effects in isolated CL/P was motivated by a number of observations. First, several studies have reported associations between orofacial clefts and periconceptional maternal cigarette smoking, alcohol intake and vitamin use (reviewed in Jugessur et al., 2009a;Dixon et al., 2011;Marazita, 2012;Rahimov et al., 2012;Leslie and Marazita, 2013;Beaty et al., 2016), but our comprehensive analysis of 334 autosomal cleft candidate genes showed little evidence of interaction with these maternal exposures despite being the largest GxE study of clefts at the time (Skare et al., 2012). Second, the lack of GxE effects may be due to a combination of limited statistical power to detect nothing but the largest GxE effects and the a priori small environmental contributions to CL/P (∼9%) and CPO (∼10%) (Grosen et al., 2011). Because most previous GxE studies have used a candidate-gene approach and are based on relatively small sample sizes, the small environmental contributions are likely to have further reduced the power to detect a GxE effect. Third, only two genome-wide studies of GxE effects have so far been published in orofacial clefts-in isolated CPO Wu et al., 2014). We thus focused on the larger sample of isolated CL/P and screened for GxE effects using the same GWAS dataset as in Wu et al. (2014). The connection between vitamins and variants in ESRRG is novel in the context of orofacial clefts. To shed more light on how vitamins, ESRRG and clefting might relate to one another, we performed in silico analyses using the Hetionet database (Figure 9). By harnessing data from several publicly available meta-databases, Hetionet generates a detailed overview of the relationships between a given disorder/disease, genes and compounds that might easily be overlooked if the focus were solely on specific aspects of disease-gene associations (Greene et al., 2015;Himmelstein et al., 2017). Our analyses revealed a rich network of genes connecting cleft lip to ESRRG and to two vitamins in particular-vitamins A and D. These genes are significantly influenced by the levels of these two vitamins, which might partially explain why the extra vitamin intake by pregnant mothers appears to protect the fetus against clefts. There is some evidence in the literature linking vitamin A   itself and genes related to vitamin A, e.g., retinoic acid receptor alpha, RARA, with the risk of orofacial clefts (Rothman et al., 1995;Mitchell et al., 2003;Bille et al., 2007;Johansen et al., 2008;Boyles et al., 2009;Skare et al., 2012). For instance, early studies in mice showed that the timing of exposure to retinoic acid (a metabolite of vitamin A) is important   FIGURE 8 | Regional association plot to assess the strength of association and regional information flanking the lead SNP in ESRRG, rs1339221, shown here in blue alongside its p-value.
in the disruption of the expression patterns of key growth factors, resulting in abnormally small palatal shelves that cannot fuse (Abbott et al., 1989(Abbott et al., , 2005Abbott and Birnbaum, 1990). Compared to vitamin A or vitamin B complex, few studies have examined the effects of high or low dose of vitamin D on clefting risk. Unfortunately, the cleft consortium (Beaty et al., 2010) did not include detailed information on the use of different types of vitamins, which would have allowed a more targeted analysis of vitamins A and D in relation to the ESRRG variants.
The impact of specific variants in ESRRG on the risk of orofacial clefts is also uncharted territory, but several lines of evidence point to a biologically plausible link between estrogens, ESRRG and craniofacial malformations. First, estrogens are a group of steroid-based sex hormones that are involved in several important developmental and physiological processes, including cartilage proliferation and growth, and formation of the craniofacial complex (Ahi, 2016). Second, sex hormones are involved in several traits associated with sexual dimorphism (Callewaert et al., 2010;Randall et al., 2013;Sanger et al., 2014). Given the consistently observed higher male-to-female ratio of isolated CL/P (∼2:1 in Caucasians), it is plausible that the skewed sex prevalence is a manifestation of opposing sex steroid actions. Third, exposure to very high or very low doses of estrogens during embryonic development results in craniofacial skeletal malformations in various animal models (Fushimi et al., 2009;Cohen et al., 2014;Morthorst et al., 2016). For instance, when zebrafish are exposed to bisphenol-A (an endocrine-disrupting chemical that mimics estrogen), they develop craniofacial malformations (Kramer et al., 1990). Signaling and dosage regulation of estrogens are finely orchestrated by estrogen receptors (ERs) and estrogen related receptors (ESRRs). These two closely-related families of receptors share target genes, coregulators and promoters (Maglott et al., 2005). ESRRG, the gene identified in our GxVitamin analyses, encodes the orphan nuclear receptor "Estrogen-related receptor γ" (ERRγ). ERRγ itself does not appear to be important for skeletal development, but it is a sex-dependent negative regulator of postnatal bone formation (Cardelli and Aubin, 2014). The majority of the genes in Figure 9 have not previously been linked to orofacial clefts, except perhaps for PDGFA. This gene belongs to the PDGF family of genes that play important roles in the PDGF receptor-alpha (PDGFR-α) signaling pathway. Compared to PDGFA, PDGFC has a well-substantiated role in palatogenesis (Eberhart et al., 2008). Mice with the Pdgfc gene knocked out (Pdgfc −/− ) exhibit a complete cleft of the secondary palate (Ding et al., 2004). The phenotype is less severe in Pdgfa −/− mice. The ligands PDGFA and PDGFC share a common pathway with the PDGFR-α receptor in regulating the development of craniofacial structures. Targeted deletion of murine Pdgfra causes neural tube closure defects, including midfacial and palatal clefts (Soriano, 1997). Retinoic acid, on the other hand, inhibits proliferation of embryonic palatal shelves in mice by downregulating Pdgfc activity (Han et al., 2006). Pedigo et al. (2007) later reported the location of a complex retinoic acid response element in a region upstream of the transcription start site of PDGFA that was previously shown to harbor basal and vitamin D-inducible enhancer activity, which Frontiers in Genetics | www.frontiersin.org FIGURE 9 | Relationships between cleft lip, vitamins and ESRRG, visualized as a heterogenic network, where each relationship evidenced in public databases is depicted as an arrow connecting two nodes. Different types of nodes (e.g., genes, diseases, etc.) are connected by different types of relationships (e.g., regulation, interaction, etc.). The nodes and arrows are colored according to the legend at the top of the figure. The thickness of the "DOWNREGULATES_AdG" arrow is less than that of "UPREGULATES_AuG," and similarly for "DOWNREGULATES_CdG" and "UPREGULATES_CuG." The numbers on the two latter arrows show the Z-score, which is a measure for the strength of the influence. The figure was generated using the Hetionet network (see the Materials and Methods section for details).
lends support to the connections seen between these entities in Figure 9.
In addition to the novel insights from the biological findings, we aimed at developing a robust method for genome-wide screening of GxE effects in GWAS data. Using a log-linear model with multiplicative dose-response is a very efficient statistical approach, and case-parent triads provide sufficient information to reconstruct haplotypes with high precision. Yet, our study may still lack power to detect effects that are small. Figure 10 depicts power simulations for different sample sizes and different proportions of exposed to unexposed mothers, reflecting those that were available in the current GWAS dataset. The figure shows that there is acceptable power at a nominal significance level of 5% to detect a single-SNP GxE effect with relative risk ratio (RRR) of 1.4 or higher in a sample consisting of 2,000 triads, and an effect with RRR of 1.6 or higher with 1,000 triads. The proportion of exposed to unexposed cases does not have a major impact on power as long as it does not deviate substantially from 1-to-1. In a GWAS setting, there is of course the extra burden of extensive multiple testing. Controlling the FDR in our study for each exposure is less taxing on statistical power and better suited for discovery than, for instance, a strict Bonferroni correction. However, it increases the need for independent verification of top hits.
Our study has limitations. Even though it was based on the largest collection of isolated CL/P triads to date, it still had limited statistical power in the smallest group of exposed mothers. This was especially true for the GxSmoke and GxAlcohol analyses in the Asian sample, where the low level of exposure prevented any meaningful GxE analysis. The lack of a replication cohort was also a major shortcoming. Our study is the first to investigate the risk of isolated CL/P using GWAS data in two major ethnicities, and there are currently no other similar studies on isolated CL/P that can be used to confirm our findings. Further, our assumption of a detectable RRR of 1.4 with a sample size of 2,000 triads and an RRR of 1.6 with 1,000 triads may be optimistic in the context of birth defects and complex traits in general. The current scan for GxE effects in isolated CL/P should therefore be considered exploratory and hypothesisgenerating at this stage, assuming that other researchers who in the future have access to comparable cleft data and more detailed information on different types of vitamins (including vitamins A and D) would be interested in replicating these findings.
To summarize, we identified significant interactions between variants in ESRRG and vitamin use in the pooled analysis. Our in silico analyses revealed an intricate network of genes linking cleft lip, ESRRG and two vitamins in particular: vitamins FIGURE 10 | Single-SNP power for different sample sizes and proportions of exposed to unexposed mothers. The figure shows the power for detecting a GxE effect with increasing relative risk ratios (RRR). The relative risk in the smallest exposure group was increased, but the relative risk in the largest exposure group was set equal to 1. Additionally, we used a minor allele frequency (MAF) of 0.2, a nominal significance level of 5%, and a total sample size of 1,000 triads and 2,000 triads with different proportions of exposed/unexposed mothers. Note that the minor allele was used as the risk allele.
A and D. These GxE effects are novel and warrant further investigations to unravel the potential interplay between vitamins and ESRRG variants. If confirmed in other cleft samples, they could provide prospects for exploring the impact of estrogens and vitamins on orofacial clefts, with potential translational applications.

AUTHOR CONTRIBUTIONS
AJ, HG, RL, ØH, MG, and JR: Conception of the work, study design, data analysis and interpretation, manuscript preparation and final manuscript approval; HG, MG, JR, and ØH: Statistical modeling and software design; AJ, HG, and RL: Funding acquisition.