Maternal effect genes as risk factors for congenital heart defects.

Maternal effect genes (MEGs) encode factors (e.g., RNA) in the oocyte that control embryonic development prior to activation of the embryonic genome. Over 80 mammalian MEGs have been identified, including several that have been associated with phenotypes in humans. Maternal variation in MEGs is associated with a range of adverse outcomes, which, in humans, include hydatidiform moles, zygotic cleavage failure, and offspring with multi-locus imprinting disorders. In addition, data from both animal models and humans suggest that the MEGs may be associated with structural birth defects such as congenital heart defects (CHDs). To further investigate the association between MEGs and CHDs, we conducted gene-level and gene-set analyses of known mammalian MEGs (n = 82) and two common groups of CHDs: conotruncal heart defects and left ventricular outflow tract defects. We identified 14 candidate CHD-related MEGs. These 14 MEGs include three (CDC20, KHDC3L, and TRIP13) of the 11 known human MEGs, as well as one (DNMT3A) of the eight MEGs that have been associated with structural birth defects in animal models. Our analyses add to the growing evidence that MEGs are associated with structural birth defects, in particular CHDs. Given the large proportion of individuals with structural birth defects for whom etiology of their condition is unknown, further investigations of MEGs as potential risk factors for structural birth defects are strongly warranted.

Introduction

Maternal effect genes (MEGs) encode factors (e.g., RNA) that are present in the oocyte and required for embryonic development. As the embryonic genome is initially silent, these maternal factors control the earliest developmental events, including cell division, epigenetic remodeling, imprinting, and determination of higher-order chromatin structure. In addition, MEGs encode factors that drive the transition from maternal to embryonic genomic control of development.

The importance of MEGs in the development of non-mammalian species has been established.2, 3, 4, 5 Although less is known about the role of MEGs in the development of mammals, the first mammalian MEG was reported in 2000, and over 80 mammalian MEGs have subsequently been identified. Most evidence implicating these genes as MEGs comes from animal studies. However, the first human MEG was reported in 2006, and potentially damaging variants in at least 10 additional genes have been identified in women with a range of adverse post-fertilization reproductive outcomes, including zygotic cleavage failure, recurrent pregnancy loss, and hydatidiform moles.9, 10, 11, 12 The known mammalian MEGs are enriched for genes that map to biological pathways and processes (e.g., chromatin, imprinting, and methylation) that are essential for early embryogenesis and set the stage for subsequent developmental events.

Pathogenic variants in MEGs and absence of MEG products are often associated with early lethality (e.g., zygotic cleavage failure, pre-implantation developmental arrest). However, later phenotypes have also been observed. For example, in mouse models, absence of the maternal Ezh2 gene product is associated with severe growth retardation in neonates, and, in humans, pathogenic variants in several MEGs (NLRP2, NLRP5, NLRP7, and PADI6) have been identified in women with children affected by multi-locus imprinting disorders. Further, in humans, pathogenic variants in several MEGs are associated with a spectrum of outcomes with variability observed both within and between women. Understanding of the phenotypic spectrum associated with individual MEGs is, however, likely to be incomplete, given the small number of women carrying damaging MEG genotypes that have been described in the literature.

It has been suggested, based on theoretic arguments, that the phenotypic consequence of MEGs could include structural birth defects., In particular, a role for MEGs in defects of right-left asymmetry, which include defects of heart development, has been hypothesized based on the observation that chirality in snails is under the control of the maternal genome. Although studies of MEGs have not routinely included the evaluation of structural birth defects, eight MEGs (Dnmt1, Dnmt3a, Dnmt3l, Kdm1b, Nlrp2, Tet3, Trim28, and Zfp57; approximately 10% of the known mammalian MEGs), have been associated with birth defects, including congenital heart defects (CHDs), in animal models. In addition, in our genome-wide analysis to identify maternal genes associated with conotruncal heart defects, several of the most significant associations involved genes that have been implicated as MEGs. Based on this observation, we compiled a list of MEGs from review articles published between 2014 and 2016,,, and conducted a post hoc analysis to determine whether associations with p < 0.05 in our analyses were enriched for the genes on this list. This analysis identified a 2.3-fold enrichment of MEGs (enrichment p = 0.057), but was driven by the gene-level study results (i.e., it was not planned a priori) and used an outdated list of MEGs. Consequently, we completed a comprehensive review of the literature and systematically curated an up-to-date list of mammalian MEGs. Building on this work, we then undertook a new study to assess the hypothesis that mammalian MEGs (either individually or as a set) are associated with CHDs.

Subjects and methods

Datasets

Our analyses were based on four CHD case-parent trio datasets derived from study populations recruited through the Children’s Hospital of Philadelphia (CHOP) and the Pediatric Cardiac Genomics Consortium (PCGC). Study participants were recruited under protocols approved by the institutional review boards at CHOP or the PCGC clinical study sites. CHD cases included individuals with a conotruncal defect (CTD) or left ventricular outflow tract defect (LVOTD). The four datasets are referred to as CHOP-CTD, CHOP-LVOTD, PCGC-CTD, and PCGC-LVOTD.

Cases with CTDs included individuals with tetralogy of Fallot, D-transposition of the great arteries, ventricular septal defects (conoventricular, posterior malalignment, and conoseptal hypoplasia), double-outlet right ventricle, isolated aortic arch anomalies, truncus arteriosus, or interrupted aortic arch. Cases with LVOTDs included individuals with hypoplastic left heart syndrome, coarctation of the aorta with or without bicuspid aortic valve, and aortic valve stenosis and excluded individuals with variants of hypoplastic left heart syndrome, such as mal-aligned atrioventricular canal defects or double-outlet right ventricle with mitral valve atresia. Medical records were reviewed to ensure accuracy of the cardiac phenotype. Trios in which the case had a known or suspected genetic syndrome were excluded.

Genetic data

Genotype data were generated using Illumina arrays: CHOP-CTD, HumanHap550 v1, v3, and 610 BeadChip; CHOP-LVOTD, HumanHap550 v3; and PCGC-CTD and LVOTD, 1M and 2.5M arrays. For the current analyses, additional genotypes were imputed using the Michigan Imputation Server with the TOPMed reference panel build 38 (version R2), which is based on high-coverage, whole-genome sequence data.20, 21, 22 Prior to imputation, the following exclusions were applied to each dataset: case-parent trios with a Mendelian error rate >1%; duplicate samples (i.e., samples with pairwise identity by descent >0.6); individuals with genotyping rates <95%; non-autosomal variants; and variants with a minor allele frequency (MAF) <1%, genotyping rate <90%, or deviation from Hardy Weinberg equilibrium (p ≤ 1 × 10−5) in parents.

The two CHOP datasets were combined and imputed together. Similarly, the two PCGC datasets were combined and imputed together. For both imputations, we applied all of the standard Michigan Imputation Server variant-level, chunk-level, and sample-level quality control measures, and used the Eagle v.2.4 algorithm and the Haplotype Reference Consortium panel to estimate haplotype phase. For the current analyses, variants with quality scores ≤0.8 were also excluded.

Statistical analyses

Characteristics of the study participants and the genetic data used in these analyses were summarized using counts and proportions.

Gene-level analyses

Maternal genetic effects were evaluated using an unmatched case-control approach in which mothers and fathers from the case-parent trios were considered as cases and controls, respectively. Genes were defined by their transcription start-stop coordinates, based on the Genome Reference Consortium Human Genome Build 38. For each gene, we specified an annotation window that included 1 kb up- and downstream of the start-stop coordinates, and all genotyped and imputed variants in the window were mapped to the gene. The four datasets (CHOP-CTD, CHOP-LVOTD, PCGC-CTD, and PCGC-LVOTD) were analyzed individually and results from the individual datasets were combined in meta-analyses. Since the effect of any given MEG may be specific to a particular type of CHD or more broadly associated with CHDs, meta-analyses were conducted for the CTD datasets only, the LVOTD datasets only, and both CTD and LVOTD datasets. Analyses of the individual datasets and the meta-analyses were conducted using Multi-marker Analysis of Genomic Annotation (MAGMA), a stand-alone program for gene and gene-set analyses.

In MAGMA, gene-level p values can be estimated using a multiple linear principal components regression model, in which case-control status is regressed on principal components derived from the variants in a gene and an F-test is used to compute the gene-level p value. Gene-level p values can also be estimated from variant-level p values that are calculated within MAGMA. These gene-level p values can be estimated by the average p value across all variants in the gene or by the lowest variant-level, permutation-based p value within the gene. As implemented in MAGMA, all three approaches have been shown to have appropriate type-I error rates. In addition, aggregate gene-level p values can be obtained by combining two or more of the gene-level estimates. We used an aggregate statistic based on all three estimates, because it provides a more even distribution of power and sensitivity for a wider range of genetic models than the other gene-level statistics. In MAGMA, the aggregate statistics are transformed to Z scores using the probit function, and the resulting scores are used as input for gene-level meta-analyses and gene-set analyses.

For each dataset, we conducted gene-level analyses using only data for common variants (MAF ≥ 0.05), only data for rare variants (MAF < 0.05), and data for both common and rare variants. For each gene, rare variants were included using a single burden score, with the included variants weighted by their allele frequency. To control for population stratification (i.e., differences in the genetic ancestry of cases [mothers] and controls [fathers]), all models were adjusted for genotypic principal components. For each dataset, principal component analyses were conducted using PLINK (version 1.09) and data from mothers and fathers of cases, and included only variants with an MAF ≥5% and linkage disequilibrium R2 ≤ 0.20. Pairwise combinations of the first 20 principal components were plotted in ascending order (e.g., the first plot was between the first and the second principal components, the second plot between the third and fourth). Based on visual inspection of these plots, principal components from pairs that showed structure (i.e., clustering versus random distribution) were included in the models. For each dataset, observations that were ≥6 standard deviations from the mean of any principal component that was included in the models were excluded from all analyses. Scatterplots were constructed and outliers were identified using R (version 3.4). Because a large proportion of the variants included in these analyses are likely to be benign, all gene-level analyses were repeated using only variants predicted to fall in the top 10% of deleterious substitutions. This subset of variants was identified using combined annotation dependent depletion (CADD) scores, which combine multiple annotations (e.g., conservation metrics, regulatory information) into a single deleteriousness score., Only variants with CADD scores (version 1.6) of at least 10 (corresponding to the top 10% most deleterious substitutions) were included.

Meta-analyses were conducted to combine data from the two CTD datasets, the two LVOTD datasets, and all four datasets. These meta-analyses were conducted with MAGMA using Stouffer’s Z method to combine gene-level Z scores across datasets, with weights set to the square root of the sample size. Meta-analyses were conducted using results from the common variant only, rare variant only, and common plus rare variant analyses based on all variants (referred to as unfiltered analyses) and based on only variants with CADD ≥10 (referred to as CADD10 filtered analyses). To aid in the interpretation of results across analyses, we calculated the Pearson correlation coefficient for the gene-level meta-p values across analyses (e.g., common variant only versus common + rare; unfiltered versus CADD10 filtered).

Although our analyses were focused on MEGs, the gene-set analyses (described below) required gene-level p values for both MEGs and other genes. Hence, for completeness, we assessed the association of each gene with case/control status. For these genome-wide gene-level analyses, an association was considered significant when the gene-level meta-p-value was less than the Bonferroni p value corrected for the total number of evaluated genes, and was considered suggestive when meta-p < 10−5. Genes for which the meta-p value was at least suggestive and less than the p values for each contributing dataset (i.e., indicating that the evidence for association was stronger in the combined data than in each of the contributing datasets) were considered candidate maternal CHD-related genes. As we have previously conducted a maternal gene-level meta-analysis of the CHOP-CTD and PCGC-CTD datasets using a different analytic approach, results from the current CTD meta-analysis were compared with those from our prior study. Details of our prior analytic approach and comparison of results from the two approaches are provided in the Supplemental data and Table S1.

For the MEG-specific analyses, we considered the list of 82 mammalian MEGs reported by Mitchell (Table S2). An association with an MEG was considered to be significant when the gene-level meta-p value was less than the Bonferroni p value corrected for the number of evaluated MEGs, and was considered suggestive when meta-p < 0.05. Genes for which the meta-p value was at least suggestive and less than the p values for each contributing dataset were considered candidate CHD-related MEGs.

MEG gene-set analyses

We used MAGMA to conduct competitive MEG-gene-set analyses using linear regression. Competitive gene-set analyses test whether genes in a set (e.g., MEGs) are more strongly associated with the outcome (e.g., CHDs) than are genes that are not in the set. The dependent variable in these analyses was the gene-level Z score and the primary independent variable was a binary variable (S) indicating whether a gene is (S = 1) or is not (S = 0) in the MEG gene set. Additional covariates were included, using the default options in MAGMA, to control for gene size, mean minor allele count in the gene, and within-gene linkage disequilibrium. To account for linkage disequilibrium between genes in close proximity, MAGMA models the residuals as a multivariate normal, with correlations set to the gene-gene correlations estimated as part of the gene-level analyses. In MAGMA, gene-gene correlations are estimated for pairs of genes within 5 Mb of each other and are otherwise set to zero. The gene-set statistic tests the null hypothesis that the mean association of case-control status with the genes in the set is greater than that of genes that are not in the set (i.e., H0, βS = 0 versus H1, βS > 0). Gene-set analyses were conducted using gene-level Z scores from each meta-analysis. A gene-set p value < 0.05 was considered significant.

Results

The distribution of CHD phenotypes in cases from the CTD and LVOTD case-parent trios are provided in Table 1. Sample sizes, total numbers of analyzed genes, and number of analyzed MEGs for the four analytic datasets are provided in Table 2. Pearson correlation coefficients for gene-level meta-p values, across analyses, are summarized in Table 3. With a single exception, the gene-level meta-p values from the different analyses were significantly correlated (p < 0.003). The correlations ranged from quite modest to very strong. The strongest correlations were observed between meta-p values from analyses of common variants only and analyses of common and rare variants combined (analyses of unfiltered variants, range 0.95–0.96; analyses of CADD10 filtered variants, range 0.80–0.81).

Table 1

CHD phenotypes in the four study datasets

Phenotype	CHOP, n (%)	PCGC, n (%)
Conotruncal heart defects
Tetralogy of Fallot	267 (39.9)	106 (29.9)
D-transposition of the great arteries	128 (19.1)	70 (19.7)
Ventricular septal defect	136 (20.3)	45 (12.7)
Double outlet right ventricle	76 (11.3)	57 (16.1)
Isolated aortic arch anomalies	30 (4.5)	7 (2.0)
Truncus arteriosus	20 (3.0)	13 (3.7)
Interrupted aortic arch	13 (1.9)	9 (2.5)
Other	0 (0.0)	48 (13.5)
Left ventricular outflow tract defects
Hypoplastic left heart syndrome	157 (49.5)	85 (44.3)
Coarctation of the aorta	97 (30.6)	67 (34.9)
Aortic stenosis	63 (19.9)	28 (14.6)
Other	0 (0.0)	12 (6.3)

CHOP, Children’s Hospital of Philadelphia; PCGC, Pediatric Cardiac Genomics Consortium.

Table 2

Summary of the datasets used in the gene-level analyses of the association between maternal genes and CHDs in offspring

Summary variables	CHOP-CTD	PCGC-CTD	CHOP-LVOTD	PCGC-LVOTD
Mothers (n)	650	309	301	190
Fathers (n)	543	312	253	190
Analyses with no CADD filter
# of genes analyzed (C/R/C + R)a	18,441/18,963/18,438	18,657/18,986/18,654	18,426/18,893/18,419	18,656/18,953/18,651
# of MEGs analyzed (C/R/C + R)a	79/79/79	79/79/79	79/79/79	79/79/79
Analyses with CADD ≥ 10 filter
# of genes analyzed (C/R/C + R)a	14,170/18,481/14,133	14,404/18,521/14,366	14,102/18,195/14,003	14,378/18,379/14,307
# of MEGs analyzed (C/R/C + R)a	57/79/57	59/78/59	58/78/58	58/78/58

CADD, combined annotation dependent depletion; CTD, conotruncal heart defect; LVOTD, left ventricular outflow tract defect; MEGs, maternal effect genes.

C, common variants (MAF ≥ 0.05); R, rare variants only (MAF < 0.05); C + R, common and rare variants.

Table 3

Gene-level meta-p-value correlations across analyses

Phenotype	Analysis	Correlations
		C Versus C + R	R Versus C + R	C Versus R	Unfiltered versus CADD10
CTD	unfiltered	0.96	0.20	0.05	–
	filtered	0.81	0.56	0.03	–
	common only	–	–	–	0.60
	rare only	–	–	–	0.20
	common + rare	–	–	–	0.44
LVOTD	unfiltered	0.95	0.21	0.07	–
	filtered	0.81	0.54	0.02	–
	common only	–	–	–	0.61
	rare only	–	–	–	0.22
	common + rare	–	–	–	0.45
CTD + LVOTD	unfiltered	0.95	0.20	0.05	–
	filtered	0.80	0.55	0.01a	–
	common only	–	–	–	0.61
	rare only	–	–	–	0.22
	common + rare	–	–	–	0.44

p = 0.16; all other p values < 0.0003.

Genome-wide analyses

The genome-wide, gene-level meta-analyses identified three candidate maternal genes for CHDs (i.e., genes with meta-p < 10−5 and less than the p value for each contributing dataset). For each of these genes, the statistical criteria for candidacy were met in the analyses of common variants only and in the analyses that included both common and rare variants. The gene-level meta-p values for both analyses are provided in Table 4. For simplicity, and because the gene-level meta-p values from these two analyses are highly correlated, we present only the lower of the two meta-p values in the text below.

Table 4

Candidate maternal genes for CHDs identified in genome-wide, gene-level analyses

Analysis	Genea	Filter	Variants	Bonferroni thresholdb	Meta-p	p values for contributing datasets
						CTD		LVOTD
						CHOP	PCGC	CHOP	PCGC
LVOTD	GABRB3	CADD10	common + rare	<2.7 × 10−6	4.0 × 10−6	–	–	4.0 × 10−3	1.1 × 10−4
			common	<3.4 × 10−6	6.5 × 10−6	–	–	4.0 × 10−3	4.0 × 10−4
LVOTD	SORCS3	unfiltered	common + rare	<2.6 × 10−6	4.8 × 10−7	–	–	6.0 × 10−4	2.4 × 10−3
			common	<2.7 × 10−6	5.3 × 10−7	–	–	8.5 × 10−4	2.0 × 10−3
CTD + LVOTD	SPRED3	unfiltered	common + rare	<2.6 × 10−6	2.5 × 10−6	3.6 × 10−3	2.6 × 10−3	1.6 × 10−1	2.8 × 10−3
			common	<2.7 × 10−6	3.4 × 10−6	3.0 × 10−3	3.0 × 10−3	1.7 × 10−1	2.6 × 10−3

Genes with meta-p values <10−5 and lower than the p values in each dataset included in the meta-analysis.

0.05/number of analyzed genes.

Two candidate maternal genes were identified in analyses of LVOTDs only: gamma-aminobutyric acid type a receptor subunit beta 3 (GABRB3), which encodes a member of the ligand-gated ionic channel family (CADD10 filtered variants, meta-p = 4.0 × 10−6); and sortilin-related VPS10 domain containing 3 (SORCS3), which encodes a type-I receptor transmembrane protein that is a member of the vacuolar protein sorting 10 receptor family involved in protein trafficking and cellular signaling in neuronal and non-neuronal cell (unfiltered variants, meta-p = 4.8 × 10−7). The third candidate was identified in the analyses that included both CTDs and LVOTDs: sprout-related EVH1 domain containing three (SPRED3), which encodes a member of a family of proteins that negatively regulates mitogen-activated protein signaling, particularly during organogenesis (unfiltered variants, meta-p = 2.5 × 10−6).

MEG analyses

Our analyses provided no evidence of an association between CHDs and the MEG gene set as a whole (gene-set p values: 0.13–0.95, Table S3). However, the MEG-specific, gene-level meta-analyses identified 14 candidate CHD-related MEGs (i.e., MEGs with meta-p < 0.05 and less than the p value for each contributing dataset). Details of the specific data and statistical evidence underlying the associations with these 14 genes are provided in Table 5. The 14 candidate CHD-related MEGs include three (CDC20, KHDC3L, TRIP13) known human MEGs and one (DNMT3A) MEG that has been associated with birth defects in a mouse model (Table 6). Only one of the 14 candidate CHD-related MEGs (KDM4A) was implicated as a candidate in analyses of both CTDs only and LVOTDs only. However, an additional five MEGs were implicated in either the analyses of CTDs only (KHDC3L, RNF2, and TRIP13) or the analyses of LVOTDs only (KMTD2 and PUM1), as well as in the combined CTD + LVOTD analyses. In general, the p values for these five MEGs were lower in the combined compared with the defect-specific analysis.

Table 5

Candidate CHD-related MEGs

Analysis	MEGa	Filter	Variants	Meta-p	p values for contributing datasets
					CTD		LVOTD
					CHOP	PCGC	CHOP	PCGC
CTD	EZH2	CADD10	common	3.40 × 10 −2	1.22 × 10−1	6.50 × 10 −2
	KDM4A	CADD10	common	3.72 × 10 −2	1.70 × 10−1	4.70 × 10−2
	KDM4A	CADD10	common + rare	3.47 × 10−2	1.59 × 10−1	4.10 × 10−2
	KHDC3L	unfiltered	rare	3.38 × 10−2	6.40 × 10−2	1.47 × 10−1
	RNF2	CADD10	common + rare	1.03 × 10−2	3.40 × 10−2	9.70 × 10−2
	TRIP13	CADD10	rare	1.28 × 10−2	2.00 × 10−2	1.29 × 10−1
LVOTD	CDC20	unfiltered	common	3.34 × 10−2			8.30 × 10−2	1.15 × 10−1
	CDC20	unfiltered	common + rare	3.58 × 10−2			1.27 × 10−1	7.20 × 10−2
	CDX2	unfiltered	common	1.02 × 10−2			2.00 × 10−2	1.69 × 10−1
	CDX2	unfiltered	common + rare	1.65 × 10−2			2.60 × 10−2	2.24 × 10−1
	DNMT3A	CADD10	rare	1.05 × 10−2			5.60 × 10−2	4.30 × 10−2
	IGF2BP2	CADD10	rare	4.67 × 10−2			1.06 × 10−1	1.16 × 10−1
	KDM4A	unfiltered	common	1.69 × 10−2			6.20 × 10−2	6.30 × 10−2
	KDM4A	unfiltered	common + rare	1.29 × 10−2			6.50 × 10−2	4.50 × 10−2
	KDM4A	CADD10	rare	6.34 × 10−3			2.00 × 10−2	8.60 × 10−2
	KMT2D	unfiltered	common	3.93 × 10−2			1.04 × 10−1	1.02 × 10−1
	KMT2D	unfiltered	common + rare	4.52 × 10−2			1.25 × 10−1	9.70 × 10−2
	PUM1	CADD10	rare	4.19 × 10−3			3.60 × 10−2	3.40 × 10−2
	SMARCA4	CADD10	common + rare	1.89 × 10−2			5.90 × 10−2	9.00 × 10−2
CTD + LVOTD	BCAS2	unfiltered	common	2.04 × 10−2	2.60 × 10−2	4.64 × 10−1	1.03 × 10−1	4.78 × 10−1
	CTCF	CADD10	rare	3.20 × 10−2	6.40 × 10−2	2.20 × 10−1	4.17 × 10−1	1.35 × 10−1
	KDM4A	unfiltered	common	2.53 × 10−2	5.19 × 10−1	6.60 × 10−2	6.20 × 10−2	6.30 × 10−2
	KDM4A	unfiltered	common + rare	2.21 × 10−2	5.19 × 10−1	6.50 × 10−2	6.50 × 10−2	4.50 × 10−2
	KDM4A	CADD10	common + rare	1.27 × 10−2	1.59 × 10−1	4.10 × 10−2	1.37 × 10−1	2.28 × 10−1
	KDM4A	CADD10	rare	1.32 × 10−2	2.24 × 10−1	2.74 × 10−1	2.00 × 10−2	8.60 × 10−2
	KHDC3L	unfiltered	rare	4.07 × 10−2	6.40 × 10−2	1.47 × 10−1	3.12 × 10−1	4.91 × 10−1
	KMT2D	unfiltered	common	3.03 × 10−2	1.80 × 10−1	2.80 × 10−1	1.04 × 10−1	1.02 × 10−1
	KMT2D	unfiltered	common + rare	3.67 × 10−2	1.93 × 10−1	2.94 × 10−1	1.25 × 10−1	9.70 × 10−2
	PUM1	CADD10	rare	3.58 × 10−3	2.58 × 10−1	5.60 × 10−2	3.60 × 10−2	3.40 × 10−2
	RNF2	CADD10	common	9.60 × 10−3	1.70 × 10−2	7.00 × 10−2	1.10 × 10−1	8.12 × 10−1
	TRIP13	CADD10	rare	2.52 × 10−3	2.00 × 10−2	1.29 × 10−1	3.20 × 10−2	3.72 × 10−1

MEGs with meta-p values <0.05 and lower than the p values in each dataset included in the meta-analysis.

Table 6

Description of candidate CHD-related MEGs

Gene symbol (mouse/human)	Gene name	Impact/function of maternal gene in mouse models	Human phenotype
Bcas2/BCAS2	BCAS2 pre-mRNA processing factor	developmental arrest at two- to four-cell stage; compromised DNA damage response in early embryos; accumulation of DNA damage and micronuclei28
Cdc20/CDC20	cell division cycle 20	needed for first mitotic division29	EEA, FF, OMA30
Cdx2/CDX2	caudal type homeobox 2	55% of embryos from maternal knockouts arrest at morula stage; 69% of embryos from maternal-zygotic knockouts arrest at morula stage; associated with cell death and specification of the trophectoderm31
Ctcf/CTCF	CCCTC-binding factor	increased zygotic lethality and increased methylation of the H19 differentially methylated domain32; increased meiotic and mitotic errors, differential zygotic gene expression, and apoptosis33
Dnmt3a/DNMT3A	DNA methyltransferase 3 alpha;	hypomethylation and aberrant expression of maternally imprinted genes in the embryo34; embryonic death before E11.5 with pericardial edema, lack of brachial arches, and open neural tube defects.34,35
Ezh2/EZH2	enhancer of zeste 2 polycomb repressive complex 2 subunit	severe growth retardation in neonates14; required for establishment of H3K27me3 in zygotes36
Igf2bp2/IGF2BP2	insulin-like growth factor 2 mRNA binding protein 2	subfertilility; 94% of embryos die prior to blastocyst stage; downregulation of transcription during zygotic genome activation37
Kdm4a/KDM4A	lysine demethylase 4a	majority of embryos fail to develop past the two- to four-cell stage; mediates H3K9me3 demethylation at broad domains of H3K4me3 in oocytes; associated with decreased expression of genes linked to zygotic genome activation38,39
Khdc3/KHDC3L	KH domain containing 3 like	delays pre-implantation development; high incidence of aneuploidy; inactivation of spindle assembly checkpoint40; member of the subcortical maternal complex in mouse and human41	EEA, HM, RPL42, 43, 44; DNA methylation loss and imprinting defects45
Kmt2d/KMT2D	lysine methyltransferase 2D	majority of embryos arrest at one- to two-cell stage; controls promoter-specific chromatin modification during oogenesis and early development46
Pum1/PUM1	Pumilio RNA binding family member 1	abnormal pre-implantation development; regulation of maternal mRNA47
Rnf2/RNF2	ring finger protein 2	double-maternal knockout of Rnf2 and Ring2 (components of polycomb repressive complex 1) results in embryonic arrest at 2-cell stage and impaired zygotic genome activation48
SMARCA4	SWI/SNF-related, matrix-associated, actin dependent regulator of chromatin, subfamily A, member 4	embryonic arrest at two-cell stage; impaired zygotic genome activation49
Trip13/TRIP13	thyroid hormone receptor interactor 13	oocyte loss; defect in double-strand break repair	OMA, ZCF50

EEA, early embryonic arrest; FF, fertilization failure; HM, hydatidiform mole; IF, implantation failure; OMA, oocyte maturation arrest.

Description of genes in Table 5.

Discussion

In mammals, MEGs have been associated with both early embryonic (e.g., zygotic cleavage failure) and later offspring (e.g., multi-locus imprinting disorders) phenotypes. In addition, in mouse models, a subset of MEGs has been associated with structural birth defects, including CHDs. Our prior work also suggested a potential association between MEGs and CHDs in humans. However, to our knowledge, this is the first study to address the a priori hypothesis that, in humans, maternal variation in MEGs is associated with the risk of CHDs in offspring.

Both our genome-wide and MEG-specific analyses provide evidence that maternal genetic variation may be associated with CHDs in offspring. In the genome-wide assessment, three candidate maternal genes for CHDs were identified (GABRB3, SORCS3, and SPRED3). These genes were not on our pre-specified list of mammalian MEGs. However, according to the database of the transcriptome in mouse early embryos, they are expressed in the oocyte and early (one-, two-, and four-cell) embryo, which is consistent with the expectation for MEGs.

In the MEG-specific analyses, 14 candidate CHD-related MEGs were identified. These candidates include three (27%) of the 11 known human MEGs: CDC20, KHDC3L, and TRIP13. Evidence that CDC20 and TRIP13 are human MEGs comes from studies in women with infertility due to oocyte maturation abnormalities, zygotic cleavage failure, or early embryonic arrest.,, For KHDC3L, which encodes a member of the subcortical maternal complex (a multi-protein complex that is uniquely expressed in the oocyte and early embryo), biallelic mutations have been identified in approximately 5% of women with a history of recurrent hydatidiform moles. Mutations in KHDC3L have also been reported in women with a history of recurrent pregnancy loss.

The candidate CHD-related MEGs also included one (12%) of the eight MEGs that have been associated with structural birth defects: DNMT3A. In mouse oocytes, disruption of Dnmt3a results in hypomethylation of maternally imprinted genes. This hypomethylation is inherited by the embryo, and results in aberrant expression of these maternally imprinted genes. Embryos from females homozygous for an oocyte conditional knockout of Dnmt3a die before embryonic day 11.5 and present with open neural tube defects, lack of brachial arches, and pericardial edema.,

The lysine demethylase, KDM4A, was also identified as a candidate CHD-related MEG. This was the only MEG for which there was evidence of association in both the CTD-only and the LVOTD-only analyses. Although this gene has not been associated with birth defects in mouse models, such an association has been reported for another lysine demethylase, Kdm1b. Mice homozygous for a Kdm1b knockout develop normally, but offspring of female knockouts fail to develop past mid-gestation and present with placental defects, neural tube defects, pericardial edema, and growth impairment. In addition, the liveborn offspring of female mice homozygous for a hypomorphic mutation in yet another lysine demethylase, Kdm1a, have a high neonatal death rate (26% versus 5% in offspring of wild-type females) that is consistent with the animals having subtle developmental defects. Further, several studies indicate that absence of the maternal product for lysine demethylase genes impairs embryonic genome activation.,, Consequently, maternal variation in these genes could affect heart development by altering the expression of relevant embryonic genes.

Most studies of MEGs in humans have focused on rare, predicted deleterious variants in protein-coding genes. However, common variants in protein-coding genes have also been associated with reproductive outcomes including recurrent implantation failure and recurrent pregnancy loss., Our analyses identified associations between both common and rare variants in MEGs and CHDs. For example, evidence for KDM4A as a candidate in the combined analysis of CTDs and LVOTDs was obtained from both the analysis of common variants only (unfiltered, meta-p = 0.025) and the analysis of rare variants only (CADD10 filtered, meta-p = 0.013). Our analyses also identified associations with the combined phenotypes of CTDs and LVOTDs. Although some of the candidate MEGs were associated with only one phenotype, eight of the 14 candidates were identified in the combined analyses, suggesting that at least some MEGs may be associated with CHDs in general rather than with specific CHD phenotypes.

As for any study, the results of our analyses should be interpreted within the context of the existing body of knowledge. However, the role of the maternal genotype in determining reproductive and offspring outcomes has not been extensively studied. Given this gap, our ability to draw conclusions from our analyses is limited. In particular, although associations between the maternal genotype and offspring phenotype would be expected if the gene acts as an MEG, there are alternative explanations. For example, in addition to the possibility of false-positive findings, such associations could result from an effect of the maternal genotype on the uterine environment, which in turn affects embryonic development, or from a parent-of-origin effect whereby the consequences of the inherited allele depend on the sex of the transmitting parent.

For genes that are known MEGs, it is reasonable to assume that an observed association (if true) reflects the gene’s action as an MEG. Such an assumption is, however, less justified for associations that do not involve a known MEG (e.g., the associations with GABRB3, SORC3, and SPRED3 identified in our genome-wide analyses). We found no evidence in the literature that helped to differentiate between the potential explanations for the observed associations of CHDs with SORC3 or SPRED3. Further, while there is some evidence that maternal GABRB3 may be associated with autism spectrum disorder in offspring, the mechanism underlying this association is also unclear. Consequently, while our studies are suggestive of associations between the maternal genotypes for these three genes and CHDs, further studies are required to confirm these associations and, if confirmed, to identify the mechanisms that underlie the associations.

The results of our analyses must also be interpreted in light of the relatively large number of comparisons that were performed. Because consideration of MEGs as potential risk factors for birth defects is a new area of inquiry, we opted for a comprehensive assessment of both phenotype (i.e., consideration of CTDs and LVOTDs, separately and combined) and genotype (i.e., consideration of rare and common variants, separately and combined, and with and without CADD10 filtering). These comprehensive analyses allowed us to draw preliminary conclusions regarding the nature of MEG-CHD associations that would not have been possible under a more focused analysis plan. For example, the results of our analyses suggest that both common and rare variation in MEGs may be associated with CHDs and that some MEGs may be associated with a broad array of CHD phenotypes. To address concerns arising from multiple testing, we required MEG candidates to have meta-p values lower than our pre-specified criteria for suggestive evidence and for the meta-p value to be less than the p values from each of the contributing datasets. The latter criterion ensures that the meta-p values are not being driven by the results from a single dataset. Nonetheless, our criteria for declaring candidates were lenient. Specifically, we used a Bonferroni correction to account for the number of genes within an analysis, but did not account for the number of different analyses that were performed. Consequently, the results from these analyses should be considered preliminary.

Conclusions

Our analyses add to the growing literature on the role of the maternal genotype in mammalian development. In particular, we provide preliminary evidence that MEGs are associated with CHDs in humans. Further study of the association of the maternal genotype with CHDs and other structural birth defects may provide an explanation for the large proportion of birth defects for which etiology is unknown. Importantly, understanding the role of the maternal genotype could significantly affect genetic and reproductive counseling. For example, for women carrying damaging MEG genotypes, oocyte donation may be an appropriate risk-reducing reproductive option. Given the personal and public health impact of CHDs and other structural birth defects and the large proportion of cases that remain unexplained, research aimed at enhancing our understanding of the role of the maternal genotype in the etiology of these conditions is strongly warranted.

Data and code availability

This study did not generate datasets or code.

The genotype data used in these studies are available at:

PCGC CTDs and LVOTDs:

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs001194.v2.p2.

CHOP CTDs:

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000881.v1.p1.

CHOP LVOTDs:

https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000781.v1.p1.