Gene-by-gene interactions associated with the risk of conotruncal heart defects.

BACKGROUND: The development of conotruncal heart defects (CTDs) involves a complex relationship among genetic variants and maternal lifestyle factors. In this article, we focused on the interactions between 13 candidate genes within folate, homocysteine, and transsulfuration pathways for potential association with CTD risk.
METHODS: Targeted sequencing was used for 328 case-parental triads enrolled in the National Birth Defects Prevention Study (NBDPS). To evaluate the interaction of two genes, we applied a conditional logistic regression model for all possible SNP pairs within two respective genes by contrasting the affected infants with their pseudo-controls. The findings were replicated in an independent sample of 86 NBDPS case-parental triads genotyped by DNA microarrays. The results of two studies were further integrated by a fixed-effect meta-analysis.
RESULTS: One SNP pair (i.e., rs4764267 and rs6556883) located in gene MGST1 and GLRX, respectively, was found to be associated with CTD risk after multiple testing adjustment using simpleM, a modified Bonferroni correction approach (nominal p-value of 4.62e-06; adjusted p-value of .04). Another SNP pair (i.e., rs11892646 and rs56219526) located in gene DNMT3A and MTRR, respectively, achieved marginal significance after multiple testing adjustment (adjusted p-value of .06).
CONCLUSION: Further studies with larger sample sizes are needed to confirm and elucidate these potential interactions.

INTRODUCTION

Congenital heart defects (CHDs) are the most common type of birth defects, affecting 48.9 million people globally in 2015 (Vos et al., 2016). The prevalence of disease is approximately 1% (Vos et al., 2015), contributing to ~300,000 infant deaths per year (Naghavi et al., 2015; Wang et al., 2016). Children with CHDs are usually at increased risks of other health conditions, such as developmental disabilities (Limperopoulos et al., 2000; Razzaghi, Oster, & Reefhuis, 2015) and cognitive disorders (Shillingford et al., 2008). CHDs comprise various cardiac anomalies that are etiologically heterogeneous. Conotruncal heart defects (CTDs), such as tetralogy of Fallot, truncus arteriosus, double outlet right ventricle, transposition of the great arteries, pulmonary atresia, malalignment ventricle septal defect, and interrupted aortic arch type B, are a major and severe subtype of CHDs, accounting for 20%–30% of all CHD cases (Botto, Lin, Riehle‐Colarusso, Malik, & Correa, 2007; Kuehl & Loffredo, 2005). With the advance of surgical therapy and medical care, the survival rate of CTD‐affected infants has improved substantially over the past 50 years (Oster et al., 2013). However, those infants who survived may require repeated surgeries and continuing medical care into their adulthood (Gilboa et al., 2016), imposing enormous emotional and financial burden to their families and the health care system (McClung, Glidewell, & Farr, 2018; Simeone et al., 2014). Understanding the pathogenesis of CTDs is of crucial importance to reduce the public health impact of the disease.

Previous studies have identified a number of risk factors, such as maternal obesity (Brite, Laughon, Troendle, & Mills, 2014; Gilboa et al., 2010), diabetes (Lisowski et al., 2010; Stavsky et al., 2017), tobacco use (Alverson, Strickland, Gilboa, & Correa, 2011; Malik et al., 2008), and medication use (Li et al., 2014a, 2014b) during pregnancy. Genetic variants have also been identified for association with disease susceptibility. For example, the microdeletion on the long arm of chromosome 22 (i.e., DiGeorge syndrome) substantially increases the disease risk, and may account for ~12% of conotruncal malformations (Hacihamdioglu, Hacihamdioglu, & Delil, 2015; Liu et al., 2014). Other chromosomal abnormalities and copy number variations, such as 1q21, 8p23, and11q25 (Edwards & Gelb, 2016; Thienpont et al., 2007), may also increase the risk of CTDs. In addition, previous studies have shown that multiple genes, such as NKX2.5, GATA4, GATA6, TBX1, TBX5, CITED2, HAND2, NOTCH1, JAG1, ACTC1, and THRAP2, are associated with the transcription factors, ligands‐receptors, contractile and miscellaneous proteins that are involved in essential cardiac developmental processes (Hu et al., 2013; Huang et al., 2010; Li, Pu, Liu, Xu, & Xu, 2017; Wessels & Willems, 2010; Zhang, Hong, et al., 2018a). In recent years, genetic association studies have identified additional genes from various biological pathways, such as folic acid, homocysteine, and transsulfuration, for potential association with CTD risk (Hobbs et al., 2014; Shaw et al., 2005, 2009). However, these genes generally have small to intermediate effect sizes, and collectively may only explain a small proportion of CTDs cases (Fahed, Gelb, Seidman, & Seidman, 2013). The etiology of CTDs is highly complex, involving multiple genetic variants, maternal environmental exposures and lifestyle factors (Gelb & Chung, 2014; Pierpont et al., 2007).

To date, most studies have adopted a single‐locus analysis strategy, by testing each single genetic variant for potential association with the disease. However, accumulating evidence has demonstrated that gene‐by‐gene interaction exists pervasively in biological pathways and may be a major source of the unexplained disease heritability (Eichler et al., 2010; Moore, 2003). For example, an animal study analyzed the transcriptome data among mouse embryos with various combinations of wild‐type and mutant genotypes. The results suggested that the Tbx5/Osr1/Pcsk6 gene interaction was associated with atrial septal defects through the signaling network (Zhang et al., 2016). Studies have also shown that gene‐by‐gene interactions can alter the risk of CTDs among human populations. For example, physical interactions were found between particular combinations of GATA4 and TBX5 molecules, with varying frequencies between families with and without cardiac septal defects (Garg et al., 2003). It was also suggested that additional gene‐by‐gene interactions might exist between MYH6 and TBX5 (Chen et al., 2017), MYH6 and NKX2‐5 (Granados‐Riveron et al., 2012), and ZFPM2 and GATA4 (Zhang et al., 2014). These studies might have provided potential explanations in terms of the complex pathophysiological and etiological processes that result in CTDs, and also support that gene‐by‐gene interactions are a common mechanism contributing to the disease development.

In this study, we evaluate the potential gene‐by‐gene interactions among 13 candidate genes selected from three biological pathways (i.e., folate, homocysteine, and transsulfuration) that have altered metabolites in CHD‐affected pregnancies (Hobbs, Cleves, Zhao, Melnyk, & James, 2005; Hobbs et al., 2006; Kapusta et al., 1999; Obermann‐Borst et al., 2011). A total of 515 common variants from 13 genes were analyzed in two independent samples, including a sequencing study with 328 case‐parental triads and a microarray study with 86 case‐parental triads. The analysis results from two independent samples were further integrated by a fixed‐effect meta‐analysis.

MATERIALS AND METHODS

Ethics statement

The study was approved by the Institutional Review Boards at University of Arkansas for Medical Sciences and Indiana University Bloomington, 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 consent. For minors, informed written consent was obtained from their legal guardian for DNA collection.

Study population

The National Birth Defects Prevention Study (NBDPS) is a large population‐based, case‐control study of birth defects. The detailed study design and procedure has been described elsewhere (Gallagher et al., 2011; Rasmussen et al., 2002; Yoon et al., 2001). Briefly, since 1997, subjects have been recruited from the population‐based birth defects surveillance system in ten states (i.e., AR, CA, GA, IA, MA, NC, NJ, NY, TX, and UT). The study covers an annual birth population of 482,000 (i.e., 10% of U.S. births). Case infants were identified if presenting at least 1 major birth defect. Only non‐syndromic cases were included, and cases with known single gene defects or chromosomal abnormalities were excluded. Case information and medical records were obtained by trained specialists. All diagnostic tests were reviewed by pediatric cardiologists to ensure uniform criteria for diagnoses. Study information was collected from recruited mothers via computer‐assisted telephone‐interviews with detailed questions about environmental exposures. After the interview, a buccal cell collection kit was sent to collect biological samples of both biological parents and the infant. The collected buccal cells were sent back to NBDPS laboratories for DNA extraction and genotyping.

In our current study, we adopt the case‐parental triad design. A case‐parental triad is defined as a singleton live‐born infant with CTDs and his/her biological parents. Only complete triads with genetic data available for three family members were included in our study. Subjects in the sequencing study was genotyped using Illumina's targeted sequencing technology, while those in the microarray study were genotyped in one of our earlier studies using Illumina® GoldenGate™ platform (Hobbs et al., 2014). The human reference genome is version GRCh37.p13.

Sequencing study

The sequencing study included 328 case‐parental triads. Each subject was sequenced for 13 targeted regions, including exons, introns, and 1kb upstream and downstream regions for the following folate‐related genes: GCLC (OMIM: 606857), SHMT1 (OMIM: 182144), NOS2A (OMIM: 163730), MTRR (OMIM: 602568), MTHFS (OMIM: 604197), GLRX (OMIM: 600443), DNMT3A (OMIM: 602769), SOD2 (OMIM: 147460), MGST1 (OMIM: 138330), GPX4 (OMIM: 138322), TCN2 (OMIM: 613441), MTHFD2 (OMIM: 604887), and TYMS (OMIM: 188350). These genes were selected for sequencing based on a preliminary analysis of data from an earlier study (Hobbs et al., 2014). The 13 selected genomic regions are described in Table 1. DNA was extracted in accordance with well‐established NBDPS protocols. Extracted DNA (50 ng per sample) was prepared for targeted sequencing using Illumina Nextera hybridization enrichment. Target libraries were pooled at 48x and sequenced on a HiScan‐SQ or NextSeq500.

Table 1

Candidate genes or regions sequenced in this studya

Genes	Region of SNPs	Total (kb)	Number of SNPs	Pathways
DNMT3A	chr2: 25,451,420–25,568,539	117.1	45	Homocysteine
MTHFD2	chr2: 74,423,130–74,449,302	26.2	10	Folate
MTRR	chr5: 7,846,305–7,908,359	62.0	46	Homocysteine
GLRX	chr5: 95,148,650–95,164,023	15.4	9	Transsulfuration
GCLC	chr6: 53,357,360–53,414,706	57.3	37	Transsulfuration
SOD2	chr6: 160,091,691–160,119,427	27.7	29	Transsulfuration
MGST1	chr12: 16,492,829–16,533,128	40.3	41	Transsulfuration
MTHFS	chr15: 80,130,515–80,198,784	68.3	78	Folate
SHMT1	chr17: 18,222,647–18,274,323	51.7	91	Folate
NOS2A	chr17: 26,075,524–26,131,931	56.4	31	Transsulfuration
TYMS	chr18: 655,999–680,838	24.8	34	Folate
GPX4	chr19: 1,100,976–1,109,213	8.2	10	Transsulfuration
TCN2	chr22: 31,001,072–31,029,898	28.8	54	Folate

Human reference genome version is GRCh37.p13.

Microarray study

The microarray study included 86 case‐parental triads that are non‐overlapping subjects with the sequencing study. These subjects are NBDPS participants in one of our previous studies (Hobbs et al., 2014). Each subject was initially genotyped for about 1,500 genetic variants in 62 candidate genes, including the 13 genes that were included in the sequencing study. Illumina's GoldenGate Custom DNA microarray was used for genotyping. The detailed process of SNP selection, DNA extraction and genotyping can be found elsewhere (Hobbs et al., 2014).

Quality control

Thirty‐one samples from the sequencing study were also genotyped by Illumina HumanOmni 5M BeadChip as part of another project (unpublished data). Genotype calls were compared at 196 overlapping sites, and the concordance rate across platforms was 99.0%, indicating high accuracy of variant calls. All samples were screened for Mendelian inconsistencies using PLINK (Purcell et al., 2007). A cell‐line from the 1,000 Genome project (Siva, 2008) was sequenced as a positive control and was used as an external quality assurance (EQA) measure for NBDPS. In the sequencing study, all SNPs had a call rate greater than 95%, and all families and SNPs showed Mendelian error rates less than 5%. All genotypes with Mendelian inconsistency were set to missing.

Imputation and SNP selection

The targeted sequencing provided genotypes for 9,402 variants in the sequencing study. However, only 74 of those variants were readily available in the microarray study. In order to maximize the overlap of variants between two studies, we further conducted genotype imputation for each study by using 1,000 Genome Project Phase 3 Genome Build b37 as reference panels. Software IMPUTE version 2 was used for the imputation with a “multi‐population reference panels” option (Howie, Donnelly, & Marchini, 2009), which chooses a “custom” reference panel for each individual based on all available reference haplotypes. The advantages for such a strategy has been discussed elsewhere (Howie, Marchini, & Stephens, 2011). After the imputation and quality control, the sequencing study and the microarray study included 9,588 and 3,505 variants, respectively. A total of 2,945 variants overlapped between two studies. In this article, we only considered 515 common variants with at least 5% minor allele frequencies in both studies for interaction analysis. We excluded variants with an allele frequency of less than 5% due to the limited statistical power for detecting potential interactions among rare variants.

Statistical methods

Single‐locus association test

We first evaluated the association between each genetic variant and CTD risk with a genotypic transmission disequilibrium test (TDT) (Schaid, 1996, 1999). Conditional logistic regression model was fitted assuming additive inheritance. The conditional probability of disease risk was modeled by contrasting affected case infants with the corresponding pseudo‐controls whose allele combinations are theoretically possible but not observed based on Mendelian transmission of given parents’ genotypes (Falk & Rubinstein, 1987). A major advantage of using matched pseudo‐controls is its robustness to confounding causes of association, such as population stratification (Cordell, 2002). Such a strategy has been commonly used for studies of other birth defects, such as cleft palate (Liu et al., 2018, 2019; Xiao et al., 2017). We used genotypic TDT instead of conventional allelic TDT to have consistent modeling with the interaction analysis described below. The genotypic TDT analysis was performed in R version 3.43 using Bioconductor package “trio” (Schwender et al., 2014).

Gene‐by‐gene interaction

Our study included 13 candidate genes (Table 1), leading to a total of 78 gene‐by‐gene combinations. For each given gene‐by‐gene combination, we estimated the interaction effect for all possible SNP pairs within two respective genes. We denote as the conditional probability of being a case infant rather than the pseudo‐controls for the i‐th family; and and to be the genotype of two SNPs from two respective genes, coded as minor allele counts. A conditional logistic regression model was fitted as follows:

The model conducted matched case‐control analysis by contrasting the affected infants with all possible pseudo‐controls assuming additive inheritance. For single bi‐allelic variant, two heterozygous parents can produce four inheritance patterns of genotypes, suggesting 1 case versus 3 matched pseudo‐controls. Similarly, two bi‐allelic variants lead to 16 two‐SNP inheritance patterns. The ratio between cases and their pseudo‐controls is 1:15. The interaction effect between two variants can thus be tested against the null hypothesis of 0. The same model was fitted for both the sequencing study and the microarray study. The analysis was conducted in R version 3.43 using Bioconductor package “trio” (Schwender et al., 2014).

Meta‐analysis

We further conducted a fixed‐effect meta‐analysis in order to integrate the results from the sequencing and the microarray studies. Top associations were selected based on the significance of meta‐analysis (i.e., p‐values) and the consistency of testing results across two studies, that is, achieving nominal significance level in both studies and having the same direction of effect. The analysis was conducted in R version 3.43 using package “rmeta.”

Multiple testing adjustment

A total of 118,757 SNP pairs were tested for possible interaction effects. Because of the strong linkage‐disequilibrium (LD) structure within the same candidate genes (Fig. S1), the conventional Bonferroni correction would be overly conservative. Alternatively, we adopted a SimpleM method for estimating the effective number of independent tests among highly correlated data (Gao, Starmer, & Martin, 2008) and applied a Bonferroni correction based upon this estimate. The details and advantages of SimpleM has been discussed elsewhere (Gao, Becker, Becker, Starmer, & Province, 2010). In particular, SimpleM uses principal component analysis to estimate the effective number of independent tests while accounting for the majority of genetic variation (e.g., 95%). All analyses were performed using R version 3.4.3.

RESULTS

The demographic information of our study populations is summarized in Table 2. A total of 328 and 86 case‐parent triads were included in the sequencing and the microarray studies, respectively. Most of the maternal characteristics, including race, body mass index, education level, folic acid supplementation, alcohol use and smoking status, were similar between two studies (i.e., p‐values > .05). The average maternal age at delivery was 28.31 and 30.03 years in the sequencing and the microarray studies, respectively, which were statistically different but not of clinical importance. Family income levels were also statistically different between two studies, which may largely due to the high level of missing values in the sequencing study. In both studies, the most frequent maternal characteristics were Caucasian race, normal body mass index (BMI), and over 16 years of education. The majority of the families had more than $50,000 annual household income. In terms of maternal lifestyle factors, more women were taking folic acid supplementation regularly during the pregnancy, and less women smoking or drinking.

Table 2

Characteristics of case‐parental triads in sequencing and microarray study

	Sequencing study (n = 328)	Microarray study (n = 86)	p‐value
Maternal Age
Mean (SD)	28.31 (6.00)	30.03 (6.21)	.022
Maternal race, n (%)			.604
Caucasian	233 (71.04%)	68 (79.07%)
African American	19 (5.80%)	3 (3.49%)
Hispanic	52 (15.85%)	10 (11.63%)
Others	24 (7.32%)	5 (5.81%)
Maternal BMI, n (%)			.161
Underweight	14 (4.27%)	2 (2.32%)
Normal	147 (44.82%)	53 (61.73%)
Overweight	81(24.70%)	16 (18.60%)
Obese	76 (23.17%)	15 (17.44%)
Missing	10 (3.05%)	0 (0.00%)
Maternal Education, n (%)			.466
0–11 years	35 (10.67%)	7 (8.14%)
12 years	72 (16.77%)	19 (22.09%)
13–15 years	100 (30.49%)	21 (24.42%)
16 or more years	121 (36.89%)	36 (41.86%)
Missing	0 (0.00%)	3 (3.49%)
Family Income, n (%)			.004
< 10,000	34 (10.36%)	9 (10.46%)
10,000 to 30,000	55 (14.78%)	15 (17.44%)
30,000 to 50,000	41 (12.50%)	13 (15.12%)
> 50,000	62 (18.90%)	46 (53.49%)
Missing	136 (41.46%)	3 (3.49%)
Folic Acid Supplementation, n (%)	193 (58.84%)	49 (56.98%)	.755
Maternal Alcohol Use, n (%)	92 (28.05%)	33 (38.37%)	.063
Maternal Smoke, n (%)	44 (13.41%)	16 (18.60%)	.211
Missing	1 (0.30%)	1 (1.16%)

Single‐locus association test

Each of 515 SNPs was examined using a genotypic TDT as described above. The testing p‐values are shown in Figure 1. The multiple testing significance thresholds were calculated as 0.05/144 = 3.5*10−4, where the effective number of independent tests estimated by SimpleM was 144. In addition, results for each gene were summarized in Table 3. None of the associations remained significant after multiple testing adjustment.

Figure 1

Barplot of single‐locus genotypic TDT test from meta‐analysis. No significant findings after multiple testing adjustment. Each bar represents the negative 10‐based logarithm of p‐value for a SNP, ordered by its physical location in the genome. The dashed line represents the threshold of statistical significance after multiple testing adjustment

Table 3

Results of single‐locus genotypic TDT test for each genea

Gene	SNP	Chro	Position	Study	Estimate (SE)	95% CI	p‐value	Adjusted p‐value
MTHFD2	rs1723285	2	74,426,038	Sequencing Microarray Meta‐analysis	1.42 (0.15) 1.66 (0.22) 1.49 (0.13)	(1.06, 1.92) (1.07, 2.57) (1.16, 1.91)	2.04*10−2 2.42*10−2 1.45*10−3	.21
NOS2A	rs2072324	17	26,116,896	Sequencing Microarray Meta‐analysis	0.65 (0.49) 0.76 (0.26) 0.68 (0.12)	(0.49, 0.86) (0.46, 1.27) (0.53, 0.86)	2.47*10−4 0.30 1.63*10−3	.23
MGST1	rs10744119	12	16,522,636	Sequencing Microarray Meta‐analysis	0.68 (0.14) 0.90 (0.26) 0.72 (0.12)	(0.52, 0.90) (0.54, 1.51) (0.57, 0.92)	6.08*10−3 0.69 9.06*10−3	.99
MTRR	rs1801394	5	7,870,973	Sequencing Microarray Meta‐analysis	1.40 (0.12) 0.98 (0.22) 1.29 (0.10)	(1.11, 1.77) (0.63, 1.50) (1.05, 1.58)	4.38*10−3 0.91 1.41*10−2	1.00
SOD2	rs5746105	6	160,112,638	Sequencing Microarray Meta‐analysis	0.66 (0.12) 1.23 (0.23) 0.76 (0.11)	(0.52, 0.85) (0.79, 1.92) (0.62, 0.95)	1.05*10−3 0.36 1.47*10−2	1.00
TCN2	rs4820886	22	31,016,539	Sequencing Microarray Meta‐analysis	0.81 (0.18) 0.41 (0.34) 0.70 (0.16)	(0.57, 1.14) (0.21, 0.81) (0.52, 0.96)	0.22 1.01*10−2 2.43*10−2	1.00
SHMT1	rs8073885	17	18,273,850	Sequencing Microarray Meta‐analysis	0.67 (0.22) 0.69 (0.39) 0.68 (0.20)	(0.43, 1.05) (0.32, 1.48) (0.46, 0.99)	7.92*10−2 0.34 4.56*10−2	1.00
MTHFS	rs35919462	15	80,131,692	Sequencing Microarray Meta‐analysis	1.43 (0.16) 0.90 (0.26) 1.26 (0.14)	(1.05, 1.94) (0.54, 1.50) (0.97, 1.64)	2.47*10−2 0.70 8.60*10−2	1.00
GCLC	rs2300420	6	53,376,551	Sequencing Microarray Meta‐analysis	0.66 (0.20) 1.17 (0.39) 0.74 (0.18)	(0.44, 0.98) (0.54, 2.52) (0.52, 1.06)	3.81*10−2 0.70 9.66*10−2	1.00
TYMS	rs2741184	18	679,660	Sequencing Microarray Meta‐analysis	1.36 (0.18) 1.08 (0.27) 1.27 (0.15)	(0.96, 1.92) (0.63, 1.84) (0.95, 1.69)	8.29*10−2 0.78 0.11	1.00
DNMT3A	rs62131064	2	25,568,528	Sequencing Microarray Meta‐analysis	0.60 (0.25) 1.22 (0.45) 0.71 (0.22)	(0.37, 0.98) (0.51, 2.95) (0.46, 1.09)	4.28*10−2 0.66 0.12	1.00
GLRX	rs871775	5	95,158,768	Sequencing Microarray Meta‐analysis	1.02 (0.19) 2.58 (0.34) 1.28 (0.17)	(0.70, 1.49) (1.33, 5.03) (0.92, 1.78)	0.92 5.25*10−3 0.14	1.00
GPX4	rs1808194	19	1,102,323	Sequencing Microarray Meta‐analysis	0.92 (0.12) 1.00 (0.26) 0.94 (0.11)	(0.73, 1.18) (0.60, 1.67) (0.75, 1.17)	0.54 1.00 0.58	1.00

Estimates are exponential of β coefficients, representing relative risks. SNPs were ordered based on adjusted p values.

Gene‐by‐gene interaction

Each SNP pair was examined using a conditional logistic regression model described above. Under the null hypothesis of no association, the testing p‐values were expected to follow a uniform distribution, and the distribution was evaluated by a quantile‐quantile (QQ) plot (Fig. S2). No early departure from the uniform distribution was found, with genomic inflation factors (λ) close to 1, suggesting no inflated type I error rate. The top 5 SNP pairs are summarized in Table 4. In particular, the interaction effect between rs4764267 and rs6556883, located in genes MGST1 and GLRX, respectively, had an adjusted p‐value of .04, which remained significant after multiple testing adjustment. The interaction effect and corresponding genotypic relative risks varied largely between two studies, which we attributed mainly to the small sample size of our microarray study. These two SNPs may form 9 possible two‐SNP genotype combinations. We estimated the relative risks among these genotype combinations formed by rs4764267 and rs6556883 (Table 5 ). To illustrate the interaction effect of these two SNPs, we further plotted these estimated relative risks in Figure 2. When genotype GG was observed at SNP rs6556883 (i.e. red line), allele T at SNP rs4764267 was estimated to increase the disease risk compared to allele G. However, when genotype AG or AA was observed at SNP rs6556883 (i.e. green and blue lines), allele T at SNP rs4764267 was expected to decrease the disease risk but with varying effect size. Such a pattern remained consistent across the two studies. SNPs rs4764267 and rs6556883 were located within genes MGST1 and GLRX, respectively. The other four SNP pairs in Table 4 were from the same gene pairs of DNMT3A and MTRR, and were in strong LD. We believe these four SNP pairs were representing the same gene‐by‐gene interaction. In particular, the SNP pair of rs11892646 and rs56219526 was marginally significant after multiple testing adjustments (p‐value = .06). We presented the analysis results in Table 6, and illustrated the relative risk among genotype combinations in Figure 3.

Table 4

Top 5 gene‐by‐gene interaction by conditional logistic regressiona

Gene Pair	SNP pair	Chro	Position	Study	Estimate (SE)	95% CI	p‐value	Adjusted p‐valueb
MGST1 GLRX	rs4764267 rs6556883	12 5	16,523,580 95,152,085	Sequencing Microarray Meta‐analysis	0.49 (0.18) 0.24 (0.50) 0.45 (0.17)	(0.34, 0.70) (0.09, 0.65) (0.32, 0.63)	1.23*10−4 4.78*10−3 4.62*10−6	.04**
MTRR DNMT3A	rs56219526 rs11892646	5 2	7,891,210 25,494,474	Sequencing Microarray Meta‐analysis	2.66 (0.27) 4.53 (0.52) 2.98 (0.24)	(1.56, 4.53) (1.62, 12.69) (1.85, 4.78)	3.34*10−4 4.03*10−3 6.55*10−6	.06*
MTRR DNMT3A	rs3776454 rs11892646	5 2	7,896,604 25,494,474	Sequencing Microarray Meta‐analysis	2.53 (0.27) 4.53 (0.52) 2.85 (0.24)	(1.49, 4.29) (1.62, 12.69) (1.78, 4.57)	5.96*10−4 4.03*10−3 1.25*10−5	.12
MTRR DNMT3A	rs326125 rs11892646	5 2	7,888,001 25,494,474	Sequencing Microarray Meta‐analysis	2.37 (0.27) 4.50 (0.50) 2.74 (0.24)	(1.40, 4.03) (1.68, 12.05) (1.72, 4.37)	1.39*10−3 2.76*10−3 2.31*10−5	.22
MTRR DNMT3A	rs10380 rs11892646	5 2	7,897,191 25,494,474	Sequencing Microarray Meta‐analysis	2.26 (0.24) 4.34 (0.52) 2.54 (0.22)	(1.40, 3.66) (1.56, 12.06) (1.65, 3.93)	8.98*10−4 4.87*10−3 2.62*10−5	.24

Estimates are exponential of β coefficients, representing relative risks.

Adjusted p‐values were marked with **if less than .05 and marked with *if less than .1.

Table 5

Relative risk for all possible genotype combinations of rs4764267 and rs6556883a

rs4764267	rs6556883	Estimate (Sequencing)	Estimate (Microarray)
GG	GG	1.00 (ref)	1.00 (ref)
TG	GG	1.70 (1.23, 2.36)	1.73 (0.77, 3.91)
TT	GG	2.89 (1.51, 5.56)	3.01 (0.59, 15.28)
GG	AG	1.55 (1.04, 2.31)	4.38 (1.59, 12.08)
TG	AG	1.29 (0.83, 2.02)	1.84 (0.54, 6.24)
TT	AG	1.08 (0.54, 2.18)	0.78 (0.11, 5.40)
GG	AA	2.41 (1.08, 5.35)	19.19 (2.53, 145.84)
TG	AA	0.98 (0.49, 2.00)	1.96 (0.27, 14.20)
TT	AA	0.40 (0.14, 1.19)	0.20 (0.01, 4.67)

Estimates are relative risks and 95% confident intervals.

Figure 2

Estimated relative risk of CTD among all genotypic combinations of rs4764267 and rs6556883: (a) is for sequencing study and (b) is for microarray study

Table 6

Relative risk for all possible genotype combinations of rs11892646 and rs56219526a

rs11892646	rs56219526	Estimate (Sequencing)	Estimate (Microarray)
TT	GG	1.00 (ref)	1.00 (ref)
CT	GG	0.78 (0.52, 1.16)	0.52 (0.21, 1.29)
CC	GG	0.61 (0.27, 1.35)	0.27 (0.04, 1.66)
TT	AG	0.48 (0.32, 0.73)	0.79 (0.38, 1.65)
CT	AG	1.00 (0.60, 1.64)	1.86 (0.71, 4.85)
CC	AG	2.07 (0.82, 5.20)	4.35 (0.76, 24.80)
TT	AA	0.23 (0.10, 0.54)	0.63 (0.14, 2.74)
CT	AA	1.28 (0.53, 3.09)	6.69 (1.20, 37.34)
CC	AA	7.02 (1.34, 36.71)	70.95 (3.11, 1619.84)

Estimates are relative risks and 95% confident intervals.

Figure 3

Estimated relative risk of CTD among all genotypic combinations of rs11892646 and rs56219526: (a) is for sequencing study and (b) is for microarray study

DISCUSSION

In this study, we conducted a two‐phase investigation to evaluate the possible interactions among 13 candidate genes for association with CTD risk. One SNP pair (i.e., rs4764267 and rs6556883) was identified with a significant p‐value after multiple testing adjustment, and another SNP pair (i.e., rs11892646 and rs56219526) was marginally significant. All identified genes are functionally involved in the transsulfuration pathway, supporting its possible contribution to the genetic susceptibility of CTDs.

Our results are consistent with previous work. Gene MGST1 encodes microsomal glutathione S‐transferase 1, which contributes to the antioxidant system by catalyzing glutathione binding to toxic electrophilic compounds. Decreased expression of MGST1 has been linked structural cardiac abnormalities in Nos3 deficient mice (Campbell, Li, Biendarra, Terzic, & Nelson, 2015). GLRX encodes glutaredoxin‐1, which catalyzes the glutathione regeneration and provides reactant for MGST1, and is also involved in the response to oxidative stress. Overexpression of GLRX2 in transgenic mice has been found to attenuate doxorubicin (DOX)‐induced heart defects (Diotte, Xiong, Gao, Chua, & Ho, 2009). Others and we have previously reported the association of polymorphisms in MGST1 and GLRX with CHDs in other NBDPS samples although the molecular mechanism remains unclear (Chowdhury et al., 2012; Nembhard et al., 2017, 2018). Meanwhile, accumulated evidence suggests that the oxidative stress/transsulfuration pathway plays an essential role in cardiac pathogenesis. For example, others and we have shown differences in markers of oxidative stress, including glutathione level, among women with infants affected by CHDs compared to those with healthy infants (Hobbs et al., 2005). In addition, decreased activity of antioxidant enzyme such as glutathione peroxidase and lower levels of molecular antioxidants were observed among children affected by CHDs compared to healthy controls (Mukhopadhyay, Gongopadhyay, Rani, Gavel, & Mishra, 2015; Rokicki, Strzalkowski, Klapcinska, Danch, & Sobczak, 2003). Moreover, animal studies showed that oxidative stress might explain the increased risk of CHDs among female mice with diabetes (Wang, Reece, & Yang, 2015; Wu et al., 2016). In terms of the genes DNMT3A and MTRR, evidence suggest their association with CHD, respectively, in both animal models (Deng, Elmore, Lawrance, Matthews, & Rozen, 2008; Feng et al., 2013) and population‐based case‐control studies (van Beynum et al., 2006; Shaw et al., 2009; Sheng et al., 2013). It is also biologically plausible that MTRR and DNMT3A may interact functionally since they work sequentially in homocysteine‐methionine cycle (Ly, Hoyt, Crowell, & Kim, 2012).

In our study, a few assumptions were made while evaluating the possible interactions. First, by contrasting the affected infants with pseudo‐controls, we assumed that these theoretically possible but unobserved siblings were healthy individuals. Such an assumption is in general reasonable since CTDs are a relative rare disease with around 0.1% prevalence in the general population (Zhang, Li, et al., 2018b). Second, additive mode of inheritance was assumed when estimating the genetic effects, which is commonly used in genetic association studies.

The current analysis has several strengths. First, it is robust to population stratification and other environmental confounders as a result of case‐parent design and use of conditional logistic regression model. Second, it integrates data from two studies to replicate results and achieve greater power. Third, we have optimized statistical power by including only common variants in the analysis. The study is limited by modest sample size, which reduces statistical power for epistasis detection. However, this is the one of the largest studies to investigate the gene‐by‐gene interaction and the risk of conotruncal heart defects.

CONCLUSION

We identified two gene‐by‐gene interactions that are potentially associated with conotruncal heart defects. The genes are involved in transsulfuration pathways. Further studies with larger sample sizes are needed to confirm the associations. While we provide a possible explanation of their interactions, these explanations are our speculations based on existing knowledge. Additional molecular studies are needed to elucidate the molecular mechanisms of this hypothesis.

CONFLICT OF INTEREST

The authors confirm that no conflicts of interest to disclose.

AUTHOR CONTRIBUTION

CL and ML conceptualized the study, conducted the data analysis, and wrote the manuscript. DW, CH, and SM designed the sequencing study. CH and SM designed the microarray study. CH, SM, and NBDPS collected the samples and provided study oversight. All authors read and revised the manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.