Identification of CAVS genomic loci

A GWAS meta-analysis for CAVS was performed in 14,819 cases (61% men) and 927,044 controls of European ancestry from six cohorts (Supplementary Fig. 1, Supplementary Fig. 2, and Supplementary Data 1). In total, 17,130,226 variants were included in the meta-analysis, 7,159,971 of which were available in all six cohorts. The genomic inflation factor was 1.099 and the LD score regression intercept was 1.020. Thirty-two loci reached genome-wide significance (P < 5 × 10^−8), including 20 novel loci (Table 1, Fig. 1 and Supplementary Data 2). Among the novel associations, variants located in LPL and LDLR, two major regulators of circulating blood lipids, were identified. The 32 lead variants at the genome-wide associated loci had a concordant direction of effect in the recent meta-analysis by Chen et al.¹⁴. Thirty showed nominal association with CAVS (P < 0.05), of which 29 remained significant when using a threshold of false discovery rate <5%. In the deCODE cohort (n = 2457 CAVS cases and 349,342 controls)¹¹, 28 out of the 32 lead variants had the same direction of effect (Supplementary Data 3). Variants at seven other loci previously reported in European populations were nominally associated (P < 0.05) with a concordant direction of effect, including six that remained significant when using a threshold of false discovery rate <5% (Supplementary Data 4).

The nearest gene at each genome-wide significant locus is indicated, in black for known loci and in red for novel loci. The association of each variant with calcific aortic valve stenosis was obtained from an inverse-variance weighted fixed-effect meta-analysis combining the effect per allele in the cohorts with available data. A p-value below 5 × 10^−8 was considered significant (genome-wide threshold). The genomic inflation factor was 1.099 and the LD score regression intercept was 1.020. The quantile-quantile plot (inset) illustrates the distribution of p-values for each variant tested.

Variant annotation and prioritization

Among variants in linkage disequilibrium with independent significant SNPs, there was a significant enrichment for intronic and exonic functional annotations (Supplementary Data 5). Among the lead SNPs, two were missense variants located in STARD9 (rs28744617) and ASCC2 (rs61736786). Missense variants in ALDH1A2 (rs4646626), MUC4 (rs2246901), and CDAN1 (rs4265781) were in high LD (r²â€‰â‰¥â€‰0.8) with a lead SNP. The missense variants in ASCC2, ALDH1A2 and MUC4 had Combined Annotation Dependent Depletion (CADD) scores of 14.52, 17.37, and 21.50, respectively, ranking them in the 5% most deleterious variants in the genome. A nonsense variant in LPL previously associated with circulating triglyceride levels, rs328 (S447X), was in high LD (r²â€‰= 0.958) with the lead variant at this locus (Supplementary Data 6). The 95% credible set, determined using a Bayesian approach implemented in CAVIAR, included five or less variants for five loci (nearest genes: PALMD, PRRX1, TMEM44, LPA, PDE3A and HMGA2). The complete sets for each locus are reported in Supplementary Data 7. Conditional analyses led to the identification of independent genome-wide significant signals at the LPA (n = 3) and IL6 (n = 2) loci (Supplementary Data 8).

Gene mapping

We mapped the protein-coding genes located within 10 kilobases (kb) of a genome-wide significant SNP and the genes located on either side for lead intergenic SNPs, identifying a list of 79 nearby genes (Supplementary Data 2). A gene-based analysis using MAGMA identified 158 genes (117 additional genes) at a false discovery rate threshold <5%, for a total of 196 prioritized genes (Supplementary Data 9).

Expression in human aortic valves

To identify the genes most likely to be involved in CAVS, we completed extensive transcriptomic profiling of the most relevant tissue, the aortic valve. We performed RNA-sequencing on 500 human aortic valve samples of various disease stages (81.6% with CAVS) and valve morphology (44% bicuspid) from patients recruited in our institution to create the QUEBEC-CAVS-RNA dataset (Supplementary Data 10). First, the expression level of the 196 genes prioritized by positional mapping or MAGMA was evaluated. Among the 187 genes with available gene expression quantification in aortic valves, 22 had a median expression level above the 90^th percentile of all the protein-coding genes, including seven genes located nearby the meta-analysis lead SNPs: RPLP2, CCND1, ACTR2, TALDO1, PDGFRA, PRRX1, and UQCR10 (Supplementary Data 11).

To identify genes for which there is specific expression in the aortic valve, we compared gene expression with 43 GTEx tissues by calculating expression specificity scores (ESS). Among the genes of interest, nine had an ESS > 0.1 for the aortic valve, which corresponded approximatively to the top 5% of the empirical distribution for all the protein-coding genes in the aortic valve (Supplementary Data 12). Five of them were located nearby the meta-analysis lead SNPs: CHST6, CCND1, ALDH1A2, ADAMTS7 and PRRX1. The maximum ESS was observed in the aortic valve for all five genes except PRRX1, for which there was a slightly higher ESS in fibroblasts in GTEx. Other genes with a high ESS in a relevant tissue included HMGA2 (ESS = 0.80 in fibroblasts), NKX2-5 (ESS = 0.42 in atrial appendage), PDE3A (ESS = 0.14 in tibial artery) and PDGFRA (ESS = 0.13 in fibroblasts) (Supplementary Data 12 and Supplementary Fig. 3).

Expression quantitative trait loci

Genome-wide expression quantitative trait loci (eQTL) analyses in 484 human aortic valves identified 4,671,347 significant SNP-gene pairs. Among the 32 meta-analysis lead SNPs, 48 significant SNP-gene pairs were identified in the aortic valve, located in 20 loci (Supplementary Data 13). For these 32 lead SNPs, there was a significantly higher proportion of significant SNP-gene pairs (48/915, 5.25%) compared to all the SNP-gene pairs tested (4,671,347/182,925,823, 2.55%, P = 4.25 × 10^−7). A Wilcoxon rank sum test also confirmed an enrichment towards stronger associations for these 915 SNP-gene pairs (P = 1.71 × 10^−7) (Supplementary Fig. 4).

In 43 tissues from GTEx, significant eQTLs were found for 25 loci (Supplementary Data 14). The number of significant SNP-gene pairs ranged from one (brain anterior cingulate cortex) to 37 (tibial nerve). Taken together, the most significant eQTL was in the aortic valve for 12 loci, including four genes for which there was no reported eQTL in the GTEx tissues: rs6702619-PALMD (P = 7.1 × 10^−119), rs7804522-TWIST1 (P = 5.7 × 10^−22), rs1965668-NKX2-5 (P = 7.9 × 10^−12) and rs72854462-TEX41 (P = 7.8 × 10^−6). For 16 other loci, the most significant eQTL was found in a GTEx tissue. However, these were distributed among ten different tissues, which each had between one and three of the most significant SNP-gene pairs for a given locus (Supplementary Fig. 5). Among these, associations include rs1706003-ATP13A3 in the left ventricle (P = 1.7 × 10^−26), rs771264-RNF144A in fibroblasts (P = 1.9 × 10^−12) and rs11330858-MECOM in the aorta (P = 7.5 × 10^−5).

Transcriptome-wide association study, colocalization and Mendelian randomization

A transcriptome-wide association study (TWAS) was performed to identify genes for which genetically predicted expression in the aortic valve is associated with CAVS. Thirty-five genes reached statistical significance at a false discovery rate <5%, including PALMD and NAV1, identified in a previous TWAS⁹ (Supplementary Data 15). To obtain further evidence on the potential causal role of these genes in CAVS, colocalization and Mendelian randomization (MR) analyses were performed. After removing genes with a low probability for colocalization of the GWAS and eQTL signals (PP4 < 0.75, n = 20), with no MR instrument (n = 3) or with high heterogeneity of the instrument (n = 2), ten genes remained (PALMD, NAV1, PRRX1, ATP13A3, BCL10, TWIST1, RAD9A, NRBP1, FES and AFAP1). All genes showed a significant association in MR with the inverse-variance weighted and weighted median approaches (P_IVW and P_WM adjusted for false discovery rate <5%), suggesting a causal association between gene expression in the aortic valve and CAVS. There was no evidence of pleiotropy with the Egger intercept test (P_intercept > 0.05 for all genes). Among the novel genes identified, three were located at a genome-wide significant locus: PRRX1, ATP13A3 and TWIST1, for which gene expression in the aortic valve was positively associated with CAVS (P_IVW = 4.5 × 10^−17, 9.4 × 10^−9, and 5.7 × 10^−4, respectively) (Fig. 2). The lead SNPs in the meta-analysis at these three loci were strong eQTLs in the aortic valve (all P_eQTL < 1 × 10^−15). Five other genes were identified; the association between genetically predicted expression and CAVS was positive for three: BCL10, RAD9A and NRBP1 (P_IVW = 0.0066, 5.6 × 10^−5, and 9.4 × 10^−4, respectively) and negative for two: FES and AFAP1 (P_IVW = 7.9 × 10^−4 and 5.5 × 10^−4, respectively) (Supplementary Figs. 6 and 7).

a, d, g LocusCompare plots at the PRRX1, ATP13A3, and TWIST1 loci. P for calcific aortic valve stenosis was obtained from the inverse-variance weighted fixed-effect GWAS meta-analysis. P for valve eQTL was obtained from the nominal association between genotype and normalized gene expression. b, e, h Boxplots showing normalized gene expression in human aortic valves according to the genotype at the lead SNP at the PRRX1, ATP13A3, and TWIST1 loci. The center mark in the box represents the median, the bounds of the box represent the 25^th and 75^th percentiles and the whiskers are the most extreme data point, which is no more than 1.5 times the interquartile range. P GWAS was obtained from the inverse-variance weighted fixed-effect GWAS meta-analysis. P eQTL was obtained from the nominal association between genotype and normalized gene expression. The allele in red is the risk allele. c, f, i Scatterplot representing the effect of each SNP selected in the instrument for the Mendelian randomization (MR) analysis on gene expression in human aortic valves (n = 484) and risk of calcific aortic valve stenosis (n = 14,819 cases and 927,044 controls) at the PRRX1, ATP13A3, and TWIST1 loci. Data are presented as the effect and 95% confidence interval (+/−1.96*standard error). Red line: inverse-variant weighted (IVW) MR; Dotted red lines: 95% confidence interval for IVW MR; Green line: Weighted median MR; Pink line: Egger MR.

To evaluate a potential impact of valve morphology on gene expression regulation for these genes, we calculated eQTL separately for tricuspid and bicuspid aortic valves. For all 10 genes identified by the transcriptomic analyses, the effect of the lead variant on gene expression was similar in tricuspid and bicuspid valves (all P < 0.05 for heterogeneity) (Supplementary Data 16).

Differential expression according to genotype at the TWIST1 locus

Considering the potential causal relationship between genetically determined expression of TWIST1 in the aortic valve and CAVS risk, as well as the fact that the lead GWAS SNP, rs7804522, is only associated with gene expression in the aortic valve (valve-specific eQTL), we explored further the impact of this variant on the transcriptomic profile in the aortic valve. After selecting individuals with severe aortic stenosis, we compared gene expression between the 66 individuals homozygous for the risk allele (rs7804522-C) and the 138 individuals homozygous for the other allele. We found 509 differentially expressed genes at a false discovery rate <5%, indicating a profound effect of this genotype on transcriptome-wide gene expression (Fig. 3 and Supplementary Data 17). For comparison, a differential expression analysis performed using the same method for the four other genome-wide significant loci corroborated by transcriptomic analyses in the aortic valve (PALMD, NAV1, PRRX1 and ATP13A3) showed only five or less differentially expressed genes between individuals homozygous for the risk allele and those homozygous for the other allele (Supplementary Data 18). At the TWIST1 locus, a total of 148 genes were up-regulated while 361 genes were down-regulated in the group homozygous for the risk allele (Fig. 3). The dysregulated genes were enriched for several metabolic pathways including neutrophil degranulation, cell cycle, cell division and apoptosis (Fig. 3 and Supplementary Data 19).

a Volcano plot representing the differentially expressed genes between individuals homozygous for rs7804522-C and individuals homozygous for rs7804522-G. Red points represent up-regulated genes (n = 148). Blue points represent down-regulated genes (n = 361). Statistical significance was set at P ≤ 0.001 corresponding to false-discovery rate <5% in the differential expression analysis. The top 10 up and down-regulated genes are labeled. b Lead independent enriched terms for the dysregulated genes. The statistical significance of the association for each term was obtained from a hypergeometric test and is illustrated by the color (p-value). The number of overlapping genes is illustrated by the size of the bubble.

Prioritization of causal genes

To prioritize potential causal genes, we combined the evidence from ten different features, including lead variants location and annotation, association at the gene level, expression in the aortic valve, colocalization, TWAS and MR. Taken together, 17 genes located at genome-wide significant loci and eight genes identified by TWAS had four or more features that suggested their implication in CAVS (Fig. 4). Among the genes located at genome-wide significant loci, PRRX1, NAV1, ALDH1A2 and MUC4 had seven or more features while ACTR2, SMAD9, PALMD, ATP13A3, FADS2 and TWIST1 had four or more features including significant associations in the TWAS and MR analyses. Out of the 25 prioritized genes, seven encoded druggable human proteins (MUC4, ALDH1A2, NRBP1, FES, PDGFRA, LPL and NPC1). Existing drugs interacting with these genes according to the drug-gene interaction database are reported in Supplementary Data 20.

Nearest gene: gene closest to a lead SNP in the GWAS meta-analysis; Intronic: annotation of the lead SNP in the meta-analysis; Missense or nonsense: lead GWAS SNP is in linkage disequilibrium (r²â€‰â‰¥â€‰0.8) with a missense or nonsense variant for the gene; MAGMA: significant in MAGMA analysis at P < 0.00039 corresponding to false discovery rate <5%; High expression in valve: above the 90th percentile of all protein-coding genes; Valve-specific expression: expression specificity score >0.1; COLOC: colocalization PP4 > 0.75; Valve eQTL: significant eQTL; TWAS: significant in transcriptome-wide association study at P < 0.00017 corresponding to false discovery rate <5%; MR: significant in Mendelian randomization analyses (P < 0.05). Gold squares indicate a significant positive association; Blue squares indicate a significant negative association.

Pathway enrichment

A pathway analysis was performed with the Metascape tool by including the genes identified using MAGMA as well as the 35 genes significant in the TWAS analysis. Among the top significantly enriched terms, we identified regulation of interleukin-6 production, embryonic development, regulation of osteoblast differentiation, response to growth factor and plasma lipoprotein assembly, remodeling, and clearance (Supplementary Fig. 8 and Supplementary Data 21).

Cross-phenotype analyses

We used the interactive cross-phenotype analysis of GWAS (iCPAG) database to identify other phenotypes that share genetic associations. After adjustment for multiple testing, there were 35 significantly associated phenotypes (P_Bonferroni < 0.05), including plasma lipids (lipoprotein (a), LDL-cholesterol, triglycerides, apolipoprotein B, total cholesterol, lipoprotein-associated phospholipase A2), blood pressure traits (pulse, diastolic and systolic blood pressure), other cardiovascular diseases (coronary artery disease, peripheral arterial disease, carotid atherosclerosis, metabolic syndrome), but also aortic root size, coronary artery calcification, bone mineral density, resting heart rate and leukocyte count (Supplementary Data 22).

The association of the lead GWAS and TWAS variants with 44 other cardiovascular traits was then evaluated in UK Biobank (Fig. 5). The CAVS risk allele for several variants was positively associated with circulating lipids such as LDL-cholesterol and apolipoprotein B or blood pressure traits. Notably, 14 variants had a significant association with pulse pressure while 11 were associated with coronary artery disease with a direction of effect concordant with CAVS.

The nearest gene is indicated on the x-axis; * denotes TWAS loci. The traits are from the following categories: blood lipids (orange), blood pressure (brown), anthropometric traits (gray), cardiovascular diseases (red), diabetes (blue) and others (green). The effect size and statistical significance of the genetic association with the phenotype in UK Biobank are illustrated by the color (Z-score) and the size of the bubble (p-value).

Accordingly, there was a significant positive genetic correlation between CAVS and lipids, blood pressure traits and coronary artery disease. A positive genetic correlation was also observed for abdominal aorta calcification, ischemic stroke, peripheral artery disease, body-mass index, diabetes and C-reactive protein (Fig. 6 and Supplementary Data 23).

Data are presented as genetic correlation (r_g) and 95% confidence interval (+/−1.96*standard error). The traits are from the following categories: blood lipids (orange), blood pressure (brown), anthropometric traits (gray), cardiovascular diseases (red), diabetes (blue) and others (green). LDL low-density lipoprotein, HDL high-density lipoprotein, BMI body mass index, WHR waist-to-hip ratio, WHRadjBMI waist-to-hip ratio adjusted for body mass index.