A Three-Way Interaction among Maternal and Fetal Variants Contributing to Congenital Heart Defects.

Congenital heart defects (CHDs) are the most prevalent and severe type of birth defects, occurring in approximately 1 of every 100 live births ( Moller et al. , 1993 , Botto et al. , 2001 , Hoffman & Kaplan, 2002 , Hoffman et al. , 2004 ). A substantial portion of infants born with serious CHDs die in infancy, and many of those who survive may require repeated surgeries and lengthy hospitalizations ( Cleves et al. , 2003 , Gilboa et al. , 2010 , Nembhard et al. , 2001 ). CHDs are a major concern for pediatric morbidity and mortality, and there are currently few strategies for reducing the public health impact of these conditions ( Jenkins et al. , 2007 ). The genetic susceptibility of CHDs has been recognized for decades ( Hobbs et al. , 2002 ). As demonstrated by twin studies, the concordance rates of CHD phenotypes were significantly higher among monozygotic twins (10.0%) than dizygotic twins (2.5%) ( Berg et al. , 1989 ). Animal studies have identified several biological pathways controlling the development of the fetal heart. For example, two transcription factors cytoplasma 1 (the nuclear factor of activated T-cells; NF-ATc) and Smad6 were suggested to be involved in the formation of cardiac valves ( Galvin et al. , 2000 ). Mice lacking NF-ATc exhibit fatal defects in valve formation and disruption of Smad6 leads to abnormally thickened valves ( High & Epstein, 2008 ). Meanwhile, maternal environmental exposures may have an effect modification on key developmental genes that shape cardiac development ( Ashworth & Antipatis, 2001 , Keen et al. , 2003 , Doolin et al. , 2002 ). It is believed that more than 85% of CHDs result from complex interactions between genetic variants, epigenetic modifications and maternal environmental exposures ( Botto & Correa, 2003 ). To identify genetic and environmental factors that foster the development of CHDs, extensive candidate-gene-based and genome-wide association studies have been conducted ( Hobbs et al. , 2006 , Goldmuntz et al. , 2008 , Hobbs et al. , 2011 , Wessels & Willems, 2010 , Pediatric Cardiac Genomics et al. , 2013 , Zaidi et al. , 2013 ). Multiple genes and genetic variants are associated with CHDs. However, most of these studies, including our own, have adopted a single-locus analysis strategy, by testing each single genetic variant individually. As discussed above and elsewhere, CHDs are usually caused by the interplay between multiple genetic variants and environmental factors. If a gene operates primarily through a complex mechanism involving multiple other genes, the effect may often be missed if one examines it in isolation, without allowing for its potential interactions with other genes ( Cordell, 2009 ). Particularly for maternal and perinatal research, two types of gene-gene interactions are possible during pregnancy: intra-generational interaction within either maternal or fetal genes, and inter-generational interaction between maternal and fetal genes ( Sinsheimer et al. , 2010 ). Both maternal and fetal genes may interact with other genes within each individual genome, while the maternal and fetal genes may also interact across genomes, leading to either conflicting or beneficial environment for fetal growth, which may influence the phenotypes of both mothers and babies ( Haig, 2004 ). To our knowledge, only one intra-generational interaction between variants from maternal genes MTHFR and CBS has been reported for association with CHDs ( Lupo et al. , 2010 ). We and others have proposed to detect inter-generational interactions with a penalized logistic regression model, and applied it to the studies of CHDs ( Li et al. , 2010 , Li et al. , 2014b , Li et al. , 2014a ). However, those studies are limited to interactions between genes from the same genomic region. Further, all existing interaction studies have adopted a two-way interaction strategy, by testing the interaction effect between two genetic variants or interaction effect between one genetic variant and one environmental factor ( Lupo et al. , 2010 , Li et al. , 2014a , Li et al. , 2014b , Tang et al. , 2014 , Tang et al. , 2015 ). Although convenient and easy to interpret, such a strategy cannot consider higher order interactions (e.g., three-way interactions), or the complex gene-by-gene (GXG) and gene-by-environment (GXE) interactions between maternal variants, fetal variants and environmental exposures. In reality, higher-order interactions (e.g., three-way interactions) could exist ( Moore, 2003 , Ritchie et al. , 2001 ). By not limiting the search to two-way interactions, we may have a better chance of identifying potentially important GXG/GXE interactions, which may elucidate how maternal and fetal genetic variants interact with one another and interplay with environmental factors to cause CHDs. However, detecting high-order GXG/GXE interactions underlying complex human diseases have remained to be a major challenge in genetic studies ( Cordell, 2009 , Thomas, 2010 ). A main difficulty for the conventional regression-based methods, such as logistic regression model, is the specification of underlying disease models. By specifying a disease model, certain inheritance of disease risks is implicitly assumed. When such assumption is violated, the effect estimate could be biased and the type I errors be inflated. Further, the number of parameters to be specified in the model increases exponentially with the order of interaction being analyzed. It may become impractical to specify such a model for high dimensional data ( Carlborg & Haley, 2004 , Peduzzi et al. , 1996 ). To address these difficulties, non-parametric dimension-reduction methods have gained popularity for GXG/GXE research in the past few years, especially for detecting high-order GXG/GXE interactions ( Ritchie et al. , 2001 , Lou et al. , 2007 , Lu et al. , 2012 , Li et al. , 2011 ). These proposed methods have been widely applied to population-based studies of complex human diseases, such as breast cancer and drug addiction. It would be of great interest to extend these methods for detecting GXG/GXE interactions in maternal and prenatal research. Our current study is motivated by this gap in the existing literature, aiming to detect joint actions among a large number of maternal variants, fetal variants and maternal environmental exposures, while allowing for high-order statistical interactions. The current study included 872 maternal SNPs, 872 fetal SNPs selected from 62 candidate genes, and 19 maternal environmental exposures. We applied a recently developed likelihood-ratio Mann-Whitney (LRMW) method to search for GXG/GXE combinations among a large number of genetic and environmental factors. Finally, we explore possible biological mechanisms with respect to the GXG/GXE combinations that jointly alter CHD risk.

We applied the LRMW method to our dataset, aiming to search for GXG/GXE combinations among 872 maternal SNPs, 872 fetal SNPs, and 4 maternal environmental exposures, including folic acid supplementation, alcohol drinking, smoking and obesity. These maternal lifestyle factors have been shown to be associated with CHD development ( Hobbs et al. , 2011 ). For example, maternal intake of folic acid containing supplements may reduce the risk of CHDs ( Shaw et al. , 1995 , van Beynum et al. , 2010 ), while smoking, drinking and obesity may increase the risk ( Jenkins et al. , 2007 , Malik et al. , 2008 , Stothard et al. , 2009 ). After tenfold cross validation, LRMW identified a final model with 3 SNPs, which gave the maximum LRMW statistic, indicating a joint association among 3 SNPs. None of the environmental factors was selected into the final model. The 3 SNPs are fetal rs625879, maternal rs2169650 and maternal rs8177441, located within gene BHMT 2, GSTP 1 and GPX 3, respectively. Eight GXG groups were formed by the combination of three SNPs. The three SNPs and corresponding grouping strategies are summarized in Table 2 . Three SNPs were identified through forward selection, and the corresponding GXG groups had an increasing area under the curve (AUC), consistent with an increasing ability to predict cardiac defects (i.e. from 0.54, 0.58 to 0.60). It is worthwhile to note that LRMW method is entirely non-parametric. It allows complex interactions without specifying any interaction model explicitly. Therefore, it is possible for the identified SNPs to have a joint action through either additive or interaction effects. In order to examine the underlying genetic mechanism, we looked into each SNP pair to identify patterns of risk for CHDs among the corresponding GXG groups. For each SNP pair, all mother-offspring pairs can be partitioned into 4 groups based on their GXG combinations. For example, fetal SNP rs625879 and maternal SNP rs2169650 together form four two-SNP genotypes: G 1 ={fetal rs625879 = CC and maternal rs2169650 = GG }, G 2 ={fetal rs625879 = CC and maternal rs2169650 = AA or AG }, G 3 ={fetal rs625879 = AA or AC and maternal rs2169650 = GG }, and G 4 ={fetal rs625879 = AA or AC and maternal rs2169650 = AA or AG }. We found that the effect of fetal SNP rs625879 was modified differently by the genotypes of maternal SNP rs2169650 (see. Figure 1 ), indicating a potential interaction effect between the two SNPs. In particular, it should be noted that this interaction is “essential” and not completely removable by a monotonic transformation of the data ( Wu et al. , 2009 ). The examination of the other two SNP pairs also indicated an additive effect between fetal SNP rs625879 and maternal SNP rs8177441 (see Figure 2 ), and an interaction effect between maternal SNP rs2169650 and maternal SNP rs8177441 (see Figure 3 ). The proposed LRMW method is non-parametric, and does not assume any mode of inherence of disease. However, for the same reason, it does not provide any estimation of effect sizes. In order to estimate the effect sizes among identified genotype groups, we fit a logistic regression model based on the observed patterns of disease risk. The final model included main effects of three SNPs, and two interaction effects: 1) between fetal SNP rs625879 and maternal SNP rs2169650; 2) between maternal SNP rs2169650 and maternal SNP rs8177441. The results are summarized in Table 3 . The results showed that a mother-offspring dyad had the highest risk of disease (OR=2.47; 95% C.I. [1.20, 5.08]), when genotype CC , AG or AA , and GG were observed for fetal rs625879, maternal rs2169650 and maternal rs8177441, respectively. The logistic regression attained a significant joint association for three SNPs ( P -value=1.13e-07). It should also be noted that this is a nominal p -value without accounting for the selection of SNPs in LRMW method. Permutation test of LRMW with 5,000 replicates obtained an empirical p- value of 0.11, which did not reach the conventional 0.05 threshold. However, we think it may still suggest a potential significant association, considering the fact we had limited sample size (i.e. 1,498 samples) but a very high-dimensional search space (i.e. 1,748 genetic and environmental factors together and their high-order combinations).

Our study builds on previous hypotheses and studies demonstrating that the etiology of nonsyndromic CHDs is a process involving multiple genetic variants, both maternal and fetal, maternal environmental exposures as well as complex GXG and GXE interactions. Existing studies typically have examined each single genetic variant for its association with CHDs. Relatively few studies have been conducted to examine potential CHD-related GXG and GXE interactions. Our understanding of the complex GXG/GXE interaction underlying CHDs is still in its infancy, especially for high-order interactions. Further, detecting those GXG/GXE interaction is especially challenging in maternal and prenatal research, since both intra-generational interaction and inter-generational interactions may occur, and both maternal and fetal genes may further interact with maternal environmental exposures ( Sinsheimer et al. , 2010 ). Our current study is motivated by the limitations and challenge in existing studies, aiming to detect complex GXG/GXE interactions among maternal variants, fetal variants and maternal environmental exposures, with the consideration of possible high-order interactions. Although, none of the environmental factors was identified to be significant in this study, to our knowledge, our finding is the first three-way interaction among maternal and fetal variants, involving two types of interactions: an intra-generational interaction within maternal genome (i.e. maternal GSTP 1 and GPX 3) and an inter-generational interaction between maternal and fetal genome (i.e. maternal GPX 3 and fetal BHMT 2). Findings from previous studies indicated that the interactions between maternal genes in the folate metabolic pathway may influence the risk of CHDs ( Lupo et al. , 2010 ), and interactions between maternal and fetal genes may also influence birth defects, such as neural tube defects ( Lupo et al. , 2014 ). The findings in our study are in line with previous findings, further suggesting that the interaction effects among maternal and fetal genes could be complex, and genes from homocysteine and transsulfuration may interact with one another to influence CHD development jointly. The current study also extends our previous investigation of maternal-fetal genotype interactions to variants from different genomic regions ( Li et al. , 2014b , Li et al. , 2014a ). For comparison purpose, we also applied a conventional logistic regression model to the same dataset. Similar to our LRMW method, a forward selection strategy was used for logistic regression by minimizing Akaike Information Criteria (AIC). By using a logistic regression model, a total number of 358 SNPs were selected, forming a highly complex final model. It will require a large number of parameters to further investigate their interactions. In addition, it is less straightforward for interpretation. In our study, we identified a joint action among three variants, fetal SNP rs625879, maternal SNPs rs216950 and rs8177441. It has been suggested that SNP rs625879 was associated with various complex diseases and birth defects, such as coeliac disease ( Hozyasz et al. , 2012 ), endometriosis-associated infertility ( Szczepanska et al. , 2011 ) and orofacial clefts ( Mostowska et al. , 2010 ). Further, the three identified SNPs are located in gene BHMT 2, GSTP 1 and GPX 3, respectively. Gene BHMT 2, or betaine--homocysteine S-methyltransferase 2, encodes one of two methyl transferases that can catalyze the transfer of the methyl group from betaine to homocysteine, which plays a crucial role in methylation reactions. Anomalies in homocysteine metabolism have been implicated in disorders ranging from vascular disease to birth defects, such as spina bifida and CHDs ( Giusti et al. , 2010 , McGeachie et al. , 2009 , Weisberg et al. , 2003 , Shaw et al. , 2009 ). Gene GSTP 1, or glutathione S-transferase pi 1 , belongs to a family of enzymes that play an important role in detoxification by catalyzing the conjugation of many hydrophobic and electrophilic compounds with reduced glutathione. GSTP 1 proteins were suggested to function in xenobiotic metabolism and play a role in susceptibility to various diseases, such as esophageal squamous cell carcinoma and gastric cardia cancer ( Zendehdel et al. , 2009 ). Gene GPX 3, or glutathione peroxidase 3 (plasma) , belongs to the glutathione peroxidase family, which functions in the detoxification of hydrogen peroxide. GPX 3 protects cells and enzymes from oxidative damage, by catalyzing the reduction of hydrogen peroxide, lipid peroxides and organic hydroperoxide. It was also suggested that decreased GPX 3 activity may lead to inadequate nitric oxide (NO) levels, which disrupts platelet inhibitory mechanisms and increases arterial thrombosis ( Kenet et al. , 1999 ). GPX 3 was also found to be associated with ischaemic stroke among children and young adults ( Voetsch et al. , 2007 ). While it is biologically plausible that those three identified genes may jointly alter the risk of CHDs, it is necessary for further studies to validate and replicate this finding. A few limitations should also be noted. First, the permutation test of LRMW method attained an empirical P -value of 0.11, which did not reach the commonly used 0.05 threshold. We think this can be largely due to our limited sample size (i.e.1,498 samples) and high-dimensional search space (i.e. 1,748 genetic and environmental factors together with consideration of their high order combinations). However, NBDPS is the largest study ever conducted in the U.S., and we have included all the available samples collected to date. We expect power to increase as NBDPS moves further with more samples. Second, no environmental factors were identified in our study. One possible reason is that our LRMW method utilized a forward selection strategy, which search for GXG/GXE combinations forwardly. A major advantage is that the forward search is computational efficient and feasible for high dimensional data, such as genome-wide data. However, to detect an interaction effect between two factors, our method requires at least one of them to have a main effect. In our framework, a genetic variant or environmental factor without a main effect may only be introduced into the final model if it interacts with other factors having significant main effects ( Lu et al. , 2012 ). An interaction effect without any main effect of all factors will require exhaustive search among all factors, which may substantially increase the computational intensity and possibly reduce power ( Li et al. , 2011 ). We also expect “pure interaction” to be a rare scenario in disease etiology. Third, we have combined conotruncal heart defects (332 case-mother pairs) and obstructive heart defects (291 case-mother pairs) as cases in our study. It is possible that disease heterogeneity may exist between two subtypes of CHDs. However, previous studies have demonstrated that various types of pediatric heart disease may share same causal genes. For example, cardiac transcription factor, NKX2.5, has been associated with a wide variety of cardiac phenotypes such as, atrial septal defect, atrioventricular block, double outlet right ventricle, tetralogy of fallot, and ventricular septal defect ( Benson, 2010 ). Further, stratified analysis with each particular subtype will significantly reduce our sample size, limiting our power to detect interactions, especially for high-order interactions.

The study was approved by University of Arkansas for Medical Sciences’ Institutional Review Board and the National Birth Defects Prevention Study (NBDPS), with protocol oversight by the Centers for Disease Control and Prevention (CDC) and National Center on Birth Defects and Developmental Disabilities. All study subjects gave informed consent. For minors, informed consent was obtained from their legal guardian or themselves, if emancipated. All subjects were participants of the National Birth Defects Prevention Study, the largest population-based case control study conducted in the US of nonsyndromic birth defects, covering an annual birth population of 482,000, or 10 % of U.S. births. CHD cases were ascertained from birth defect registries in ten participating states that had similar inclusion criteria: Arkansas, California, Georgia, Iowa, Massachusetts, New Jersey, New York, North Carolina, Texas, and Utah. NBDPS began in 1996, and a detailed description of the study can be found elsewhere ( Yoon et al. , 2001 , Rasmussen et al. , 2002 , Gallagher et al. , 2011 ). In the current study, we included all available mother-offspring pairs with estimated dates of delivery between 1997 and 2010, which is based on the due date reported by the mother during interview. Case pairs were defined as those in which the child had at least one type of conotruncal heart defect or obstructive heart defect. Control pairs were defined as those in which the child had no structural birth defect. Control families were randomly selected from birth certificate and/or birth hospital records and thus represent a random sample from the general population ( Yoon et al. , 2001 ). For the current analyses, the study population comprised 1,498 mother-offspring pairs, including 623 case pairs and 875 control pairs. The maternal characteristics are summarized in Table 1 . The cases and control families were balanced with respect to maternal age, race and ethnicity, maternal education, household income, folic acid intake, alcohol consumption, cigarette smoking and maternal obesity (all p -values>0.05). Our research team commissioned a custom 1,536 SNP panel covering 62 genes in the homocysteine, folate, and transsulfuration pathways potentially related to the development of CHDs, using the Illumina® GoldenGate custom genotyping platform, as detailed elsewhere ( Chowdhury et al. , 2012 ). Laboratories that use NBDPS archival DNA samples must demonstrate their proficiency in genotyping techniques by passing an External Quality Assessment (EQA). The EQA for High-throughput genotyping platforms, such as the Illumina Golden Gate assay employed in this project, consists of genotyping of a single SNP in one externally supplied blood-buccal trio. Results from the laboratory must be concordant (99%) for paired blood-derived and buccal-derived DNA and between genomic and whole genome amplified DNA. In addition, inter-lab results must be concordant as well as concordance between pre-characterized DNA and third-party results. Negative controls must not yield results and genotypes of trios must be consistent with Mendelian inheritance. Arkansas Center for Birth Defects Research and Prevention has consistently scored 100% on NBDPS EQA. GenomeStudio, the software toolkit developed by Illumina, was used for the initial genotype calling. We found that the quality of genotype clustering varied substantially among SNPs, which we attribute to the in silico design of the custom SNP panel without the subsequent quality checks that would be applied to a standard commercial array. To ensure the data quality, we updated the genotype clustering by using SNPMClust, a bivariate Gaussian model-based genotype clustering and calling algorithm, which is available as an R package on the Comprehensive R Archive Network (CRAN). After running SNPMClust, clustering and classification plots for all SNPs were visually inspected, leading to dropping a SNP from analysis or running SNPMClust under non-default settings in some cases. To ensure high-quality genotypes, we applied stringent quality control measures and excluded SNPs with obviously poor clustering behavior, no-call rates > 10%, Mendelian error rates > 5%, minor allele frequencies < 5%, or significant deviation from Hardy-Weinberg Equilibrium in at least one racial group ( P -value < 10 −4 ). After the quality assessment, the analytical data had 872 variants in total. Recently, we and others have proposed a likelihood ratio Mann-Whitney (LRMW) method for detecting genetic variants and environmental factors that are jointly associated with disease phenotypes ( Lu et al. , 2012 ). The LRMW method searches for joint action among a large number of genetic variants, such as genome-wide data, and environmental factors. It extends the traditional univariate Mann-Whitney test to assess the joint association of multiple genetic variants and environmental factors simultaneously, while allowing for high-order statistical interactions. We have demonstrated through simulations that the LRMW method has improved power over the conventional Mann-Whitney method when higher order interaction exists among loci. The application of the LRMW method also identified a four-locus interaction among SNPs for type II diabetes, which was replicated in an independent dataset ( Lu et al. , 2012 ). The method has been implemented in a C++ software package, referred to as genome-wide gene-gene interaction analyses (GWGGI) ( Wei & Lu, 2014 ). In this article, we extend the LRWM method in the context of mother-offspring pair data, aiming to detect the joint action among maternal variants, fetal variants and maternal environmental exposures. Following the same notation, we briefly describe our method below. The theoretical details can be found elsewhere ( Lu et al. , 2012 ). Assuming that multiple SNPs and environmental factors may be jointly associated with the disease outcome, we expect all mother-offspring pairs can be partitioned into various groups based on their gene-by-gene and gene-by-environment (GXG/GXE) combinations, among which the likelihood of disease may not be all equal. Determining GXG/GXE groups on high-dimensional data could be statistically complicated and computationally intense, since a large number of GXG/GXE groups can be formed. To search for GXG/GXE groups that are most significantly associated with disease phenotype, we adopt a forward selection algorithm, which has shown promising performance in previous studies ( Li et al. , 2011 , Lu et al. , 2012 , Lu & Elston, 2008 ). The algorithm starts with a null model by including all mother-offspring pairs as one group, and then uses forward selection to introduce SNPs or environmental factors into the model to partition all mother-offspring pairs into various GXG/GXE groups. In the first step, all SNPs and environmental factors were examined one at a time, and the one factor with the maximum LRMW statistic described below was selected. For example, if a bi-allelic SNP with alleles denoted as A and a , is selected in the first step, all pairs can be partitioned into two groups in two possible ways, assuming either dominant or recessive inheritance, respectively: G 1 ={ AA or Aa } and G 2 ={ aa } (carrying one copy of risk allele will increase disease risk) G 1 ={ AA } and G 2 ={ Aa or aa } (carrying two copies of risk alleles will increase disease risk) After that, an environmental factor with two levels (e.g. smoker vs non-smoker) can further partition all pairs into 4 groups (i.e. G 1 /smoker, G 1 /non-smoker, G 2 /smoker, and G 2 /non-smoker). Or, a second SNP with alleles denoted as B and b can further partition all individuals into 4 groups in two possible ways: G 1 / BB or Bb , G 1 / bb , G 2 / BB or Bb , G 2 / bb (Assuming a dominant inheritance of allele B ) G 1 / BB , G 1 / Bb or bb , G 2 / BB , G 2 / Bb or bb (Assuming a recessive inheritance of allele B ) Furthermore, because the first SNP and the second SNP may be a maternal SNP and a fetal SNP (or vice versa), respectively, the selected GXG/GXE groups may represent joint gene actions within maternal or fetal genomes or the joint gene actions between maternal and fetal genomes. In practice, the SNPs and environmental factors are selected forwardly, forming GXG/GXE groups adaptively to maximize a LRMW statistic defined below. Suppose the selected SNPs and environmental factors can form R GXG/GXE groups, G 1 , G 2 ,…… G R , we first define a likelihood ratio ( LR ) for each risk group: L R ( G r ) = P ( G r ∣ D ) / P ( G r ∣ D ¯ ) ; 1 ≤ r ≤ R ; where D represents cases and D̄ represents controls. The LRMW statistic can be defined as: U LRMW = ∑ i = 1 R ∑ j = 1 R n G i D n G j D ¯ φ [ L R ( G i ) , L R ( G j ) ] ; where n G i D and n G j D ¯ are the number of case-mother pairs carrying genotype G i and the number of control-mother pairs carrying genotype G j , respectively; and φ[.] is a kernel function to compare the LR risk score between two groups. A commonly used kernel function is φ [ L R ( G i ) , L R ( G j ) ] = { 1 i f L R ( G j ) L R ( G i ) ; which leads to a Mann-Whitney statistic comparing the difference in LR risk scores between cases and controls ( Mann & Whitney, 1947 , Wilcoxon, 1945 ). In our study, we model each mother-offspring pair as a unit, and build a LRWM statistic for case-mother and control-mother pairs. In order to avoid over-fitting of data, we conduct a 10-fold cross-validation, i.e., randomly dividing the study population into 10 subsets, to determine the most parsimonious statistical model based on GXG/GXE groups. The joint association is also evaluated by permutation testing, which randomly shuffles the phenotype and applies the same procedure, including forward selection of GXG/GXE groups and ten-fold cross-validation. An empirical p -value can be obtained by repeating the permutation a large number of times (e.g. 5,000 times).