Common Genetic Variants Contribute to Risk of Transposition of the Great Arteries.

RATIONALE: Dextro-transposition of the great arteries (D-TGA) is a severe congenital heart defect which affects approximately 1 in 4,000 live births. While there are several reports of D-TGA patients with rare variants in individual genes, the majority of D-TGA cases remain genetically elusive. Familial recurrence patterns and the observation that most cases with D-TGA are sporadic suggest a polygenic inheritance for the disorder, yet this remains unexplored.
OBJECTIVE: We sought to study the role of common single nucleotide polymorphisms (SNPs) in risk for D-TGA. METHODS AND
RESULTS: We conducted a genome-wide association study in an international set of 1,237 patients with D-TGA and identified a genome-wide significant susceptibility locus on chromosome 3p14.3, which was subsequently replicated in an independent case-control set (rs56219800, meta-analysis P=8.6x10-10, OR=0.69 per C allele). SNP-based heritability analysis showed that 25% of variance in susceptibility to D-TGA may be explained by common variants. A genome-wide polygenic risk score derived from the discovery set was significantly associated to D-TGA in the replication set (P=4x10-5). The genome-wide significant locus (3p14.3) co-localizes with a putative regulatory element that interacts with the promoter of WNT5A, which encodes the Wnt Family Member 5A protein known for its role in cardiac development in mice. We show that this element drives reporter gene activity in the developing heart of mice and zebrafish and is bound by the developmental transcription factor TBX20. We further demonstrate that TBX20 attenuates Wnt5a expression levels in the developing mouse heart.
CONCLUSIONS: This work provides support for a polygenic architecture in D-TGA and identifies a susceptibility locus on chromosome 3p14.3 near WNT5A. Genomic and functional data support a causal role of WNT5A at the locus.

In This Issue, see p

Meet the First Author, see p

Editorial, see p

Complete transposition of the great arteries (TGA), also referred to as dextro-TGA (D-TGA), is a developmental cardiac outflow tract (OFT) defect.[1] With an incidence of 20 to 30 per 100 000 live births, it is among the most common forms of severe, cyanotic congenital heart disease (CHD).[1,2] The anatomy of D-TGA is characterized by ventriculoarterial discordance, and about one-third of cases present additional lesions such as ventricular septal defect (VSD) and pulmonary outflow tract obstruction.[3] If untreated, D-TGA leads to severe hypoxemia and neonatal mortality. Notably, advances in surgical treatments have improved patient survival and have led to a substantial increase in adults with surgically corrected D-TGA, albeit with significant long-term morbidity in some patients.[4] Given that most D-TGA patients now reach reproductive age, there is a growing need for a better understanding of the genetic architecture of D-TGA to allow for effective reproductive counseling.

While there are several reports of D-TGA patients with rare variants in individual genes,[5-10] the majority of D-TGA cases remain genetically elusive. Familial recurrence patterns and the observation that most cases with D-TGA are sporadic suggest a polygenic inheritance,[11-16] with genetic variants of different frequency and effect size contributing to risk. In recent years, genome-wide association studies (GWAS) have uncovered common genetic variants contributing to susceptibility for specific CHD lesions (eg, atrial septal defect,[17] tetralogy of Fallot,[18] and left-sided lesions[19,20]). Although some of these studies included patients with diverse CHD lesions, including D-TGA,[17,21-23] thus far no GWAS has been realized selectively in patients with D-TGA. We here conducted a case-control GWAS on 1237 cases with D-TGA recruited in Europe, North America, and Australia. We identified a susceptibility locus at 3p14.3 at genome-wide statistical significance and demonstrated that 25% of variance in D-TGA susceptibility was attributable to common genetic variation, supporting an important role for such variants in susceptibility to D-TGA. Functional analysis of the 3p14.3 susceptibility locus provided compelling evidence for a causal role of WNT5A, encoding the Wnt Family Member 5A protein, known for its critical role in cardiac development.[24]

Methods

Data Availability

A detailed, expanded Methods description can be found in the Supplemental Material. For details on all essential research materials used in this study, please see the Major Resources Table in the Supplemental Material. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Study Samples

We included unrelated individuals with D-TGA as the main cardiac lesion. For prespecified subgroup analyses, we subcategorized the cases into D-TGA in the absence of a VSD and pulmonary stenosis, referred to as “simple D-TGA,” and D-TGA with associated VSD and/or pulmonary stenosis, referred to as “complex D-TGA.” Cases for the genome-wide association analysis in the discovery set were included from centers in the Netherlands, France, the United Kingdom, Canada, Germany, Australia, and Belgium. For the purpose of replication, we included independent cases from the Netherlands and the United States. Previously genotyped controls were included from various studies on the basis of country of origin and genotyping platform to match the cases. An overview of genotyping arrays used in the discovery and replication sets are provided in Table S1. All cases and controls gave written informed consent. The study was approved by the relevant institutional ethical review boards.

Association Analysis

After quality control (QC) and genotype imputation, the association of alternate allele dosage with D-TGA in the discovery set was tested using logistic regression assuming an additive genetic model, adjusting for sex and the first 5 principal components (PCs) using PLINK 1.9.[25] Single nucleotide polymorphisms (SNPs) associated with a P<5×10−8 were considered genome-wide significant in the discovery cohort and were taken forward for replication, with the threshold for statistical significance set to P<0.05/number of tests. An exploratory genome-wide meta-analysis was performed, where SNPs with P<5×10−8 were considered genome-wide significant and those with a P<10−5 were considered suggestive.

Estimation of SNP Heritability and Polygenic Score

Single nucleotide polymorphism (SNP)-based heritability of D-TGA was estimated in the discovery set after additional stringent post-imputation QC as suggested[26] and using Phenotype-Correlation-Genotype-Correlation regression,[27] assuming a prevalence of 0.03%.[2] Furthermore, we derived a genome-wide polygenic risk score (PRS) in the discovery set using LDpred[28] and assessed the association of this PRS with D-TGA in the replication set using logistic regression, correcting for sex.

LacZ Enhancer Screen in Mice

A predicted regulatory element (RE1) overlapping the genome-wide association signal (in mice: chr14:29,229,674-29,231,629; mm9) was cloned into Hsp68-LacZ,[29] and mouse in vivo transgenic reporter assays were performed by pronuclear injection as reported previously.[30] Embryos were harvested at E10.5 and subjected to X-gal staining as reported previously.[30] Embryos were genotyped using LacZ primers (forward, 5′-CAGCCTGAATGGCGAATGGCGCTT; reverse, 5′-CCCGTTGCACCACAGATGAAACGC). All embryos with X-gal staining were included for subsequent analysis, sex was not determined in the embryos. Only patterns that were observed in at least 3 different embryos resulting from independent transgenic integration events of the same construct were considered reproducible.[29,31]

Enhancer Screen in Zebrafish

An amplified fragment of human RE1 was cloned into a custom DNA vector containing an E1b-mCherry-sv40pA reporter cassette and an attB recombination site for phiC31-mediated targeted transgenesis. The constructs were then microinjected in presence of phiC31 mRNA in 1-cell zebrafish embryos carrying a single genomic attP landing site. RE1-driven mCherry expression was assessed in 2dpf zebrafish embryos. The minimal number of injected analyzed embryos was 20, based on previous work in our laboratory.[32,33]

Deletion of the GWAS Locus in Mice

Mice lacking the region including RE1 were generated using CRISPR/Cas9 genome editing. Whole-mount nonradioactive in situ hybridization analysis was performed to assess differences in expression of Wnt5a between wild-type and knock-out mice. Transcript levels of candidate target genes were determined by quantitative polymerase chain reaction. Sample size was not predetermined statistically and is based on previous work with a minimum of 5 animals per condition for each experiment.[34,35]

TBX20 Chromatin Immunoprecipitation Sequencing

TBX20 chromatin immunoprecipitation sequencing data from E11.5 mouse hearts was generated previously.[30] These data are accessible from the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress, accession number E-MTAB-3967[30]). Histone modification data was generated previously by the ENCODE consortium.[36]

Cardiac Tbx20 Conditional Knock-Out Mice

RNA sequencing after fluorescence-activated cell sorting, on cardiomyocytes from conditional cardiac Tbx20 knock-out mice (Tnnt2-rtTA;TetO-Cre;Tbx20) and control cardiomyocytes was generated previously.[37] These data are accessible from the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress, accession number E-MTAB-5596). Analysis was performed utilizing the limma software package (3.26.8).[38] Raw read counts were first quantile normalized and transformed to the log-scale using the voom function.[39] Next, differential expression analysis was performed, and P were multiplicity corrected using Benjamini-Hochberg correction.

Transient Transfection and Luciferase Assay In Vitro

To generate reporter constructs for the 3 common haplotypes, the region including the putative regulatory element (RE1) was amplified from genomic DNA of homozygous carriers and cloned into the pGL4-SCP1-luc vector. HL-1 cells[40] were transfected with the reporter plasmid DNA, and luciferase activity measurements were performed. Statistics were performed using a Kruskal-Wallis nonparametric test.

Results

Description of Study Population

A total of 1555 samples from patients with D-TGA were included, originating from 8 countries (Table S1): (1) a discovery set that consisted of 1307 patients with D-TGA that were recruited at multiple centers in the Netherlands (until June 2018), Germany, France, Belgium, the United Kingdom, Canada, and Australia; and (2) a replication set that consisted of 248 patients with D-TGA recruited in the Netherlands (after June 2018) and the United States. Approximately two-thirds (n=1056) had D-TGA in the absence of a VSD and pulmonary stenosis, referred to as “simple D-TGA”, and the remainder (n=499) had an associated VSD and/or pulmonary stenosis, referred to as “complex D-TGA.” Controls were obtained from previously genotyped datasets that best matched the cases in terms of ancestry and genotyping platform (total n=13 023). Although no systematic cardiac evaluation was performed in this control population, it is expected that a negligible portion of controls have D-TGA given its low prevalence and severe presentation.

Genome-Wide Association Analysis

Genotypic data of the discovery set was either generated in this study or was obtained previously[17] (see Supplemental Methods and Table S1). QC exclusions were applied at both the SNP and at the sample levels. We restricted our analyses to cases and controls of European descent (as this comprised the majority of cases) and assigned each case to 3 ancestry-matched controls (see Supplemental Methods). Details on QC and principal component analyses are presented in Tables S2 and S3, and Figures S1 through S3. After QC and imputation using the Haplotype Reference Consortium panel,[41] a total of 5 013 439 high-quality common SNPs (minor allele frequency ≥5%) in 1094 cases and 3282 controls remained for analysis in the discovery set. The GWAS quantile-quantile plot is shown in Figure S4 and summary results after genomic correction (λ=1.03) in the discovery set are shown as a Manhattan plot in Figure 1A. A locus on chromosome 3p14.3 reached genome-wide statistical significance in the logistic regression (lead SNP rs56219800; P=1.3×10−8; imputation quality score R2=0.97; Figure 1C). The alternate allele C conferred an odds ratio (OR) of 0.69 (95% CI, 0.61–0.79; minor allele frequency 0.19 in cases and 0.25 in controls).

Figure 1.

Manhattan plots of the genome-wide association results. P of all tested single nucleotide polymorphisms (SNPs) on a –log10 scale from the logistic regression analysis in the discovery set including 1094 cases and 3282 controls (A) and the meta-analysis of the discovery and replication set including 1237 cases and 3711 controls (B). Dashed lines represent genome-wide significance P thresholds of 5×10−8. C, Regional association plot of the 3p14.3 locus and surrounding genes (±500 kb). The y-axis represent the −log10(P) values from the meta-analysis. Round points represent SNPs in this region, colored by degree of linkage disequilibrium with rs56219800, the top associated SNP (purple point). The blue highlighted region represents the size of the locus, as defined by all SNPs with an r2≥0.5 with the lead SNP and that showed an association P of <10−3.

The association of rs56219800 with D-TGA was replicated in an independent replication set that consisted of 143 cases and 429 controls after QC (P=0.02; OR=0.66 [95% CI, 0.46–0.93]). In a meta-analysis, combining the discovery and replication case-control sets, rs56219800 remained the only signal reaching genome-wide statistical significance (P=8.6×10−10; OR=0.69 [95% CI, 0.61–0.77]; Figure 1B). Conditional analysis, correcting for the lead SNP, did not reveal any additional independent significant (P<1×10−5) signal at the 3p14.3 locus. Genome-wide case-control meta-analysis separately in simple D-TGA (N=841; Figure S5A) and complex D-TGA (N=396; Figure S5B) showed similar association effect sizes of rs56219800 (OR=0.69 [95% CI, 0.60–0.80], P=8.8×10−7 in simple D-TGA; and OR=0.68 [95% CI, 0.55–0.83], P=1.6×10−4 in complex D-TGA) and did not detect any associations at genome-wide statistical significance.

A total of 21 loci displayed a suggestive association (P<1×10−5) with D-TGA in the meta-analysis (Table S4). Notably, none of the 7 SNPs previously associated with other CHD lesions[17-20,22] were associated with D-TGA (Table S5).

SNP Heritability Estimates and Polygenic Risk Score

To assess the overall contribution of common variants to the risk of D-TGA, we estimated the SNP heritability of D-TGA in our case-control discovery set. We used Phenotype-Correlation-Genotype-Correlation regression, designed to account for the bimodal liability distribution in an ascertained cohort.[27] Assuming a prevalence of 0.03% for D-TGA,[2] we estimated SNP heritability on a liability scale of 0.25 ([95% CI, 0.16–0.35], P=1.12×10−14). Similar results were obtained when varying the prevalence from 0.01% to 0.06%, as well as when restricting the analysis to the Dutch subset genotyped on the same platform (n=455 cases and 1348 controls) (Figure S6).

To further explore the polygenicity of D-TGA, we derived a genome-wide PRS using the summary statistics of the case-control discovery set using the Gibbs sampler implemented in LDpred[28] and confirmed its association in the replication set of 141 cases and 418 controls (P=4×10−5, C-statistic=0.63, R2=0.06; for ρ=0.0003). Adding the first 20 genotypic PCs to the model still resulted in a significant association of PRS with D-TGA (P=1.25×10−5). The association was consistent across different levels of ρ (tuning parameter of the Gibbs sampler) and when using P thresholding and R2 pruning for PRS derivation (Table S6). Exclusion of chromosome 3 (including our genome-wide locus) from the PRS still resulted in a significant association with D-TGA in the replication set (P=5×10−5, C-statistic=0.62, R2=0.06; for ρ=0.0003), indicating that the genome-wide common variant polygenicity signal was not restricted to the 3p14.3 locus.

Annotation of Associating Loci

We assigned candidate genes to a locus if they were located in the same topologically associating domain (TAD) and/or when their expression levels were correlated with genetic variants within the association signal (cis-expression quantitative trait locus, expression quantitative trait locus [eQTL]) (see Supplemental Methods). The genome-wide locus at 3p14.3 is located within the last intron of ERC2 and lies ≈70 kb from the transcription start site of WNT5A (Figure 1C). These 2 genes are the only ones located within the same TAD as the association signal (Figure 2). In the human genotype tissue expression dataset (GTEx V8), multiple SNPs, including the lead SNP, were eQTLs for WNT5A in esophageal mucosa, with the D-TGA risk allele being associated with higher expression of WNT5A (P=2×10−5 for rs56219800). There was no significant cis-eQTL for the other gene in the TAD (ie, ERC2) or any other neighboring protein-coding gene. We used EMERGE, an in silico tool incorporating over 50 human cardiac epigenetic datasets to identify sequences within the associating locus with gene regulatory potential in the heart.[42] This demonstrated that the 3p14.3 association signal colocalized with an EMERGE-predicted regulatory element (located at chr3:55,605,229-55,607,027; GRCh37), hereafter referred to as RE1 (Figure 2). Virtual chromosome conformation capture-on-chip data (4C), as derived from Hi-C data,[43] obtained in H9 human embryonic stem cells during cardiomyocyte differentiation (see Supplemental Methods for more detail) showed chromatin interaction between the RE1 region and the promoter of WNT5A (Figure S7). WNT5A encodes Wnt family Member 5A, a transcription factor involved in organ development. Wnt5a null mice demonstrate perinatal lethality with multi-organ abnormalities[44] and exhibit outflow tract defects, including D-TGA.[45] Collectively, these data point to WNT5A as the causal gene at this locus and suggest that it might be regulated by RE1.

Figure 2.

Functional annotation of the 3p14.3 locus. From top to bottom: Hi-C data from the human lymphobastoid line GM12878[79]; TADs derived from iPSC-CMs[80]; RefSeq genes; results from the TGA genome-wide association studies meta-analysis on a –log10(P) scale; human EMERGE track[42]; Genomic Evolutionary Rate Profiling (GERP) conservation track.[81] The blue highlighted region represents the genome-wide associated 3p14.3 locus near WNT5A.

In the additional 21 loci that passed the suggestive significance threshold (P<10−5), a total of 93 genes were localized within the same TAD as the association signal and/or displayed a cis-eQTL effect with the association signal (Supplemental Methods and Table S7). At each locus at least one gene was implicated in either abnormal heart morphology in mouse and/or humans (Table S8) and/or was present in a list of CHD (candidate) genes[46-48] (Table S7). Thirteen of the 21 suggestive association signals colocalized with an EMERGE signal and could therefore contain regulatory elements regulating the expression of genes involved in the pathogenesis of D-TGA (Table S4).

Functional Characterization of the 3p14.3 Locus

To determine if the human EMERGE-predicted regulatory element at the 3p14.3 locus (ie, RE1) could drive gene expression during cardiac development in mice, we surveyed histone modification marks in the developing mouse heart, and found that RE1 displays hallmarks of a bona fide regulatory element, including H3K27ac associations (Figure 3A). In line with this, similar to the observation made in human above, the corresponding chromosomal region in the mouse showed a strong EMERGE signal indicative of a cardiac regulatory element (Figure S8A). This region is homologous to human RE1, indicative of functional conservation of this element. To test the functionality of the mouse RE1, we performed enhancer reporter assays in vivo and observed that RE1 drives robust cardiac expression in all transgenic mouse embryos (n=5, Figure 3B). Notably, expression was observed at E10.5 in outflow tract myocardium (5/5 lacZ positive embryos) and ventricles (5/5), regions that normally express Wnt5a at this stage.[49] The ability of RE1 to drive gene expression in the heart (among others) was confirmed in zebrafish using the human RE1 fragment (Figure 3C).

Figure 3.

Characterization of RE1 in mouse and zebrafish. A, Genome browser view from the mouse orthologous Wnt5a genomic region. ChIP-seq tracks of TBX20 and H3K27Ac are displayed below. B, Representative embryos showing that RE1-lacZ reporter expression is observed in the developing murine heart, including outflow tract myocardium and ventricles. In total, 5 embryos were analyzed. A, atrium; V, ventricle; OFT, outflow tract. C, Representative image of 2dpf zebrafish showing RE1 activity in various organs and tissues, including the heart. In total, 296 embryos were analyzed. White dashed lines represent boundaries of ventricle and atrium. v: ventricle; a: atrium. D, TBX20 ChIP-Seq demonstrating RE1 to be bound by TBX20. E, Expression levels (average normalized read counts and SD) of Wnt5a and flanking genes in E11.5 FACS sorted cardiomyocytes from control (n=4) and cTnT-Cre;Tbx20 mutant murine hearts (n=4), *P=0.04, using limma[38] adjusted for multiple testing using Benjamini-Hochberg correction.

We next assessed whether homozygous deletion in mouse of the region homologous to the associating locus at 3p14.3, including RE1, would result in a cardiac phenotype. Two independent lines of RE1−/− mice were born at a normal Mendelian ratio and showed no CHD, nor an effect on Wnt5a expression (Figure S8). To test whether RE1 might affect other, nearby genes, transcription levels of the other 3 genes in the same or adjacent TAD (ie, Lrtm1, Cacna2d3 and Erc2) were also assessed at the same embryonic stage, and all were found to be almost undetectable in the heart (Figure S8B). We hypothesize that the lack of phenotypic differences in these RE1−/− mice may be due to enhancer redundancy, as observed for multiple other developmental (cardiac) genes.[50,51]

A closer look at the mouse locus using existing epigenetic and chromatin immunoprecipitation sequencing datasets[30,42] obtained in mouse demonstrated that RE1 is bound by TBX20 in vivo in mouse hearts of E11.5 (Figure 3D), a time point at which outflow tract septation and alignment occurs. To gain more insight into a possible role of TBX20 in the regulation of Wnt5a, we compared Wnt5a expression in cardiomyocytes obtained from conditional, cardiac-specific, Tbx20 knock-out mice, and wildtype littermates (n=4 in both groups).[37] Notably, the expression of Wnt5a was significantly increased in the absence of TBX20 (P=0.04, Figure 3E). Moreover, other genes in the vicinity of RE1 were either not expressed above threshold levels, or their expression was not changed between Tbx20 mutant and control cardiomyocytes (Figure 3E). In summary, we show that (1) RE1, which overlaps the GWAS association signal at 3p14.3, robustly drives expression in the developing cardiac outflow tract in vivo, in a pattern that recapitulates Wnt5a expression, (2) RE1, which interacts with the promoter of WNT5A, is bound by TBX20, and (3) TBX20 attenuates Wnt5a expression levels in the developing mouse heart. Together with the fact Wnt5a transgenic mouse models[45,49,52] and Tbx20 null mice[53] exhibit outflow tract abnormalities, this indicates a potential causal role of WNT5A at the GWAS locus.

The above observations pointed to RE1 as a possible location for the causal variant of the 3p14.3 GWAS association signal. We cloned the 3 major haplotypes overlapping RE1 (between them accounting for almost 75% of alleles at this locus in the European sub-population of the 1000 Genomes Project; see Supplemental Methods), harboring 7 different SNPs with varying r2 with the lead SNP (Figure 4A and Figure S9) and tested their regulatory activity. In comparison to the empty vector, all 3 human haplotypes drove robust luciferase reporter gene activity in HL-1 cells, a mouse atrial cardiomyocyte tumor lineage cell line. One of them, haplotype 3, demonstrated a significantly increased luciferase reporter activity of 2.2-fold (P=2.6×10−5) and 2.4-fold (P=2.6×10−6) compared with haplotype 1 and haplotype 2 respectively (Figure 4B). While haplotypes 2 and 3 differed at multiple sites, haplotypes 1 and 3 only differed at one SNP, that is, rs79492708. Interestingly, this SNP was predicted to abrogate a bona fide TBX20 binding site with high confidence (Figure 4C). Despite the observed increase in reporter activity for rs79492708 and the likely role of TBX20 in modulating this activity, it is unlikely that this SNP alone drives the observed association in the GWAS as the r2 between the lead SNP and rs79492708 was only 0.07 (Figure 4A). We were unable to determine the independent effect of rs79492708 on D-TGA status in patients as it has a D prime of 1 with the lead SNP from GWAS (ie, with rs56219800). Additional in-depth analysis of currently available epigenetic data did not identify other putative regulatory elements overlapping the association signal (45 kb, defined by the lead SNP and SNPs with r2≥0.5) (Figure 2), precluding the targeted analysis of additional candidate SNPs.

Figure 4.

Regulatory activity of RE1 for the 3 major human haplotypes in vitro. A, Genome browser view from the human RE1 region (ie, chr3:55,605,229-55,607,027; GRCh37), depicting the location and characteristics of the single nucleotide polymorphisms (SNPs) within RE1 and the 3 major haplotypes. Provided r2 and D-prime are in relationship to the lead SNP rs56219800. B, Luciferase activity were measured in HL1 cells transfected with constructs for RE1 haplotype 1 (n=4), haplotype 2 (n=8), and haplotype 3 (n=4). Values were normalized to the activity of empty pGL4-SCP1-luc. *P=2.6×10−5, **P=2.6×10−6 (1-way ANOVA). C, TBX20 motif sequence as identified by HOMER motif analysis showing rs79492708 abolishes a TBX20-binding site.

While the exact causal SNP(s) remains to be elucidated, in aggregate the above studies in mice, zebrafish and HL-1 cells demonstrate functionality of a noncoding regulatory element overlapping the genome-wide signal at 3p14.3 and support an effect of the locus through modulation of expression of WNT5A, a gene strongly implicated in OFT development.

Discussion

We conducted the first case-control GWAS in patients with D-TGA. We identified a genome-wide significant locus on chromosome 3 that was subsequently replicated in an independent set of patients with D-TGA. Using SNP-based heritability analysis, we demonstrated that a quarter of variance in susceptibility to D-TGA may be attributable to common genetic variation, supporting a complex genetic inheritance. Through in vitro cellular and in vivo mouse and zebrafish studies, we provided multilevel evidence that the 3p14.3 locus contains a regulatory element that robustly drives expression in the developing cardiac outflow tract. Furthermore, our data support an effect of this regulatory element on the expression of WNT5A, a gene strongly implicated in OFT development, providing evidence for a causal role of this gene in susceptibility to TGA.

We studied the contribution of common genetic variants to susceptibility for D-TGA in individuals of European descent, leading to the identification of a risk locus on chromosome 3. Although the relatively small sample size of patients with this rare disease limited our statistical power to identify additional genetic variants at the conventional genome-wide statistical significance level, the SNP-heritability analysis we conducted demonstrated that 25% of variance in risk (on the liability scale) for D-TGA in our sample was attributable to additive effects of (the studied) common variants. Furthermore, a genome-wide polygenic risk score derived from the discovery set was associated with D-TGA in the replication set, also after removal of the 3p14.3 locus. Collectively, these results support a complex genetic inheritance and suggest that an accumulation of many common genetic variants together are likely to have a substantial effect on risk of D-TGA. Low frequency or rare variants (minor allele frequency <1%) likely also contribute to risk for D-TGA, although they are hard to identify with the current cohort sizes. Moreover, genetic factors predisposing to D-TGA likely comprise other variant classes besides single nucleotide variants (eg, copy number variants[54,55]) and other genetic effects besides inherited genetic effects (eg, maternal genetic factors[21] and de novo variants).

WNT5A emerged as the most likely causal gene at the 3p14.3 locus, based on chromatin topology at the locus, eQTL analysis, and its known function in outflow tract development. WNT5A encodes one of the 19 WNT (Wingless/integrase 1) family glycoproteins and activates several Wnt signaling pathways essential for normal development.[24] In humans, autosomal dominant loss-of-function mutations in WNT5A lead to Robinow syndrome,[56,57] characterized by skeletal and facial abnormalities. While this syndrome typically does not include CHD, one report described a fetus with a homozygous WNT5A loss-of-function mutation and persistent truncus arteriosus, a severe OFT defect.[58] In mice, it has been shown that Wnt5a is expressed in secondary heart field progenitor cells in the caudal splanchnic mesoderm, adjacent to the cardiac OFT.[49] Complete absence of Wnt5a in mice resulted in multiorgan abnormalities[44] and severe cardiac OFT defects, including D-TGA.[45,49] Additionally, Wnt5a expression throughout the entire splanchnic mesoderm led to OFT shortening and cardiac looping defects.[52] WNT5a is involved in many different (cardiac) developmental pathways. For instance, WNT5a is known to activate the planar cell polarity pathway[59,60] and modulate cell cohesion[52] in the secondary heart field creating a push and pulling force, essential for proper OFT elongation and subsequent looping and alignment over the inter-ventricular septum. Furthermore, it has also been suggested that WNT5a acts as a morphogen providing directional signals from the pharyngeal mesoderm to the adjacent cardiac neural crest cells,[45] and WNT5a can inhibit the canonical (β-catenin dependent) Wnt signaling pathway.[61,62] The ablation of β-catenin, the main effector of canonical Wnt signaling, from the developing mouse heart led to disturbed secondary heart field development and subsequently to abnormal cardiac looping and right ventricle formation.[63] Interestingly, Wnt5a is regulated by Tbx1,[64] which is essential for OFT development[65] and the major player in DiGeorge syndrome. Taken together, these data suggest that altered Wnt5a expression could impact various (cardiac) developmental pathways and makes it a plausible candidate gene at the locus, although the exact mechanisms need to be further elucidated. In contrast to WNT5A, ERC2, the other gene in the TAD, is a less likely candidate gene at the locus based on the lack of expression in the developing mouse heart and the absence of a cardiac phenotype in Erc2 knock-out mice.[66]

We integrated various epigenetic datasets[42] and showed that the 3p14.3 association signal overlaps a putative regulatory element (ie, RE1). We showed that RE1 drives strong expression in mice at E10.5 in both the cardiac OFT myocardium and ventricles, regions that normally also express Wnt5a during this stage of development,[49] and that RE1 interacts with WNT5A. In addition, we established that RE1 binds TBX20 and that Wnt5a expression is increased in cardiac-specific Tbx20 knock-out mice. These observations suggest that RE1 is needed to regulate Wnt5a expression in the developing heart and that the effect of RE1 is attenuated by TBX20 binding. Together with the fact that Wnt5a transgenic mouse models[45,52] and Tbx20 null mice[53] exhibit outflow tract abnormalities, this further supports a possible causal role of WNT5A at the locus and suggests that RE1 may be in some way involved in the underlying genetic mechanism.

Surprisingly, deletion of RE1 in mice did not result in an abnormal phenotype or gene expression changes. Previous studies have shown that regulatory element redundancy is a hallmark of mammalian genomes.[50,51] This means that when an entire regulatory element is absent, as is the case for RE1−/−, redundant elements in the locus compensate, thus providing a means of protection from detrimental phenotypes caused by loss of individual regulatory elements.

We searched for the causal variant within RE1 and demonstrated that all 3 major haplotypes overlapping RE1, showed robust enhancer activity. A SNP on one of these haplotypes (rs79492708) abrogates a TBX20 binding site, and this haplotype significantly increases its activity compared with the other 2 haplotypes. However, given the low r2 (0.07) between rs79492708 and rs56219800 (lead SNP), we concluded this was not the causal SNP driving the association at the 3p14.3 locus. Future studies will be needed to show if the causal SNP is another SNP (or a combination of SNPs) within RE1, although in silico analysis did not identify alteration of putative TF-binding sites other than the TBX20 one. Alternatively, the causal SNP(s) is/are located outside of RE1, although scrutiny of currently available epigenetic data did not identify other putative regulatory elements overlapping the association signal (45 kb, defined by the lead SNP and SNPs with r2≥0.5).

With regards to the many other loci just below the genome-wide significance threshold, many of these loci contain genes that are related to cardiac morphogenesis, some of which have been previously specifically linked to cardiac OFT defects. For example, a risk locus on chromosome 3 is located within the same TAD as TBX3. TBX3 deficiency in mice results in double outlet right ventricle,[67,68] an OFT defect which, similar to D-TGA, results from an alignment defect of the great vessels.[69] A subset of these genes from other loci have also previously been shown to harbor rare variants in families with cardiac OFT defects (eg, NKX2-6[70,71]) or in probands included in large-scale studies of subjects with CHD.[46-48] Some of these genes have been shown to be involved in left-right patterning and ciliopathies (eg, EFCAB1[72] and KIAA0586[73]). Although we excluded patients with known laterality defects in this study, variants in genes associated with heterotaxia have also been identified in patients with D-TGA in the absence of other extra-cardiac defects.[7,8,10,74] Moreover, cilia(-related) genes have been strongly implicated in CHD in mice[75] as well as in a variety of CHD patients.[48,76] Interestingly, the risk locus at 16p13.11 overlaps a rare copy number variant (CNV) identified in a D-TGA patient[54] and is a region where recurrent CNVs have been reported in patients with CHD.[55,77] This CNV includes the gene MYH11, which is primarily known for its involvement in familial thoracic aortic aneurysm, but recessive mutations have also been identified in a patient with D-TGA.[47] Thus, it is likely that at least some of the identified subthreshold loci are bona fide susceptibility loci. Increasing the sample size, although challenging given the rarity of D-TGA, would thus likely uncover additional significantly associated loci that could aid in dissecting the developmental pathways contributing to D-TGA. Future studies will also be needed to show whether these findings can be generalized to patients of other ethnic backgrounds.

The results from this study should be considered in the light of some limitations. Since the majority of patients included in this study did not undergo any other genetic testing, we cannot rule out that some cases will have rare genetic variants or CNVs. However, a recent study in patients with severe neurodevelopmental disorders found no difference in the role of common genetic variants between patients with a likely diagnostic rare variant and those without such a variant.[78] Heritability estimation may cause biased results in ascertained case-control datasets, as is the case in this study. However, we performed stringent genotypic QC and used Phenotype-Correlation-Genotype-Correlation to account for case over-sampling. To minimize stratification bias because of mixed ethnicity or platform differences, we also estimated SNP-heritability in a more homogenous group (Dutch only, all genotyped on the Illumina HumanOmniExpress array), which showed similar results. In addition, this GWAS was restricted to individuals of European descent. Analysis in non-European cohorts will be necessary to explore if the results are applicable to individuals of other ancestries, and to further refine the genome-wide association signal. Despite these limitations, a main strength of this study is that we have collected the largest set of D-TGA patients thus far, which allowed us to identify and replicate a genome-wide significant locus and allowed us to perform extensive in vitro and in vivo experiments demonstrating the functional relevance of this locus.

In conclusion, this work provides the first evidence of a marked contribution of genome-wide common genetic variation to the risk of D-TGA and suggests that dysregulation of WNT5A may partly underlie susceptibility to D-TGA. Increasing knowledge about the genetic cause of D-TGA will contribute to improved genetic counseling of individuals and families affected by CHD and, specifically, D-TGA.

Article Information

Acknowledgments

We thank Sylvia Mantels and Lia C.J.M. Engelfriet-Rijk for their valuable work recruiting patients and managing the CONCOR registry. We thank Iwan Vaandrager and Wouter de Laat from the Hubrecht Institute in Utrecht for their help in generating the virtual 4C figure. The KORA-Study Group consists of A. Peters (speaker), H. Schulz, L. Schwettmann, R. Leidl, M. Heier, K. Strauch, and their co-workers, who are responsible for the design and conduct of the KORA studies.

Sources of Funding

This work has been funded by the Dutch Heart Foundation (CVON project 2014-18 CONCOR-genes) (to B.J. Bouma, C.R. Bezzina, J. Bakkers, M.W.T. Tanck, and V.M. Christoffels), a research grant from the Children’s Heart Foundation (to C.R. Bezzina) and by Foundation Leducq (to V.M. Christoffels). R. Tadros received support from the Canadian Heart Rhythm Society’s George Mines Award, the European Society of Cardiology research award, and the Philippa and Marvin Carsley Cardiology Chair and is currently a clinical research scholar of the Fonds de Recherche du Québec—Santé. F. Bajolle is funded by the Assistance Publique – Hôpitaux de Paris (APHP). G.M. Blue is supported by a Heart Foundation Postdoctoral Fellowship (101894). The work of E.M. Lodder is partly financed by the Dutch Research Council (NWO) through the NWO Talent Scheme VIDI-91718361 and the CVON RESCUED project. M. Müller-Nurasyid is supported by the department of Internal Medicine I (Cardiology), Hospital of the Ludwig-Maximilians-University (LMU) Munich, Munich, Germany. P.J. Shaw is supported as an NIHR Senior Investigator and by the Sheffield NIHR Biomedical Research Centre. E. Oechslin currently holds the Bitove Family Professorship of Adult Congenital Heart Disease. Collection of the samples provided by J. Kline was supported by NIH grant RO1HL-08146 from National Heart, Lung, and Blood Institute (NHLBI) to Dorothy Warburton and Mike Wigler. The sequencing for this work was supported by a grant to the Baylor-Hopkins Center for Mendelian Genomics (5UM1HG006542) from National Human Genome Research Institute (NHGRI) and NHLBI. A.S. Bassett received support for this study from research grants from the Canadian Institutes of Health Research (CIHR) (MOP-93722) and the University of Toronto McLaughlin Centre. A.S. Bassett holds the Dalglish Chair in 22q11.2 Deletion Syndrome at the University of Toronto and University Health Network. R.L.F. van der Palen is supported by The Dutch Heart Foundation (grant number 2014T087). J. Breckpot is supported by a Senior Clinical Investigator Fellowship of the FWO - Flanders. S. Mital is supported by Canadian Institutes of Health Research, under the frame of ERA PerMed, Heart and Stroke Foundation of Canada, Ted Rogers Centre for Heart Research. E. GoldmuntzS received support for this study from research grants from the Eunice Kennedy Shriver National Institute of Child Health and Human Development (P01HD070454); the NHLBI (P50-HL74731), including the Pediatric Cardiac Genomics Consortium (PCGC) (U01-HL098188, U01HL131003, U01-HL098147, U01-HL098153, U01-HL098163, U01-HL098123, U01-HL098162, U01-HL-09003); the NHGRI (U54HG006504); the National Center for Research Resources (M01-RR-000240, RR024134; now the National Center for Advancing Translational Sciences, grant number UL1TR000003). B. Keavney is funded by a British Heart Foundation Personal Chair. D. Bonnet is funded by Association pour la Recherche en Cardiologie du Fœtus à l’Adulte (ARCFA) and Fondation Cœur et Artères. The MINE study (J.H. Veldink) has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n° 772376 – EScORIAL). The collaboration project is co-funded by the PPP Allowance made available by Health~Holland, Top Sector Life Sciences & Health, to stimulate public-private partnerships. The KORA study was initiated and financed by the Helmholtz Zentrum München – German Research Center for Environmental Health, which is funded by the German Federal Ministry of Education and Research (BMBF) and by the State of Bavaria. Furthermore, KORA research was supported within the Munich Center of Health Sciences (MC-Health), Ludwig-Maximilians-Universität, as part of LMUinnovativ. This study makes use of data generated by the Wellcome Trust Case-Control Consortium. A full list of the investigators who contributed to the generation of the data is available from www.wtccc.org.uk. Funding for the project was provided by the Wellcome Trust under award 076113, 085475, and 090355. This work was supported by the Competence Network for Congenital Heart Defects, which has received funding from the Federal Ministry of Education and Research, grant number 01GI0601 (until 2014), and the DZHK (German Centre for Cardiovascular Research; as of 2015). A subset of the data used in the replication cohort were generated by the PCGC, under the auspices of the NHLBI’s Bench to Bassinet Program https://benchtobassinet.com (dbGaP Study Accession: phs001194.v2.p2). The PCGC program is funded by the NHLBI, National Institutes of Health, U.S. Department of Health and Human Services through grants UM1HL128711, UM1HL098162, UM1HL098147, UM1HL098123, UM1HL128761, and U01HL131003. This article was not prepared in collaboration with investigators of the PCGC, has not been reviewed and/or approved by the PCGC, and does not necessarily reflect the opinions of the PCGC investigators or the NHLBI.

Disclosures

None.

Supplemental Material

Supplemental Methods

Figures S1–S9

Tables S1–S9

References[77-103]

Major Resource Table