A genetic association study detects haplotypes associated with obstructive heart defects.

Congenital heart defects (CHDs) are the most prevalent type of birth defect and the leading cause of infant mortality attributable to birth defects ( Yang et al. 2006 ). Worldwide, 1.35 million infants are born with CHDs per year, with an estimated prevalence of 8.1 per 1,000 live births in North America, and 9.3 per 1,000 live births in Asia ( Reller et al. 2008 ; van der Linde et al. 2011 ). CHDs have many subtypes, and obstructive heart defects (OHDs) are a large subgroup of CHDs with obstructive lesions in blood vessels or valves in the right or left heart, such as pulmonic stenosis, aortic stenosis, coarctation of the aorta, bicuspid aortic valve stenosis, and subaortic stenosis. OHDs account for approximately 25 % of all CHDs, and our previous study has shown that OHDs are the most lethal subtype of heart defect in infancy ( Cleves et al. 2003 ). Almost 40 % of infants born with serious OHDs 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 ). Despite recent medical and surgical advances, substantial OHD-attributable morbidity and mortality continues to represent a global burden of disease ( Botto and Correa 2003 ). CHDs are believed to be caused by a complex interplay between environmental exposures, genetic, and epigenetic factors. Familial aggregation provides direct evidence for the heritability of CHDs ( Hobbs et al. 2002 ), while twin studies have demonstrated that the concordance rates of CHD phenotypes are significantly higher among monozygotic twins (10.0 %) than dizygotic twins (2.5 %) ( Berg et al. 1989 ). The estimated heritability may vary by the phenotypes of cardiovascular malformation; for example, 89 % for bicuspid aortic valve and 71–90 % for left ventricular outflow tract malformations ( Cripe et al. 2004 ; McBride et al. 2005 ). In the past few years, we and others have identified multiple genetic variants associated with CHDs ( Goldmuntz et al. 2008 ; Hobbs et al. 2006 , 2011 ; Wessels and Willems 2010 ). Though these findings have provided valuable insights into the genetic etiology of CHDs, the disease-susceptibility genes identified so far only account for a small fraction of CHD prevalence ( Fahed et al. 2013 ). Moreover, relatively few studies have investigated the association between OHDs and genetic variants. The genetic architecture of OHDs has remained elusive. It is still unclear how many genetic variants are associated with OHDs, how they are distributed in the population, and how they interact with one another to cause OHDs. Investigation of the genetic mechanisms underlying the development of CHDs can be especially challenging. CHDs develop during embryogenesis, with multiple maternal–fetal metabolic, genomic, and epigenomic interactions. As discussed by Haig, a maternal–fetal unit may have three distinct haplotypes at each locus: the maternally derived fetal haplotype that is shared by the mother and fetus, the paternally derived fetal haplotype, and the non-inherited maternal haplotype ( Haig 2004 ). These three types of haplotypes do not always have identical effects during pregnancy, and may not always be mutually beneficial to the mother and her fetus. Therefore, interaction effects between the mother and her fetus may lead to either an opposing or mutually beneficial environment for fetal growth ( Sinsheimer et al. 2010 ). Over the past decade, evidence has accumulated demonstrating that maternal–fetal genotype interactions may be a common mechanism for various complex human diseases and birth defects, such as neural tube defects ( Relton et al. 2004 ), schizophrenia ( Palmer et al. 2002 ), and autism ( Zandi et al. 2006 ). However, relatively few studies have been conducted to detect the MFG interaction with regard to the development of CHDs ( Lupo et al. 2010 ). Here, we report a haplotype-based analysis to dissect the maternal, fetal, and MFG interaction effects associated with OHDs in a study population from the National Birth Defects Prevention Study (NBDPS; Yoon et al. 2001 ). The current study included 877 SNPs selected from 62 candidate genes. Haplotype blocks were defined by the SNPs in linkage disequilibrium (LD). Our method utilized the least absolute shrinkage and selection operator (lasso; Tibshirani 1996 ), a machine learning technique which allows simultaneous effect estimation and variable selection, providing an automatic inference for the underlying genetic mechanisms. No individual test is required to differentiate maternal, fetal, and MFG interaction effects. Finally, we explore possible biological mechanisms with respect to the MFG combinations that jointly alter OHD risk.

Application of our method to 140 haplotype blocks identified 9 blocks with non-zero coefficients, indicating a potential genotype–phenotype association. Information for the identified haplotypes is summarized in Table 3 . The frequencies of risk haplotypes were estimated based on the entire study population, including both cases and controls. Based on the non-zero coefficients obtained by adaptive lasso estimator, the identified blocks fell into three possible categories: maternal main effect (i.e. a m , d m ≠ 0), fetal main effect (a f , d f ≠ 0), or MFG interaction effect (i.e. i aa , i ad , i da , i dd ≠ 0). Among the nine identified haplotype blocks, four were found to have a BFDP less than 0.8, and were located in the genes MGMT , MTHFS , CBS , and DNMT3L . We further investigated the likelihood of disease for each MFG combination, and those 9 blocks exhibit three possible genetic mechanisms. In Figs. 1 , 2 , and 3 , we plot one example for each possible mechanism. Supplementary materials (Figure S1–S9) provide LD plots and plots of the inferred genetic mechanisms for all identified blocks. The results are summarized in Table 4 . A haplotype block within the MGMT gene showed a significant association with the disease. This haplotype includes six SNPs covering a 14.5 kb region on chromosome 10. The seven MFG combinations at this region correspond to three levels of disease likelihood. Two MFG combinations had the lowest likelihood of disease, and were used as reference group. We denoted the maternal–fetal genotype combination in the reference group as R 1 = { HH / HH; HH̄ / HH }. Compared to the reference group, the other MFG combinations formed two groups with increasing likelihoods of diseases, denoted as R 2 = { HH / HH̄; HH̄ / HH̄; H̄H̄ / HH̄ } and R 3 = { HH̄ / H̄H̄; H̄H̄ / H̄H̄ }. The OR between R1 and R2 was estimated to be 1.47 (95 % CI 1.11, 1.96), and the OR between R1 and R3 was 2.16 (95 % CI 1.22, 3.83). As illustrated in Fig. 1 , the fetal haplotype H showed an additive effect that decreases the risk of disease, while the risk of disease was unchanged by maternal genotypes. Similarly, our results show that the other four blocks may have a dominance effect ( Figure S2–S5 ). However, only one block within gene MTHFS had a BFDP less than 0.8. One block was located within the GSTM4 gene, comprising five SNPs on chromosome 1. The result is summarized in Table 4 . For this block, the MFG combinations can be partitioned into three risk groups, according to the maternal genotypes. We illustrate the pattern in Fig. 2 . However, this block had a BFDP greater than 0.8. Three blocks were identified with MFG interaction effect (i.e. i aa , i ad , i da , i dd ≠ 0). These three blocks were located within genes DNMT3A, CBS, DNMT3L , respectively, on chromosome 2, 21, and 21. The results are summarized in Table 4 . The block within gene DNMT3L has the smallest size, comprising two highly linked SNPs 0.6 kb apart. Based on the estimated coefficients, the MFG combinations were partitioned into three risk groups. As illustrated in Fig. 3 , when the maternal genotype is HH , the disease risk shows a decreasing pattern with the fetal genotypes. However, when the maternal genotype is HH ( H̄H̄ ), the risk of disease shows an increasing pattern with the fetal genotype. This pattern of “cross-over” was an indication of the potential MFG interaction effect. Similarly, the interactive pattern of the blocks in gene CBS and DNMT3A are illustrated in Figure S8–S9 . Among the three blocks, only CBS and DNMT3L had BFDP values less than 0.8.

Our study builds on previous publications that have reported that the genetic susceptibility of CHDs may be associated with maternal genotypes, fetal genotypes, and MFG interactions ( Goldmuntz et al. 2008 ; Hobbs et al. 2006 , 2011 ; Lupo et al. 2010 ; Wessels and Willems 2010 ). Differentiating those genetic effects from each other, however, remains a challenge in maternal and perinatal research. This challenge is partly due to the correlation between maternal and fetal genomes. Maternal or fetal effects are likely to confound each other if they are tested separately. Therefore, a single model that simultaneously includes both maternal and fetal effects is preferred ( Shi et al. 2008 ). In addition, although MFG interactions are thought to exist pervasively in the development of birth defects, genetic studies of MFG interactions and our understanding of the complex interactive mechanisms are still in their infancy ( Sinsheimer et al. 2010 ). The current study is motivated by the importance of and challenges faced by maternal and perinatal research, aiming to (1) dissect the genetic effects (i.e. maternal, fetal and MFG interactions) associated with OHDs, and (2) investigate the underlying genetic mechanisms for each identified gene (i.e. how a MFG combination influences the risk of disease). By employing an innovative logistic regression framework with adaptive lasso, we identified four genes potentially associated with OHDs. Further analyses of these results suggest that the identified genes may influence the phenotype through various genetic mechanisms. In the current study, we have conducted a haplotype- or region-based association test, which is likely to improve the power to detect a causal association by aggregating collective effects of small to moderate size, offer an opportunity to capture complex interactions among variants, and reduce the burden of multiple testing ( Li et al. 2011 ; Wang et al. 2013 ). However, compared to the conventional single-variant analysis, it is also less convenient to assess which particular variants might be driving the associations. In supplementary materials (Table S1) , we also provide the single-variant analysis results for each SNP within the identified haplotype block, using a conventional logistic regression for the corresponding genetic effect model inferred by the adaptive lasso estimator. Further extension of the haplotype-based approach to a gene-based approach is likely to share the benefit of the single-variant approach, and avoid inference of phase-ambiguous haplotypes. In our results, two haplotypes within gene MGMT and MTHFS exhibited fetal main effects only, which are significantly associated with the occurrence of OHDs (BFDP = 0.77 and 0.68, respectively). MGMT is thought to be involved in the prevention of DNA damage and oxidative stress, and previous studies have shown that MGMT gene expression is related to antioxidant mechanisms ( Niture et al. 2007 ). Previous work by our research group has also found associations between multiple maternal SNPs of MGMT and maternal metabolites, such as levels of glutamylcysteine (GluCys) and plasma folate ( Chowdhury et al. 2012 ). In particular, the plasma folate levels significantly varied by the maternal genotypes of rs10764896 (located in MGMT ) in both cases and controls, indicating a potential involvement of rs10764896 in folate-related metabolite plasma concentrations. This SNP (rs10764896) is within the haplotype block identified in the current study. The MTHFS gene encodes a key enzyme involved in the folate pathway that metabolizes formyltetrahydrofolate. The folate pathway is interconnected with methionine metabolism by providing the methyl groups for the remethylation of homocysteine back to methionine. Our previous study also found maternal SNPs within MTHFS that are related to total plasma glutathione (GSH) levels ( Chowdhury et al. 2012 ). Though experimental validation is needed to elucidate the functional mechanism of SNPs in MGMT and MTHFS , our current finding has provided additional evidence for their potential involvement in the development of CHDs. Two haplotypes within the CBS and DNMT3L genes exhibited MFG interaction effects which are significantly associated with the occurrence of OHDs (BFDP = 0.46 and 0.49, respectively). The CBS gene is located on chromosome 21q22.3, and CBS deficiency is a common cause of an inherited metabolic disorder, classical homocystinuria. The CBS gene catalyzes the conversion of homocysteine to cystathionine, the first step in the transsulfuration pathway. Studies have found that the expression of fetal CBS is concentrated in the neural and cardiac tissues, supporting a hypothesis of its involvement in embryo cardiac development ( Quere et al. 1999 ; Robert et al. 2003 ). In particular, a recent study in Han Chinese population identified a functional variant (rs2850144) in CBS that is significantly associated with the genetic susceptibility to CHDs ( Zhao et al. 2012 ). This SNP (rs2850144, BP: 43370045) was not genotyped in our study, but it is located approximately 5.5 kb from the identified haplotype block (BP 43354960–43364494). It is very likely that rs2850144 is in strong LD with the haplotype block identified in our study, and our results also suggest that the maternal CBS gene may interact with fetal CBS to jointly alter genetic susceptibility to the disease. The DNMT3L gene encodes an enzyme that stimulates de novo methylation by DNA cytosine methyltransferase 3 alpha, which is thought to be required for the establishment of maternal genomic imprintings ( Hata et al. 2002 ). This enzyme also mediates transcriptional repression through interaction with histone deacetylase 1 ( Aapola et al. 2002 ; Deplus et al. 2002 ). Studies have reported association of DNMT3L with multiple conditions, such as ovarian endometriosis ( Borghese et al. 2012 ), embryonal carcinoma ( Minami et al. 2010 ), and cervical cancer ( Gokul et al. 2007 ). Further studies will be necessary to validate the potential involvement of DNMT3L in the development of CHDs. A few limitations of our study should be noted. First, although all of the identified haplotypes had p values less than the nominal significant threshold of 0.05, none of them would remain significant with the stringent Bonferroni correction. However, Bonferroni correction is usually overly conservative, and BFDP has been suggested as a good alternative, especially in candidate gene-based studies when there is strong priori probability that the genes are likely to be associated with the disease. Second, two types of gene–gene interactions are possible during pregnancy: intra-generational interaction within either the maternal or fetal genome, and inter-generational interaction between maternal and fetal genomes. Our current analysis only considers the inter-generational interactions between maternal and fetal genes from the same genomic region (LD block). Intra-generational interaction and MFG interaction between genes from different regions are not within the scope of the current study. Third, the current study only included common SNPs that have minor allele frequencies (MAFs) of 5 % or higher. Evidence from Phase III of the International HapMap Project and 1,000 Genomes Project have supported that rare variants with lower MAFs may contribute considerately to the development of complex human diseases ( Abecasis et al. 2012 ; Altshuler et al. 2010 ). However, our current analysis tests each haplotype block, which is less easy to take into account the rare variants. Fourth, the focus of our study is to investigate the genetic susceptibility of OHDs. We are able to control the potential confounding effect of environmental factors such as BMI, but considering potential gene-by-environment interactions would substantially increase the model complexity and are not considered in the current study.

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) Center for Birth Defects and Developmental Disabilities. All study subjects gave informed written consent. For minors, informed written consent was obtained from their legal guardian. All subjects were participants of the National Birth Defects Prevention Study, an ongoing large-scale case control study covering an annual birth population of 482,000, or 10 % of U.S. births. OHD 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. A detailed description of NBDPS methods can be found elsewhere ( Gallagher et al. 2011 ; Rasmussen et al. 2002 ; Yoon et al. 2001 ). In the current study, we included all available mother-offspring pairs with estimated dates of delivery between October 1997 and August 2008. Case pairs were defined as those in which the child had at least one type of left or right OHDs, including hypoplastic left heart syndrome, tricuspid atresia, pulmonary valve atresia, coarctation of the aorta, interrupted aortic arch, aortic stenosis, valvar, pulmonic valve stenosis, and ebstein anomaly. 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 ( Yoon et al. 2001 ) and thus represent a random sample from the general population. The study population included 294 case pairs and 874 control pairs. A comparison of maternal characteristics ( Table 1 ) indicates that there was no significant difference between cases and controls with regard to maternal age, race, household income, education, smoking, drinking and folic acid exposure (all p values >0.05). However, the maternal BMI of case mothers was significantly higher than that of control mothers ( p value <0.001). In order to adjust for the possible confounding effect of BMI, mother-offspring pairs with missing maternal BMI values were excluded from the analysis (9 case pairs and 26 control pairs), and the final analytical dataset includes 285 case pairs and 868 control pairs. 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 described by Chowdhury et al. (2012) . We found that the quality of genotype clustering varied substantially from SNP to SNP, 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. The initial genotype calls, along with the raw intensity data, were used as inputs to SNPMClust, a bivariate Gaussian model-based genotype clustering and calling algorithm developed in-house, currently available as an R package on the Comprehensive R Archive Network (CRAN; http://cran.r-project.org/ ). 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 (60 SNPs), no-call rates >10 % (389 SNPs), Mendelian error rates >5 % (10 SNPs), minor allele frequencies (MAFs) <5 % (192 SNPs), or significant deviation from Hardy–Weinberg Equilibrium in at least one racial group ( p < 10 −4 , 8 SNPs). The haplotype blocks were determined using software Haploview version 4.2 ( Barrett et al. 2005 ). Linkage disequilibrium (LD) was first measured by the D ′ statistic between neighboring pairs of genetic variants. The Solid Spine of LD criterion, an internally developed method by Haploview, was used to determine the haplotype blocks using a threshold of D ′ > 0.6. The LD structure was ignored if two SNPs were located more than 500 kb apart from each other, which is the default setting. A total of 140 haplotype blocks were identified for association analysis. Recently, we and others proposed using penalized logistic regression methods to dissect the maternal, fetal, and MFG interaction effects with mother-offspring pair data ( Li et al. 2009 , 2010 ). In this article, we adapt our method to detect main effects and interaction effects between haplotypes, the theoretical details of which can be found elsewhere ( Li et al. 2010 ). The proposed method was implemented in R, and the source code is freely available upon request. We briefly describe our method below. Denote H as a risk haplotype that potentially alters the likelihood of disease. The multi-locus genotypes within the haplotype block can then be mapped into three possible composite diplotypes, namely HH , HH̄ and H̄H̄ , where H̄ represents all haplotypes that are different from the risk haplotype H . A logistic regression framework is then used to model the maternal effect, fetal effect, and MFG interaction effect: (1) logit ( p ( y i = 1 ) ) = μ + a m x i , m + a f x i , f + d m z i , m + d f z i , f + i aa x i , m x i , f + i ad x i , m z i , f + i da z i , m x i , f + i dd z i , m z i , f + γ 1 U W i + γ 2 O W i + γ 3 O B i where x i , m = { 1 for G i , m = H H 0 for G i , m = H H ¯ - 1 for G i , m = H ¯ H ¯ and z i , m = { - 1 / 2 for G i , m = H H 1 / 2 for G i , m = H H ¯ - 1 / 2 for G i , m = H ¯ H ¯ . and with x i ,f and z i ,f defined similarly. This coding strategy follows Cockerham’s orthogonal partition method ( Cockerham 1954 ; Kao and Zeng 2002 ), where a m(f) and d m(f) can be interpreted as the additive and dominance effects for the risk haplotype at a maternal (fetal) block, and i aa , i ad , i da , i dd can be interpreted as the additive × additive, additive × dominance, dominance × additive, and dominance × dominance interaction effects between the maternal and fetal blocks, respectively. UW, OW and OB are three dummy variables corresponding to three BMI categories: underweight, overweight, and obese. The maternal–fetal genotype at a haplotype block has seven possible combinations, each corresponding to a likelihood of disease ( Table 2 ), and coefficients were estimated by maximizing the penalized likelihood with an adaptive lasso penalty ( Zou 2006 ). During parameter estimation, phase-ambiguous genotypes were treated as missing data, and phase was determined probabilistically via an expectation–maximization (EM) algorithm. The initial probability of haplotypes was estimated using the R package haplo. stats, which also allows for missing genotypes ( Lake et al. 2003 ; Schaid et al. 2002 ). In our model, the haplotype block may have a large number of multi-locus genotypes. However, the number of composite diplotypes is always reduced to three after the haplotype configuration, which significantly lessens data dimensionality. Such a modeling strategy has also been widely adopted in previous studies ( Lin et al. 2007 ; Liu et al. 2004 , 2011 ; Zhang et al. 2012 ). Further, adaptive lasso will simultaneously estimate parameters and perform model selection through shrinkage. Coefficients that do not significantly differ from 0 are expected to be shrunk to 0. Standard errors and confidence intervals of non-zero coefficients are empirically estimated using bootstrap resampling ( Tibshirani 1996 ). It is worthwhile to note that a risk haplotype is defined here for the purpose of dimension reduction. A risk haplotype may actually have a protective effect that corresponds to a lower likelihood of disease. In the real data application, all haplotypes with a frequency greater than 5 % were examined as potential risk haplotypes, and the haplotype with a minimum BIC was selected as the optimal risk haplotype. It should be noted that both the adaptive lasso estimator and the BIC criteria are asymptotically consistent in terms of model selection, which means that the probability of the selected model being the true model converges to 1 as the sample size increases ( Yang 2005 ; Zou 2006 ). The Bayesian false-discovery probability (BFDP) was used to evaluate the significance of associations for each haplotype block using the estimated odds ratios (OR) and corresponding 95 % confidence intervals ( Wakefield 2007 ). BFDP has become a popular strategy to access the noteworthiness of an association by balancing the costs of false discovery and non-discovery. In our analyses, we used the most widely accepted practice of pre-setting the BFDP threshold at 0.80, ( Liu et al. 2010 ; Oh et al. 2010 ; Park et al. 2010 ; Spitz et al. 2012 ; Wakefield 2007 ; Zienolddiny et al. 2013 ) which means that missing a true association is considered four times more costly than falsely reporting an association. We also assume a 5 % prior probability for the existence of an association between OHD outcome and haplotype blocks. The upper bound of the 95 % CI of prior OR was assumed to be 1.5, meaning that the prior probability of an OR being greater than 1.5 is 2.5 %.