A transcriptome-wide association study identifies PALMD as a susceptibility gene for calcific aortic valve stenosis

Calcific aortic valve stenosis (CAVS) is a common and life-threatening heart disease and the current treatment options cannot stop or delay its progression. A GWAS on 1009 cases and 1017 ethnically matched controls was combined with a large-scale eQTL mapping study of human aortic valve tissues (n = 233) to identify susceptibility genes for CAVS. Replication was performed in the UK Biobank, including 1391 cases and 352,195 controls. A transcriptome-wide association study (TWAS) reveals PALMD (palmdelphin) as significantly associated with CAVS. The CAVS risk alleles and increasing disease severity are both associated with decreased mRNA expression levels of PALMD in valve tissues. The top variant identified shows a similar effect and strong association with CAVS (P = 1.53 × 10−10) in UK Biobank. The identification of PALMD as a susceptibility gene for CAVS provides insights into the genetic nature of this disease, opens avenues to investigate its etiology and to develop much-needed therapeutic options.


Main text
Calcific aortic valve stenosis (CAVS) is the most common valvular heart disease (2% >65 years old) 5 . It is characterized by a progressive remodeling and calcification of the aortic valve, leading to stenosis and heart failure. There is a long latent period before CAVS becomes severe and symptomatic, which provides a window of opportunity for intervention 6 . Unfortunately, conventional cardiovascular drugs are unable to stop or delay the progression of CAVS 7 . If severe CAVS is left untreated, the survival of affected patients is dramatically shortened following the onset of symptoms 6 . The only treatments available for symptomatic patients with severe CAVS are invasive and costly surgical or transcatheter interventions 8 . Finding new molecular targets to halt disease progression is an urgent priority 1 .
Our molecular and genetic understanding of CAVS is currently limited. A strong genetic component is suggested by early-onset forms of the disease, familial aggregation and genetic epidemiology studies 9,10 . A number of susceptibility genes were identified 3,11-14 . Genome-wide gene expression of normal and stenotic human valves have highlighted biological processes that are involved in CAVS 15,16 . In this study, we combined GWAS and valve eQTL results to identify the molecular drivers of CAVS. Table 1 shows the clinical characteristics of 1,009 cases and 1,017 controls that passed quality controls (QC) in the QUEBEC-CAVS cohort. Mean age was 71.7±8.3 years and 64% were men.
About half (53%) of CAVS cases and the majority of controls (98%) had coronary artery disease.
An expression quantitative trait loci (eQTL) mapping study was then conducted in the most disease-relevant tissue to study CAVS, i.e. human aortic valve. Valve eQTL were calculated in 233 patients that passed QC for genotyping and gene expression. A total of 10,598 independent valve eQTL were identified at P eQTL <1x10 -8 (Supplemental Table 2). Overall, these eQTL involved 2,277 genes/probes significantly associated with one to 41 independent SNPs (r 2 <0.8) (Supplemental Fig. 3a). On average, a single SNP explained 26.5% of the gene/probe expression variance. However, 29.1% of the eQTL-SNPs explained more than 30% of the expression variance (Supplemental Fig. 3b).
A transcriptome-wide association study (TWAS) 4 was then performed using the GWAS and valve eQTL datasets. A Manhattan plot showing transcriptome-wide association in valve tissue with CAVS is shown in Fig. 1b. Only one probe-CAVS association corresponding to the PALMD gene reached genome-wide significance (P TWAS =0.00007). PALMD is located on chromosome 1p21.2 and the top GWAS SNPs were the same as the top valve eQTL-SNPs (Fig.   2a). rs6702619, a genotyped variant, was the SNP in this locus most significantly associated with both CAVS (OR=1.29, 95% CI 1.14-1.46, P GWAS =6.12x10 -5 ) and the expression of PALMD (P eQTL =5.82x10 -33 ). rs6702619 is a common SNP (minor allele frequency of 47.7% in Europeans from the 1000 Genomes Project) located 65.2 kb from the transcriptional start site of PALMD.
The valve eQTL rs6702619-PALMD indicated that the risk allele for CAVS ("G") is associated with lower mRNA expression levels of PALMD in valve tissues (Fig. 2b), suggesting that lower expression increases the risk of CAVS. Variants located within 1 Mb of PALMD that increased its expression tended to decrease the risk of CAVS (Fig. 2c). A Mendelian randomization analysis suggested a causal role of lower PALMD expression on CAVS risk, without evidence of pleiotropy (P=0.00046; Egger intercept P=0.19; Supplemental Fig. 4). Results remained significant after excluding rs6702619 from the analysis (P<0.05). Together, these results indicated that risk variants on 1p21.2 confer susceptibility to CAVS through down-regulation of PALMD in aortic valve tissues. The relationship between PALMD expression levels and CAVS severity was also consistent with this direction of effect. Lower normalized, age and sex adjusted PALMD expression was associated with smaller aortic valve area (P=0.0027), higher mean transvalvular gradient (P=0.0001), and higher peak transvalvular gradient (P=8.13x10 -5 ) (Supplemental Fig. 5).
We additionally looked for PALMD eQTL in the Genotype-Tissue Expression (GTEx) dataset 17 . rs6702619-PALMD association was not significant in 44 tissues from GTEx, suggesting that the change in expression is specific to aortic valve tissue. Evaluating other SNPs within or near PALMD revealed significant eQTL with this gene in the pancreas, tibial nerve and subcutaneous adipose tissue. However, these eQTL-SNPs were not in LD with CAVS-associated variants.
We then replicated our newly associated locus using data recently released by the UK Biobank, which included 1,391 CAVS cases and 352,195 unaffected individuals all of European ancestry.
We observed a very strong association with an almost identical effect on CAVS risk for our top SNP rs6702619 with an OR of 1.27 (95% CI 1.18-1.37) (Fig. 3a). The association was also replicated for rs7543130, a genotyped variant in perfect LD with rs6702619. At the genome-wide scale level, only two loci were genome-wide significant in UK Biobank (Fig. 3b): LPA on 6q25.3-q26 previously associated with aortic valve calcification 3 and PALMD on 1p21.2 identified in this study.
The top GWAS SNP in this study (rs6702619) was recently identified by the EchoGen consortium to be associated with aortic root diameter in a large GWAS meta-analysis including up to 46,533 individuals 18 . The EchoGen consortium also identified a SNP in CACNA1C (calcium channel, voltage-dependent, L type, alpha 1C subunit) associated with aortic root diameter 18 , which we recently identified as a susceptibility gene of CAVS 14 . As indicated in the later study, rs6702619 is located within enhancer histone marks and DNase-hypersensitive sites in Encyclopedia of DNA Elements (ENCODE) data 19 . However, rs6702619 was not found to act as cis-eQTL in whole blood, monocytes, and myocardial tissue 18 . In this study, we found PALMD consistently expressed in human aortic valve tissues and rs6702619 strongly associated with mRNA expression levels of PALMD, providing a potential functional mechanism of how genetic variants on 1p21.2 may increase CAVS susceptibility. Furthermore, a Mendelian randomization analysis suggested that PALMD expression is causally associated with the risk of CAVS.
The frequency of rs6702619 is high in our French Canadian population (MAF=47%) and the modest effect size (OR=1.29) is consistent with those observed for other complex traits. The estimated population-attributable risk indicates that more than 12.5% of cases in our population and 12.3% of cases in UK Biobank could be attributed to this common variant. Of note, the frequency of rs6702619 varies greatly between ethnic groups (Supplemental Fig. 6). In the 1000 Genomes Project, allele "G" had a frequency of 48% in Europeans, 8% in Africans, 7% in East Asians and 25% in South Asians. This may explain some of the discrepancy in risk of CAVS among ethnic groups 20,21 .
PALMD is a distant homolog of a small paralemmin protein family 22 . PALMD is expressed in many tissues, but most abundantly in cardiac and skeletal muscle 23 . It is a mainly cytosolic protein, localized predominantly in actin filaments, and may be implicated in plasma membrane dynamics and cell shape control as other members of this family 22,24 . However, its molecular and cellular functions remain largely unknown. PALMD was recently shown to promote myoblast differentiation and muscle regeneration 25 . Considering the common embryologic origin of the aortic valve and the aortic root, which both arise from the secondary heart field, a potential role on smooth muscle and valve interstitial cell differentiation could explain the association with the two phenotypes 26,27 . PALMD was also identified as a pro-apoptotic gene induced by p53 in response to DNA damage in osteosarcoma cell lines 28 . However, its role in cell death was not confirmed in other cell types 25 .
Of note, the identified SNP, rs6702619, was not associated with aortic valve calcification measured by computed tomographic scanning in a large GWAS including 6,942 individuals of European ancestry (P=0.557) 3 . This suggests that calcification is not the main mechanism by which PALMD expression modulates CAVS progression and could rather appear only in later stages of the disease. The variant identified is also expected to act via pathways not involved in the risk of coronary artery disease (CAD) since the control individuals included in our study almost all had CAD. rs6702619 was indeed not associated with CAD in a recent large GWAS meta-analysis (P=0.96) 29 .
Our study has some limitations. First, only individuals of European ancestry were included in the analysis, therefore results cannot be generalized to other ethnic groups. Second, the power in our GWAS analysis was limited by the number of available samples. However, expression data on the most disease-relevant tissue and replication in a large population-based cohort confirmed the robustness of the findings. Third, aortic valve tissue was only available in stenotic valves; patterns of expression of PALMD in normal valves might differ.
In conclusion, using a TWAS approach with expression data in aortic valve tissue, we identified a new CAVS susceptibility gene, PALMD, with replication in a large prospective populationbased cohort. Further analyses of gene expression in valve tissue and disease severity suggested that the identified variant acts by decreasing the expression of PALMD, a finding supported by Mendelian randomization. Further studies are warranted to elucidate the exact mechanism of action and evaluate the potential for targeted therapeutic interventions. patients was also matched in a 1:1 ratio with cases for age, gender, type 2 diabetes and hypertension. Patients with a history of severe valvular heart disease (at any of the 4 valves), with significant aortic valve regurgitation (grade > 2) or with end stage renal disease (eGFR<15 mL/min/1.73m 2 ) were excluded. All patients signed an informed consent for the realization of genetic studies. Demographics, anthropometric measurements, lifestyle factors, previous and current medical history, current medication and blood pressure measurements were collected. In addition, plasma lipids and creatinine were measured. CAVS cases underwent a comprehensive Doppler echocardiographic examination. The transvalvular gradient was calculated using the modified Bernoulli equation and the aortic valve area calculated with the continuity equation.

References
For this study, 1,033 cases and 1,037 controls were available.

Genome-wide association study
Blood samples were collected and DNA was extracted from frozen buffy coat. respectively. Genotypes were then imputed with the Michigan Imputation Server 32 using the Haplotype Reference Consortium version 1 (HRC.r1-1) data 33 as reference set (2016-03-31).
Variants with an r 2 value of ≤0.3 or MAF <1% were removed from further analysis. A total of 7,732,680 imputed SNPs passed quality control. Genetic association tests were performed using additive logistic regression models based on expected genotype counts (dosages) as implemented in the software SNPTEST v2.5.2 34 , adjusting for age, sex, and the first ten ancestry-based principal components. The genomic inflation factor for the main case-control analysis was 1.03.
The genome-wide significant P value cutoff was set to 5x10 -8 . Regional plots were created with LocusZoom 35 . The population-attributable risk was calculated as PAR% = 100% × P × (OR -1)/[P × (OR -1) + 1], where P is the frequency of the risk allele associated with CAVS in the control group, and OR is the odds ratio calculated in the case-control cohort.

Gene expression analysis in human aortic valves
Transcriptomic analyses were performed from 240 stenotic aortic valves collected as mentioned above. All stenotic valves were tricuspid and had a fibro-calcific remodeling score of 3 or 4 36 .
Selected valves were from cases that are part of the GWAS and included 120 men and 120 women. RNA was extracted from valve leaflets, and gene expression evaluated using the Illumina HumanHT-12 v4 Expression BeadChip. Standard microarray processing and quality control analysis was performed 37 . The raw data were quantile normalized after log2transformation with the lumi package in R 38 . Only one sample failed quality control, leaving 239 samples for subsequent analyses. Probe sequences were mapped to RefSeq B38, GENCODE v24 B38, mRNA B38, and the human genome (GRCh B38) using Bowtie, and probes not mapping to any coding region were removed. A total of 45,699 probes remained after this step. Expression data were then adjusted for age and sex using residuals obtained with the robust fitting of linear models function (rlm) in the R statistical package MASS. Residual values deviating from the median by more than 3 SD were filtered out.

Expression quantitative trait loci
Only subjects that passed genotyping and gene expression quality controls were considered for eQTL analysis leaving a sample size of 233. eQTL were identified by using linear regression model and additive genotype effects as implemented in the Matrix eQTL package in R 39 . Cis-eQTL were defined by a 2 Mb window, i.e. 1 Mb distance on either side of the SNP. eQTL were calculated on adjusted expression traits to obtain test statistics, P values and false discovery rate.
Estimates of effect sizes were obtained with PLINK.

Transcriptome-wide association study (TWAS)
The TWAS was performed using FUSION 4 . Briefly, the valve eQTL was the reference data including 233 individuals that passed QC for both gene expression and genotyping. This reference set was used to calculate gene expression weights, which were computed one probe at a time using SNP genotyping data located 500 kb on both sides of the probe using prediction models implemented in FUSION including 1) the single most significant valve eQTL-SNP as the predictor (top1), 2) LASSO regression, and 3) elastic net regression (enet). All probes that passed QC in the valve eQTL were evaluated (n=45,699). Expression weights were then combined with summary-level GWAS results to estimate association statistics between gene expression and CAVS. Genome-wide significant TWAS genes were considered at P TWAS <0.0001.

Mendelian Randomization
We first selected variants located within 1 Mb of the gene of interest identified in the TWAS, PALMD, and significantly associated with PALMD gene expression at a threshold of P<0.05.
We kept only independent variants at a linkage disequilibrium threshold of r 2 < 0.05. Mendelian randomization was performed using the Wald method by regressing genetic effect estimates for CAVS risk as determined in the GWAS analysis (dependent variable) on genetic effect estimates for PALMD gene expression as determined in the eQTL analysis. Effect estimates were adjusted for the minor allele frequency of each variant. To determine the significance of the association, a bootstrap method was used. Predicted effects on CAVS risk and PALMD gene expression were sampled from a normal distribution with mean and standard deviation as determined from the GWAS and eQTL analyses. A two-tailed P value was calculated using 100,000 random simulations. To determine the presence of unmeasured net pleiotropy, we performed Egger Mendelian randomization in which a non-zero y-intercept is allowed in order to assess violations of standard Mendelian randomization 40 .

Gene expression according to severity
CAVS disease severity was assessed by aortic valve area, mean and peak transvalvular gradients.
The influence of normalized, age and sex adjusted PALMD expression levels on disease severity was tested in 239 cases using linear regression models. 95% confidence intervals were estimated using the predict function in R. We performed additive logistic regression analysis using SNPTEST v2.5.2, adjusting for age, sex, and the first ten ancestry-based principal components to evaluate the effect of the top SNP identified in our TWAS analysis. We then performed a fixed-effect meta-analysis using the inverse-variance weighted method. To look for other significant variants associated with CAVS in UK Biobank, we performed a genome-wide association study from the genotyped data using additive logistic regression models in PLINK adjusting for age, sex and the first ten principal components. The genomic inflation factor for the main case-control analysis was 1.01. The genome-wide significant P value cutoff was set to 5x10 -8 .

Statistical analysis
Statistical analyses were performed with R version 3.2.3 unless otherwise specified. 2-sided P values below 0.05 were considered significant unless otherwise specified.