https://orcid.org/0000-0001-9634-5908
Abstract
The present study aims to characterize the genetic risk architecture of bicuspid aortic valve (BAV) disease, the most common congenital heart defect.
We carried out a genome-wide association study (GWAS) including 2236 BAV patients and 11 604 controls. This led to the identification of a new risk locus for BAV on chromosome 3q29. The single nucleotide polymorphism rs2550262 was genome-wide significant BAV associated (P = 3.49 × 10−08) and was replicated in an independent case–control sample. The risk locus encodes a deleterious missense variant in MUC4 (p.Ala4821Ser), a gene that is involved in epithelial-to-mesenchymal transformation. Mechanistical studies in zebrafish revealed that loss of Muc4 led to a delay in cardiac valvular development suggesting that loss of MUC4 may also play a role in aortic valve malformation. The GWAS also confirmed previously reported BAV risk loci at PALMD (P = 3.97 × 10−16), GATA4 (P = 1.61 × 10−09), and TEX41 (P = 7.68 × 10−04). In addition, the genetic BAV architecture was examined beyond the single-marker level revealing that a substantial fraction of BAV heritability is polygenic and ∼20% of the observed heritability can be explained by our GWAS data. Furthermore, we used the largest human single-cell atlas for foetal gene expression and show that the transcriptome profile in endothelial cells is a major source contributing to BAV pathology.
Our study provides a deeper understanding of the genetic risk architecture of BAV formation on the single marker and polygenic level.
Time of primary review: 34 days
1. Introduction
Bicuspid aortic valve (BAV) is the most common congenital heart defect (CHD) with a prevalence of 0.5–2%.1 BAV is of major clinical relevance because as many as 71% of patients develop aortic valve stenosis (AS) and ∼60% experience aortic root dilatation with risk of aortic dissection. Furthermore, 10–30% of patients with BAV develop infective endocarditis.2 As a consequence, BAV accounts for nearly 50% of all aortic valve replacements.3
Genetic factors are strongly involved in the aetiology of BAV with an observed heritability of 47–89%.4,5 This range of heritability reflects the genetically heterogeneous nature of BAV with a wide spectrum of disease-causing mutations and risk variants. On one end of this spectrum, BAV is the result of highly penetrant mutations in single genes. This refers to monogenic BAV forms, which can be syndromic or non-syndromic. Examples of syndromic BAV that are characterized by additional extra-cardiac conditions include Loeys-Dietz syndrome with mutations in TGFBR1 or TGFBR2 as well as familial thoracic aortic aneurysm and dissection (TAAD) syndrome with mutations in ACTA2.6,7 In contrast, mutations in NOTCH1, ADAMTS19, SMAD6, or ROBO4 were found in monogenic non-syndromic cases of BAV.8–11 However, the vast majority of BAV cases are sporadic12 and result from a complex genetic background.
Genome-wide association studies (GWAS) provides an unbiased approach for discovering genetic risk variants for complex genetic traits and two GWAS involving BAV patients have been published so far. The first GWAS analysed 466 European BAV cases as discovery sample and 1326 European BAV cases as replication cohort.13 This study identified one risk single nucleotide polymorphism (SNP) at the 3′ end of GATA4 with genome-wide significant BAV association as well as one independent SNP within GATA4 that represents a missense variant (p.Ser377Gly) and was suggestive BAV associated. In addition, a rare missense variant (p.Thr1221Met) within DHX38 showed genome-wide significant BAV association. In a second GWAS of AS, the authors used a follow-up sample of 1555 European BAV patients in order to assess whether associated AS risk variants occurred on the background of a tricuspid or bicuspid valve formation.14 Of all genome-wide significant associated AS variants, SNPs at PALMD, TEX41, or MYH6 were also associated with BAV. All association findings along with information on the identified BAV risk variants and used GWAS samples are summarized in Supplementary material online, Table S1.
In the present study, we further characterized the genetic architecture of BAV and present single marker and polygenic association findings from the largest GWAS analysed to date. In addition, we functionally characterized a new BAV GWAS locus in zebrafish using CRISPR/Cas9-mediated knockout and antisense oligonucleotide Morpholino- (MO-) mediated gene knockdown.
2. Methods
2.1 Study subjects
The GWAS sample consisted of two European case–control cohorts that were recruited in Germany and the UK. The German cohort comprised 1539 BAV patients that were recruited by the University Hospitals of Bonn and Lübeck as well as the German Heart Centre Munich and 4414 controls that were part of the population-based Heinz-Nixdorf Recall study.15 The UK cohort comprised 697 BAV patients that were recruited by the University Hospitals of Leicester as well as 7190 controls that were part of the UK Biobank.16 In total, we analysed 2236 BAV patients (583 females, 1653 males) and 11 604 controls (4415 females, 7189 males) for the GWAS. The replication sample consisted of 421 German BAV patients (99 females, 322 males) that were recruited by the University Hospital of Bonn and the German Heart Centre Munich as well as 3643 controls (1266 females, 2377 males) that were part of a blood donor cohort recruited at the University Hospital of Kiel. A detailed description of all BAV patients and controls analysed within this study is given as Supplementary material online, Cohort description.
2.2 Declaration
The authors state that their study complies with the Declaration of Helsinki, that the locally appointed ethics committees have approved all research protocols and that informed written consent has been obtained from all study participants before the inclusion of subjects into the study.
2.3 Genome-wide genotype data, quality control, and imputation
All patient samples (GWAS and replication) were genome-wide genotyped within this study. For the control cohorts (GWAS and replication) genotype data were obtained from previous studies. All genotyping arrays are summarized for each cohort in Supplementary material online, Methods.
The pre-imputation quality control (QC) is described in detail in Supplementary material online, Methods. After the QC, the sets of remaining SNPs were imputed using the TOPMed Imputation Server and TOPMed Reference panel17–19 for the German samples. For the UK samples, the Sanger Imputation Server and a combined reference panel using the Haplotype Reference Consortium release 1.1 (HRC r1.1)20 with additional variants from the UK10K and the 1000 Genomes phase 3 reference panel21 were used. After imputation, we excluded variants with a Hardy-Weinberg equilibrium deviation of P < 10−06 in patients or P < 10−04 in controls, a minor allele count < 20 or a minor allele frequency (MAF) in cases <0.01.
2.4 GWAS and GWAS meta-analysis
For the German cohort, trait associations were calculated with Firth logistic regression as implemented in Plink2,22 conditional on the first 10 variance standardized genotype principal components.
For the British cohort, logistic models were run in Plink222 utilizing Firth regression where models failed to converge. Analyses were adjusted for the first five genetic principle components.
The meta-analysis was performed with METAL23 using the fixed-effect inverse standard error-weighting approach, with a standard genome-wide significant threshold of 5 × 10−08.
To investigate whether independent associations are present in regions of genome-wide significance we carried out conditional analyses on the strongest associated SNP in each region with the meta-analysis approach COJO from GCTA.24,25 Here, we used the German cohort as a linkage disequilibrium (LD) reference sample. Additionally, we performed reciprocal conditional analysis with rs2293232 and rs2550262 on chromosome 3q29 as implemented in Plink2.22
2.5 LD score regression analyses
SNP-based heritability was calculated using LD score regression (LDSC).26 The reference files from GRCh37 were used as advised by the authors.
We next applied LDSC to partition the SNP-based heritability for BAV across transcriptome profiles from 16 foetal cardiac cell types that characterize heart development and were available through DESCARTES,27 the largest human single-cell atlas for foetal gene expression. Details of this analysis are provided as Supplementary material online, Methods.
2.6 FUMA analysis
We used FUMA (Functional Mapping and Annotation of Genome-Wide Association Studies)28,29 via the web-interface to gain functional insights into BAV. FUMA was used with default settings for the SNP2GENE and the GENE2FUNC modules. These tools provide information about functional annotations of associated SNPs and the possible biological role of the annotated genes.
2.7 Expression analysis at the BAV risk locus on chromosome 3q29 in embryonic/foetal and adult heart
Transcriptome data from single-cell RNA sequencing (scSeq) experiments derived from seven embryonic/foetal30 and two adult hearts31 were used. A detailed description of the tissue samples and scSeq protocols can be found in the corresponding publications.30,31 In brief, raw data of 676 single-cell transcriptomes isolated from seven human embryonic/foetal hearts (4.5–10 weeks of gestation) and three heart compartments (proximal outflow tract, ventricle, atria) were downloaded from the NCBI Sequence Read Archive (accession no. PRJNA510181). In addition, in-house data of 18 485 single-nuclei transcriptomes isolated from two adult human heart ventricles and atria were used. All transcriptomes were subjected to a standardized QC procedure followed by alignment to the Genome Reference Consortium Human Build 38 (hg38) using RNA STAR32 on the Galaxy web platform.33 This was followed by deduplication (BroadInstitue Picard tools: MarkDuplicates. https://broadinstitute.github.io/picard/Accessed August 15, 2020) and count matrix creation34 for single-cell data as well as the Cellranger pipeline35 with default settings including the –include-introns option for single-nuclei data, respectively. Seurat36 objects for all samples were created for downstream analyses from filtered count matrices. After quality filtering, the data were normalized and scaled and variable features for 669 embryonic/foetal and 16 896 adult single cells were detected using SCTransform.37 Data from embryonic and adult cardiac tissue were integrated using integration anchors and label transfer as described by Stuart et al.38 PCA and uniform manifold approximation and projection for dimension reduction were used to cluster cells into distinct biological identities and cell type identification was established based on the expression of known markers. For expression analysis of APOD, MUC20, MUC4, and TNK2, the ‘RNA’-assay of the integrated Seurat object was split into adult and embryonic cardiac cell populations retaining the clustering information of the ‘integrated’ -assay. The two populations were each normalized (NormalizeData) and the Seurat command FeaturePlot was used for visualization of gene expression with settings min.cutoff = ‘q10’ and max.cutoff = ‘q90’.
2.8 Zebrafish experiments
Handling of zebrafish was done in compliance with the guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes and with German and Brandenburg state law, carefully monitored by the local authority for animal protection (LUGV, Brandenburg, Germany). Zebrafish embryos were euthanized by rapid chilling upon submersion into ice water bath (five parts ice/one part water, 2–4°C) for at least 20 min to ensure death. The Tg(kdrl:EGFP)s84339 zebrafish line was used and maintained under standard conditions as previously described.40 MOs (Gene Tools, LLC, USA) were diluted in ddH2O to 1 mM stock solution and injected into the yolk at the one-cell stage in 1 nL volume. Details of gRNAs, CRISPR/Cas9, and MOs sequences as well as experimental conditions, whole-mount immunohistochemistry, and whole-mount in situ hybridization protocol are provided in the Supplementary material online, Methods.
3. Results
3.1 GWAS single-marker association analysis
Our combined GWAS sample consisted of 2236 patients and 11 604 controls (see Methods for details). Genetic variants at three loci were associated with BAV at genome-wide significance. Two of the loci—on chromosome 1p21 and 8p23—have previously been reported,13,14 whereas the locus on chromosome 3q29 is novel. The findings for the lead variants at each locus are summarized in Table 1. The corresponding Q-Q- as well as Manhattan-plots are shown in Supplementary material online, Figures S1 and S2.
Lead associations of genome-wide significant and replicated BAV risk loci
| SNP | Chromosome [position in bp (hg38)] | Effect allele/other allele | Frequency in % | P-value | OR | 95% CI |
|---|---|---|---|---|---|---|
| rs11166276 | 1 (99 579 683) | T/C | 0.50 | 3.97 × 10−16 | 1.33 | 1.29–1.38 |
| rs13408842 | 2 (145 086 078) | A/G | 0.48 | 7.68 × 10−04 | 1.13 | 1.09–1.17 |
| rs2293232 | 3 (195 770 272) | T/C | 0.13 | 1.96 × 10−08 | 1.33 | 1.27–1.40 |
| rs2550262 | 3 (195 763 273) | C/A | 0.30 | 3.49 × 10−08 | 1.30 | 1.24–1.36 |
| rs6996406 | 8 (11 945 922) | A/G | 0.40 | 1.61 × 10−09 | 1.61 | 1.49–1.73 |
| rs3729856 | 8 (11 757 066) | G/A | 0.15 | 1.71 × 10−02 | 1.12 | 1.07–1.18 |
| SNP | Chromosome [position in bp (hg38)] | Effect allele/other allele | Frequency in % | P-value | OR | 95% CI |
|---|---|---|---|---|---|---|
| rs11166276 | 1 (99 579 683) | T/C | 0.50 | 3.97 × 10−16 | 1.33 | 1.29–1.38 |
| rs13408842 | 2 (145 086 078) | A/G | 0.48 | 7.68 × 10−04 | 1.13 | 1.09–1.17 |
| rs2293232 | 3 (195 770 272) | T/C | 0.13 | 1.96 × 10−08 | 1.33 | 1.27–1.40 |
| rs2550262 | 3 (195 763 273) | C/A | 0.30 | 3.49 × 10−08 | 1.30 | 1.24–1.36 |
| rs6996406 | 8 (11 945 922) | A/G | 0.40 | 1.61 × 10−09 | 1.61 | 1.49–1.73 |
| rs3729856 | 8 (11 757 066) | G/A | 0.15 | 1.71 × 10−02 | 1.12 | 1.07–1.18 |
The associations are shown for the risk alleles (effect alleles). Mean allele frequencies for the risk alleles across all GWAS cohorts, P-values, odds ratios (OR), and the corresponding 95% confidence intervals (CIs) are shown.
Lead associations of genome-wide significant and replicated BAV risk loci
| SNP | Chromosome [position in bp (hg38)] | Effect allele/other allele | Frequency in % | P-value | OR | 95% CI |
|---|---|---|---|---|---|---|
| rs11166276 | 1 (99 579 683) | T/C | 0.50 | 3.97 × 10−16 | 1.33 | 1.29–1.38 |
| rs13408842 | 2 (145 086 078) | A/G | 0.48 | 7.68 × 10−04 | 1.13 | 1.09–1.17 |
| rs2293232 | 3 (195 770 272) | T/C | 0.13 | 1.96 × 10−08 | 1.33 | 1.27–1.40 |
| rs2550262 | 3 (195 763 273) | C/A | 0.30 | 3.49 × 10−08 | 1.30 | 1.24–1.36 |
| rs6996406 | 8 (11 945 922) | A/G | 0.40 | 1.61 × 10−09 | 1.61 | 1.49–1.73 |
| rs3729856 | 8 (11 757 066) | G/A | 0.15 | 1.71 × 10−02 | 1.12 | 1.07–1.18 |
| SNP | Chromosome [position in bp (hg38)] | Effect allele/other allele | Frequency in % | P-value | OR | 95% CI |
|---|---|---|---|---|---|---|
| rs11166276 | 1 (99 579 683) | T/C | 0.50 | 3.97 × 10−16 | 1.33 | 1.29–1.38 |
| rs13408842 | 2 (145 086 078) | A/G | 0.48 | 7.68 × 10−04 | 1.13 | 1.09–1.17 |
| rs2293232 | 3 (195 770 272) | T/C | 0.13 | 1.96 × 10−08 | 1.33 | 1.27–1.40 |
| rs2550262 | 3 (195 763 273) | C/A | 0.30 | 3.49 × 10−08 | 1.30 | 1.24–1.36 |
| rs6996406 | 8 (11 945 922) | A/G | 0.40 | 1.61 × 10−09 | 1.61 | 1.49–1.73 |
| rs3729856 | 8 (11 757 066) | G/A | 0.15 | 1.71 × 10−02 | 1.12 | 1.07–1.18 |
The associations are shown for the risk alleles (effect alleles). Mean allele frequencies for the risk alleles across all GWAS cohorts, P-values, odds ratios (OR), and the corresponding 95% confidence intervals (CIs) are shown.
To assess the robustness of our findings, we first compared our data with previously published GWAS findings. SNP rs11166276 on chromosome 1p21 near PALMD showed the strongest genome-wide significant BAV association in our study [P = 3.97 × 10−16; odds ratio (OR) = 1.33] (Table 1 and Figure 1A). The risk variant is in perfect LD to rs7543130 (r2 = 1), which was the lead associated SNP in the GWAS of Helgadottir et al.14 On chromosome 8p23, rs6996406 showed genome-wide significant BAV association in our study (P = 1.61 × 10−09, OR = 1.61) (Table 1 and Figure 1B). The risk variant is in perfect LD to rs118065347 (r2 = 1), which was the lead SNP in the GWAS of Yang et al.13 and was linked to GATA4 as disease-contributing gene. Furthermore, we found nominal significant association at two additional BAV risk loci that have been reported previously.13,14 On chromosome 2q22, SNP rs13408842 within Intron 4 of TEX41, which was the lead variant at this locus in the GWAS of Helgadottir et al.,14 showed BAV association in our sample (P = 7.68 × 10−04, OR = 1.13) (Table 1). In addition, the missense variant rs3729856 (p.Ser377Gly) within Exon 6 of GATA4 on chromosome 8p23 that showed independent disease association at this locus in the GWAS of Yang et al.13 was significant BAV associated in our sample (P = 1.71 × 10−02, OR = 1.12) (Table 1). Thus, all common risk loci at PALMD, TEX41, and GATA413 that have been previously been found in BAV GWAS were replicated in our data set. However, we could not test the previously described risk SNPs rs137867582 (p.Thr1221Met) within DHX3813 and rs387906656 (p.Arg721Trp) within MYH6,14 because the MAFs of the corresponding risk alleles are extremely rare in the European population (<0.001% according to gnomAD41).
![Regional association plots of BAV risk loci on chromosome 1p21 near PALMD (A), 8p23 near GATA4 (B) and 3q29 within MUC4 (C). At each locus, associations [−log10(P-values)] are shown for SNPs flanking 400 kb on either side of the lead associated SNP. The lead variant is indicated by the corresponding rs-number. Other markers at each locus are displayed by different colours, which indicate different levels of LD (r2) to the lead SNP. Furthermore, all annotated genes within each region are shown with arrows indicating their transcription direction.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/cardiovascres/119/3/10.1093_cvr_cvac099/1/m_cvac099f1.jpeg?Expires=1747892646&Signature=FZVKKusSN19WUubJwVwH-flDqx7pqsMrHQIp4vtSeiqrCyIbIF2BRMUe9YyXQlkk5E3mJa7wSRKSJAfPArJNGQ6lRhnKtb1ujFXuPzfu5bB7SFv4BUnNT7te1EexkOokwY4uoaYb7-UnDcoMmWa4EWexzWkL-1TLEqtVehiEdJaFI2z4Ub6UaC8xbS0TG~ZWHYmD~hLPgq7dIvJ1fbFDvnxNfbjmNCAhh0kS-HBA3IwQ1zJeAy3t6C~IJoN8wATY5lK2L7YObjna5HWPjNJyIdEDkeT-4DYpZQtIy-Tqa9ehI10A826iOaT-yv7Gq6YoxBMBE-rPP0fGC7sbB-YQ5g__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Regional association plots of BAV risk loci on chromosome 1p21 near PALMD (A), 8p23 near GATA4 (B) and 3q29 within MUC4 (C). At each locus, associations [−log10(P-values)] are shown for SNPs flanking 400 kb on either side of the lead associated SNP. The lead variant is indicated by the corresponding rs-number. Other markers at each locus are displayed by different colours, which indicate different levels of LD (r2) to the lead SNP. Furthermore, all annotated genes within each region are shown with arrows indicating their transcription direction.
At the novel risk locus on chromosome 3q29 rs2293232 showed genome-wide significant BAV association (P = 1.96 × 10−08, OR = 1.33) (Table 1 and Figure 1C). This SNP is located in Exon 5 of MUC4 and represents a missense variant (p.Ala4448Thr) that is predicted to be benign according to gnomAD.41 In addition, a second and independent variant at this locus showed genome-wide significant BAV association (P = 3.49 × 10−08, OR = 1.30) (Table 1). This SNP, rs2550262, is only in weak LD to rs2293232 (r2 = 0.53). Accordingly, both variants showed significant BAV association upon reciprocal conditioning (Prs2293232 = 4.96 × 10−02, Prs2550262 = 1.12 × 10−03). Of note, rs2550262 is in nearly perfect LD to rs2246901 (r2 = 0.97), a missense variant in Exon 13 of MUC4 (p.Ala4821Ser) and deleterious according to gnomAD.41 Thus, this variant represents the most likely causal SNP at this locus.
Finally, we analysed 421 independent European BAV patients and 3643 ethnically matched controls in order to replicate the newly identified risk locus on chromosome 3q29. Here, we found independent BAV association for rs2550262 (P = 1.36 × 10−02, OR = 1.22), one of the lead GWAS SNPs within MUC4. In addition, rs2246901, which is in nearly perfect LD to rs2550262 and represents the most likely causal SNP at this locus, showed significant BAV association (P = 1.31 × 10−02, OR = 1.22). In the combined sample (GWAS and replication) rs2550262 was BAV associated with P = 1.85 × 10−09 (OR = 1.28) (see Supplementary material online, Table S2).
3.2 Genetic BAV architecture on the polygenic level
We used the LDSC method, which allows the collective analysis of common GWAS variants, to estimate the SNP-based heritability of BAV and to test whether GWAS risk SNPs are enriched in regions surrounding genes with specific expression in different cell types of the foetal heart.
Based on our genome-wide data, we estimate the SNP-based heritability of BAV to be 9.55 ± 3.97% standard deviation (SD). The reported BAV heritability, which is based on family studies, ranges from 47 to 89%4,5 with the lower estimate referring to complex genetic BAV forms, while the upper estimate reflects monogenic forms of BAV. Thus, our data already explain 20% of the reported heritability (h2 = 47%).
Next, we used DESCARTES27 to assess cell type-specific GWAS enrichment. This database provide an annotation of cell type-specific transcriptome profiles for 15 organs from 121 foetal tissues. Within the foetal heart, DESCARTES identified 16 cardiac cell types with distinct transcriptome profiles during development. We applied LDSC to partition the SNP-based heritability for BAV across these 16 transcriptome profiles and found nominally significant GWAS enrichment (P < 0.05) among 11 foetal cardiac cell types (see Supplementary material online, Table S3). Of note, we observed only a weak overlap of expressed genes between these cell types (see Supplementary material online, Figure S3). Of the implicated cell types, the transcriptome profile in vascular endothelial cells also showed significant enrichment for BAV GWAS signals after a Bonferroni-correction for multiple testing (Puncorrected = 5.9 × 10−04, Pcorrected = 9.5 × 10−03).
In addition to the LDSC method, we analysed the GWAS data on the polygenic level using FUMA,28,29 a tool that uses functional and biological information as well as GWAS data to prioritize disease genes. The FUMA gene-based analysis was performed for a total of 18 830 genes. Here, two genes were implicated after correction for multiple testing (P < 0.05/18 830 or P < 2.66 × 10−06) (see Supplementary material online, Figure S4). One refers to MUC4, the BAV candidate gene on chromosome 3q29 identified by the GWAS, where FUMA showed significant BAV association with P = 9.66 × 10−08. The other gene that showed significant BAV association in FUMA is DCAF4 on chromosome 14q24 (P = 1.54 × 10−06). In contrast, all FUMA tissue expression profile analyses did not reveal any significant or biologically interesting results (data not shown), which could be because FUMA uses information only from adult heart tissue.
3.3 Expression of genes at the BAV risk locus on chromosome 3q29 in embryonic/foetal and adult heart
We used transcriptome data from scSeq experiments derived from seven embryonic/foetal and two adult hearts30,31 of which expression data were available. Supplementary material online, Figure S5A shows that the cell type-specific embryonic/foetal and adult cardiac transcriptome data lead to perfectly superimposed clusters, which have been published previously.31 We then focused on expression data of genes from the embryonic/foetal and adult heart, in which variants showed at least weak LD (r2 > 0.2) to the strongest BAV-associated SNP on chromosome 3q29. According to Supplementary material online, Figure 1B, these are the genes APOD, SDHAP2, MUC20, MUC4, and TNK2. Of them, all—except SDHAP2—showed expression in embryonic/foetal and adult cardiac cell types, while MUC4 represented the gene that—compared with the other genes—was expressed in embryonic/foetal cardiac cells but almost undetectable in the adult heart (see Supplementary material online, Figure S5C). Thus, also based on these expression data MUC4 represents an interesting BAV candidate gene at the newly identified GWAS locus, although the other genes cannot be excluded as the true risk conferring BAV genes.
3.4 Cardiac valvulogenesis is delayed upon loss of Muc4 in zebrafish
In the GWAS and replication study, SNP rs2550262 in MUC4 was significantly BAV associated and in strong LD to the missense variant rs2246901 (p.Ala4821Ser) that is deleterious according to gnomAD.41 In addition, MUC4 showed almost specific expression in embryonic/foetal cardiac cells compared with adult heart tissue. We, thus, selected MUC4 for first functional testing of its potential involvement on valve leaflet morphogenesis in zebrafish.
Previous transcriptome data on cardiac tissue revealed that the predicted zebrafish orthologous gene si:ch73-105b23.6 (here called muc4) is expressed at 21.5, 54, and 72 h post-fertilization (hpf) during stages of atrioventricular (AV) valvulogenesis.42–44 To functionally study the role of Muc4 during valvulogenesis, we used two different loss-of-function strategies. First, we used a CRISPR/Cas9-mediated knockout (Figure 2) that was based on injecting of a mixed pool of four gRNAs and Cas9 protein. We tested the efficacy of the CRISPR/Cas9-mediated gene knockout on a total of 58 embryos. Genomic PCR of the muc4 locus revealed that the mixture of four gRNAs had induced mutations in almost all embryos. Among these, we chose 15 embryos that had not only a large deletion at the muc4 locus but also lacked a WT allele (Figure 2E and E’). When we characterized these 15 muc4 crispant embryos, we found that the gene knockout had not caused any gross defects in cardiac morphogenesis at 58 hpf. However, careful analysis of endocardial cell ingression into the cardiac jelly during the morphogenesis of the AV valve leaflets revealed a delay in muc4 crispants compared with wild-type (WT) non-injected control embryos (Figure 2A–D). Whereas in WT endocardial cell ingression has produced an abluminal layer of endocardial cells at that stage (Figure 2A’, asterisks), muc4 crispants exhibited a delayed endocardial cell ingression into the cardiac jelly. We identified three phenotype classes among the 15 muc4 crispants as quantified in Figure 2B and D. As an alternative approach, we also designed three non-overlapping MO antisense oligonucleotides to target muc4 pre-mRNA splice sites (see Supplementary material online, Figure S6). For this alternative MO knockdown approach (see Supplementary material online, Figure S6), a similar delay in endocardial cell ingression during AV valvulogenesis was observed (as quantified in see Supplementary material online, Figure S7). Thus, the muc4 crispant and morphant AV valvular phenotypes correlate with a temporal delay in cardiac valve development.
Zebrafish muc4 knockout using CRISPR/Cas9. (A–D) Maximum intensity projections from confocal z-stack images of hearts at 58 hpf. (A’–D’) Single confocal z-section images of the AVC region. Endocardial cells in luminal positions (arrowheads) at the AVC are marked by kdrl:GFP and Alcam.45 Endocardial cells that have already invaded the extracellular matrix and are in abluminal positions are marked with asterisks. (A’) While in WT control embryos many endocardial cells are in abluminal positions, some of the 15 muc4 crispants with a complete truncation of the genomic locus exhibit a delayed ingression of endocardial cells into the cardiac jelly. The severity of the phenotypes varied according to the following classes: (B’) Class 1 embryos (3 of 15) lack an endocardial cell ingression by that stage; (C’) Class 2 embryos (4 of 15), have only a single endocardial cell ingressing into the cardiac jelly at that stage; (D’) in muc4 crispants of class 3 (8 of 15), a few endocardial cells have ingressed into the cardiac jelly. V, ventricle; A, atrium; AVC, atrioventricular canal. Scale bars are 20 µm (A–D) and 10 µm (A’–D’). (E) Schematic representation of the zebrafish muc4 locus and gRNA binding sites (red arrows). Black arrows indicate the positions of primers used to assess the efficacy of CRISPS-induced knockout by PCR. (E’) Single embryo PCR products from control embryos (ctrl) and embryos injected with a mixture of four gRNAs and Cas9 protein (muc4_crispants). Red boxes are indicating samples that lack any WT genomic amplification bands. Among 58 injected embryos that were genotyped, 15 had a large truncation of the genomic region and did not contain any WT band. Only these crispants were used for phenotypic characterization. MW, molecular weight.
4. Discussion
Using the largest GWAS sample analysed so far, we independently confirmed all common loci that have previously been reported to contribute to the risk for BAV disease. In addition, we identified a novel risk locus on chromosome 3q29 and revealed that MUC4 at this locus might be involved in embryonic valvulogenesis. Finally, our results underline the polygenic nature in BAV development.
Although recent GWAS have identified risk variants for BAV development near PALMD (chromosome 1p21) and TEX41 (chromosome 2q22), functional data are still lacking to determine the causal BAV genes at both loci. Helgadottir et al.14 used chromatin interaction maps from different foetal and adult heart tissues and found that the BAV risk variants at both loci are located within regulatory regions that interact with promoters on neighbouring genes. On chromosome 1p21 this involved promoters of PALMD, SNX7, PLPPR5, and PLPPR4 as well as two non-coding RNAs. On chromosome 2q22 promoters of ZEB2 and GTDC1 as well as three non-coding RNAs, including TEX41 have been implicated. In contrast, GATA4 on chromosome 8p23 is so far the only gene at a GWAS locus with convincing evidence of being involved in BAV development. Common risk variants and rare deleterious mutations in GATA4 have been repeatedly found in various CHDs, including BAV, atrial and ventricular septal defect as well as Tetralogy of Fallot.46–51
In the present study, we identified and independently confirmed a fourth BAV risk locus on chromosome 3q29. Based on our GWAS data, expression data from embryonic/foetal and adult heart tissues as well as zebrafish knockout and knockdown experiments, MUC4 represents an interesting BAV candidate gene at this locus. SNP rs2246901 in Exon 13 of MUC4 is in nearly perfect LD to the lead-associated GWAS variant and represents a missense variant (p.Ala4821Ser) that is deleterious according to gnomAD.41 In line with this finding, knockout and knockdown of this gene in zebrafish revealed that muc4 may impact cardiac AV valve leaflet morphogenesis causing a developmental delay of the multi-layering of valve leaflets. The zebrafish heart is two-chambered containing three sets of valves (AV, inflow, and outflow tract valves). In our study, we focused on the development of the AV valve between atrial and ventricular chambers because the development of this valve has been well characterized and it forms at an early developmental stage between 48 and 72 hpf. Of note, zebrafish outflow tract valve leaflets, which are related to human aortic valves affected in BAV, emerge and elongate between 96 and 144 hpf and, currently much less is known about their mode of development.52 We also found that by this stage that AV valve development had partially recovered in our Muc4 knockout and knockdown zebrafish models. On the cellular level, MUC4 also represents an interesting BAV candidate gene. MUC4 is involved in epithelial-to-mesenchymal transformation (EMT),53,54 a cellular process that is of major relevance for cardiac valve development.55–57 During EMT, the molecular phenotype of the endocardial epithelial cells changes towards a mesenchymal character and those cells migrate from the inner lining of the endocardium towards the cardiac jelly forming the so-called endocardial cushions.58 At later stages fusion of these cardiac cushions results in the formation of valve leaflets.58 Interestingly, induced molecular changes in the composition of the prevalvular mesenchyme can result in the formation of BAVs.59 A similar transition of cells in adult life can be seen during cancer and especially in tumour metastasis. In this context expression of MUC4 has been previously shown to initiate the process of EMT in cancer cells,53 whereas inhibition of MUC4 expression resulted in suppressed cell invasion and EMT.54
Although our genetic and zebrafish data indicate MUC4 as an interesting BAV candidate gene, Muc4 knockout mice are viable, fertile, and show no obvious developmental abnormalities.60,61 However, mild cardiac valve abnormalities may have been missed in the latter studies. In muc4 knockout and knockdown zebrafish models, the developmental phenotype only became apparent after focused cardiac phenotyping.
In addition, our study represents the first that examined the genetic BAV architecture beyond the single-marker level. Our findings show that a substantial fraction of BAV heritability is polygenic with an estimated SNP-based heritability of 9.55 ± 3.97% (SD) that already explains around 20% of the observed BAV heritability. In addition, we found that the transcriptome profile of foetal endothelial cells in DESCARTES27 showed significant enrichment of GWAS signals. This implies that endothelial cells, which along with interstitial cells are the most common aortic valve cell type,62 are of major relevance for BAV development. However, the transcriptome profiles of 10 additional foetal heart cell types showed nominally significant enrichment of BAV GWAS signals. This might point to a complex interaction of different cell types involved in BAV development and might also explain why BAV patients often present with other CHDs, where additional disease-relevant cell types are involved.
Our study has two main limitations. Most of our patients were symptomatic cases, mainly with an AS phenotype. Even though the vast majority of BAV patients will become symptomatic, genetic variability also influences the BAV outcome. We, thus, cannot exclude the possibility that we have also identified genetic risk factors that influence the occurrence of BAV complications rather than its development. In addition, our study did not allow the identification of rare genetic risk factors contributing to BAV. We, thus, could not examine the contribution of the previously described rare risk variants in DHX38 and MYH6. Large exome-wide sequencing studies will provide corresponding data in the future.
5. Conclusions
We independently confirm that the risk loci at PALMD, TEX41, and GATA4 contribute to BAV development. In addition, our results imply that a disturbed MUC4 function contributes to BAV development. The BAV-associated missense variant rs2246901 (p.Ala4821Ser) in MUC4 is predicted to be deleterious and Muc4 knockout as well as knockdown in zebrafish led to a persistent delay in AV valvular development. However, further studies are needed to determine the relevance of MUC4 in BAV development. The present study constitutes the basis for these future studies. Furthermore, we show that a substantial fraction of BAV heritability is polygenic and around 20% of the observed heritability can be explained by our GWAS data. In addition, our findings point to the transcriptome profile in foetal endothelial cells as a major source contributing to BAV pathology.
Supplementary material
Supplementary material is available at Cardiovascular Research online.
Authors’ contribution
Design of the work, acquisition/analysis/interpretation of data for the work, drafting the manuscript, final approval of the manuscript, being accountable for the parts of the work she/he has done: J.G., A.S., R.D., F.F., D.S., M.K., J.-M.S., G.N., M.M.N., A.P.B., S.A.-S., N.J.S., J.E., T.T., J.S. Acquisition of data/study participants for the work, final approval of the manuscript, being accountable for the parts of the work she/he has done: C.P.N., B.A.-K., A.-S.G., C.M.H.B., C.F.Z., L.L.K., P.S.B., T.R.W., S.H., S.E., B.F., S.A.M., M.S., H.K., M.S., F.A.K., P.N., L.B., M.K., V.V., M.A., S.B., K.-L.L., Y.H., M.K., U.M., L.O.C., I.B., C.L., P.U., W.-K.K., J.K., D.E., U.N.-G., P.H., F.W., S.D., H.L., M.D., M.v.S., K.K., T.K., C.H., H.S.
Acknowledgements
The authors thank all study participants who contributed to this study.
Funding
This work was supported by the following grants from the Deutsche Forschungsgemeinschaft (DFG): TRR259, SFB958, SE2016/13-1, SE2016/17-1, NO246/17-1, KR3985/12-1, EXC 2167-390884018, INST 336/104-1, and INST 336/114-1 FUGG. The Leicester samples were collected as part of the BRAVE study and supported by the NIHR Leicester Biomedical Research Centre. R.D. was supported by the Leicester Biomedical Research Centre. C.P.N. and T.R.W. are funded by the British Heart Foundation.
Data availability
The data that support the findings of this study are available on request from the corresponding author (J.S.).
References
Author notes
These authors equally contributed to the work.
Conflict of interest: None declared.
