Genetic loci associated with heart rate variability and their effects on cardiac disease risk

Reduced cardiac vagal control reflected in low heart rate variability (HRV) is associated with greater risks for cardiac morbidity and mortality. In two-stage meta-analyses of genome-wide association studies for three HRV traits in up to 53,174 individuals of European ancestry, we detect 17 genome-wide significant SNPs in eight loci. HRV SNPs tag non-synonymous SNPs (in NDUFA11 and KIAA1755), expression quantitative trait loci (eQTLs) (influencing GNG11, RGS6 and NEO1), or are located in genes preferentially expressed in the sinoatrial node (GNG11, RGS6 and HCN4). Genetic risk scores account for 0.9 to 2.6% of the HRV variance. Significant genetic correlation is found for HRV with heart rate (−0.74<rg<−0.55) and blood pressure (−0.35<rg<−0.20). These findings provide clinically relevant biological insight into heritable variation in vagal heart rhythm regulation, with a key role for genetic variants (GNG11, RGS6) that influence G-protein heterotrimer action in GIRK-channel induced pacemaker membrane hyperpolarization.

NOTE: Expected -log 10 (p-values) assuming a normal distribution are shown on the x-axis. The y-axis depicts the observed -log 10 (pvalues). Only SNPs with a minor allele frequency >1% and that were present in at least 1/3 of the sample are shown. NOTE: Significance thresholds to include HRV SNPs in the polygenic risk scores are shown on the x-axis and were p < 5.0x10 -8 ; p < 5.0x10 -7 ; p < 5.0x10 -6 ; p < 5.0x10 -5 ; p < 5.0x10 -4 ; p < 5.0x10 -3 ; p < 0.05; p < 0.5; and p < 1. The y-axis depicts the percentage of variance explained by the polygenic risk score. Figure 6: Diagnostic plot for the association between genetic risk alleles for the HRV traits and heart rate. NOTE: Each SNP is plotted by its decreasing effect on HRV per risk effect allele (x-axis) versus the estimated effect of that allele on the risk for high heart rate (y-axis). A solid red line shows the effect size estimate, called α, for the risk score on HRV. The 95% CI of α is represented by red dashed lines. The grey vertical lines indicate the 95% confidence interval (CI) of the effect on heart rate for each individual SNP. The estimated effects on heart rate are in beats per minute (bpm). Graphs from top to bottom are for SDNN, RMSSD, and pvRSA/HF, respectively. and (e-f) pvRSA/HF. NOTE: In panels (a), (c), and (e) graphical representations of DEPICT's tissue enrichment analysis are provided, highlighting annotated tissues that are enriched for expression of genes located within ±40kb of SNPs associated with p < 10 -5 to SDNN, RMSSD, and pvRSA/HF, respectively (enrichment in orange p < 0.05). In panels (b), (d), and (f) graphical representations of DEPICT's reconstituted gene set enrichment analysis are shown (p < 0.05 after Bonferroni correction) for examining SDNN, RMSSD, and pvRSA/HF, respectively. DEPICT calculates every gene's likelihood to be a member of KEGG, GEO, or REACTOME-based gene sets amongst others (N=14,461) resulting in reconstituted gene sets. Shown are the reconstituted gene sets enriched for genes in SDNN-, RMSSD-, and pvRSA/HF-associated loci (nodes), respectively, and the significant interactions between these gene sets (edges).      We performed three corrections: (1) analytical construction of a meta-analysis of the log-transformed coefficients of variation (CV) of SDNN -ln ((SDNN/mean heart period)*100%) and RMSSD-ln ((RMSSD/mean heart period)*100%) using GWIS on the HRV and heart rate GWASs summary statistics; (2) meta-analysis of SNP associations with the log-transformed CV(SDNN) and CV(RMSSD) in Lifelines, NESDA, and TRAILS-Pop; (3) mediation analysis using Sobel's test testing whether heart rate mediates the SNP association with HRV in a meta-analysis of HRV traits and heart rate in LifeLines, NESDA, and TRAILS-Pop. For the mediation analysis for pvRSA/HF only data of NESDA and TRAILS-Pop were available. Significant p-values after Bonferroni correction for 11 independent SNPs (<0.0045) are shown in bold. Chr: chromosome; HR: heart rate; CV: coefficient of variation; SE: standard error. # Lifelines (N=12,101), NESDA (N=2,118), TRAILS-Pop (N=1,191). *Sobel p-values of the meta-analyses are calculated from meta z-scores that are determined using a sqrt(sample size) weighted z-score meta-analysis. Table 9: Association between HRV SNPs and heart rate. Results in (a) are the look-up results of the 17 HRV SNPs identified in this study in the metaanalysis of heart rate (Den Hoed et al., 2013). Panels (b) and (c) present the explained variances in heart rate in the Lifelines (n=12,101), NESDA (n=2,118), TRAILS-Pop (n=1,191), and ABCD (n=1,094) cohorts by (b) the weighted multi-SNP genetic risk score based on the independent genome-wide significant HRV SNPs in the stage 1+2 meta-analysis, and (c) the optimal polygenic risk scores computed at the p-value threshold that explained the largest percentage of phenotypic variance. ΔR 2 is the difference in percentage of explained variance by the multi-SNP genetic or polygenic risk score between the models with and without the risk score while adjusting both for age, sex, and principal components. The maximum number of overlap between the HRV and heart rate meta-analysis was N=10,612. * Weighted polygenic risk score were determined based on independent SNPs in the stage 1 meta-analysis. For NESDA and TRAILS-Pop the weights (i.e. effects sizes) and p-values were adjusted by analytically extracting the cohort's effect size and standard error from the meta effect size and standard error, respectively, and recalculating the p-value based on these adjusted effect size and standard error, since these cohorts were included in stage 1. * For Lifelines, NESDA, and TRAILS-Pop the weights (i.e. effects sizes) and number of genome-wide significant SNPs included in the risk score were adjusted by analytically extracting the cohort's effect size and standard error from the meta effect size and standard error, respectively, and recalculating the p-value based on these adjusted effect sizes and standard errors, since these cohorts were included in stage 1 and/or 2.

Supplementary
* Weighted polygenic risk score were determined based on independent SNPs in the stage 1 meta-analysis. For NESDA and TRAILS-Pop the weights (i.e. effects sizes) and p-values were adjusted by analytically extracting the cohort's effect size and standard error from the meta effect size and standard error, respectively, and recalculating the p-value based on these adjusted effect size and standard error, since these cohorts were included in stage 1.  * p-value from z-score weighted meta-analysis of all cohorts using METAL and β and standard error from inverse-variance metaanalysis of only HF cohorts using GWAMA were used # Effect size is either the incremental change in the phenotype for quantitative traits or the odds ratio for binary traits when the multi-SNP genetic risk score of HRV decreases by one unit. * p-value from z-score weighted meta-analysis of all cohorts using METAL and β and standard error from inverse-variance metaanalysis of only HF cohorts using GWAMA were used. Chr: chromosome; bp: base pair position (build 36). Chr: chromosome; bp: base pair position (build 37); LD: linkage disequilibrium (r 2 ). TSS200= methylation site is located within 200 bp of a transcription start site; TSS1500= methylation site is located within 1500 bp of a transcription start site; 3'UTR= methylation site is located in the 3'untranslated region; Body= methylation site is located within a gene (gene body); Island= methylation site is located in a CpG island; S_Shelf= methylation site is located in the south shelf of a CpG island; N_Shelf= methylation site is located in the north shelf of a CpG island; S_Shore= methylation site is located in the south shore of a CpG island; N_Shore= methylation site is located in the north shore of a CpG island * Caution should be exercised when using the p-value as a measure of confidence that the highlighted genes are indeed causal for the associations identified by GWAS, or as a means to compare evidence across tools. Each tool uses different databases and algorithms to generate these p-values and as such, p-values cannot necessarily be compared across tools. We advise to use the presence of a p-value as an indication that the respective tool identifies a gene as potentially being causal.

Supplementary
# DEPICT is considered to be the most advanced tool since it uses the most up-to-date information to derive which genes are likely to be causal. For DEPICT, only identified genes that were located within ±40 kb of GWAS lead SNPs are reported.

Gene-based GWAS (VEGAS)
The results of the VEGAS gene-based GWAS analysis corroborate those of the SNP-based analysis (Supplementary Table 17). In seven of our eight loci with genome-wide significant SNPs (Table 1) the closest gene or other nearby genes also emerged from VEGAS as significantly associated with at least one of the HRV traits. Only the loci on chromosome 6 (PPIL) and 15 (NEO1) were not represented. In addition, a few new genes close to the top hits in our SNP-based analyses emerged: TFPI2 and GNGT1 (both on chromosome 7), ALG10, ALG10B, and CPNE8 (all on chromosome 12), and NRTN, FUT6, FUT3, FUT5, VMAC, CAPS, RANBP3, and RFX2 (all on chromosome 19). Only one gene in a new locus was identified: CCDC141 (on chromosome 2, next to TTN coding for the titin protein, a key determinant of myocardial passive stiffness). CCDC141 and CPNE8 have previously been reported to be associated in a GWAS for heart rate 14 .

Heritability and genetic correlations
The common SNP heritability estimated by Genomic Restricted Maximum Likelihood (GREML) analysis was 10.8% and 13.2% for SDNN and RMSSD in the Lifelines cohort (Supplementary  Table 18, panel a). Applying LD score regression analysis on the summary statistics of the stage 1 meta-analysis suggested a common SNP heritability of 11.1% for SDNN, 11.2% for RMSSD, and 11.8% for pvRSA/HF (Supplementary Table 18, panel b). Classical modeling on family data from the Oman Family Study yielded heritability estimates between 15.2 and 20.0%. Previous studies showed high phenotypic correlation between pvRSA/HF and RMSSD 31, 32 and suggested a large overlap in their genetic architecture 33,34 . We confirmed the phenotypic correlations in the TRAILS-Pop, NESDA, and Lifelines cohorts where they ranged from 0.70 (SDNN-pvRSA/HF) to 0.96 (RMSSD-pvRSA/HF). Bivariate analysis using GREML and LD score regression analysis further showed genetic correlations of unity between pvRSA/HF and RMSSD, but also with SDNN, indicating complete overlap between the genetic variants influencing all three HRV traits (Supplementary Table 18, panel a+b). Modeling of the family data from the Oman Family Study confirmed high genetic correlation between the HRV traits, including SDNN (0.71-0.90) (Supplementary Table 18, panel c). This shows that SDNN, notwithstanding its different etiology, is influenced by sets of genetic variants that also influence the RMSSD and pvRSA/HF traits.
We observed large negative genetic correlations (between -0.74 and -0.55) of SDNN, RMSSD, and pvRSA/HF with heart rate (Supplementary Table 11, panel b). High genetic overlap between heart rate and HRV fits the notion that HRV increases with a net increase of the tonic vagal effects on the sinoatrial node, which in turn reduces average heart rate. However, we note that the full relationship between HRV and heart rate is more complex because heart rate is under parallel sympathetic control, changes in which are not captured well by the HRV measures used 35,36 and because changes in heart rate can lead to lower HRV, even at unchanged vagal input 37 .

Association between heart rate SNPs and HRV
After establishing the effects of HRV SNPs on the variance in heart rate we also examined the reverse question, i.e. whether SNPs known to be associated with heart rate accounted for part of the variance in HRV. The reverse association of the 21 heart rate SNPs identified by a GWAS study for heart rate 14 with the three HRV traits showed that nine heart rate SNPs were associated with HRV, all in the expected direction, after correction for multiple testing (Supplementary  Table 10, panel a)). Using gtx, a multi-SNP genetic risk scores for heart rate was associated with the three HRV traits (Supplementary Table 10, panel b).
In the Lifelines, NESDA, TRAILS-Pop, and ABCD cohorts, a multi-SNP genetic risk score for heart rate explained 0.2 to 0.9% of HRV variance in the four cohorts (Supplementary  Table 10, panel c). The full polygenic risk score for heart rate explained a similar amount of the variance in the three HRV traits: 0.2 to 0.8% (Supplementary Table 10, panel d).

Network and functional enrichment analyses
Based on our 17 HRV SNPs, we selected candidate genes as input for the interaction network analysis with the GeneMANIA algorithm to prioritize potentially causal genes 38 with additional input from the VEGAS analyses (Supplementary Table 17) and four publicly available bioinformatics tools (Supplementary Table 19). As one of the genes (LINC00477) could not be found by GeneMANIA, 24 query genes were used as input. Functional network and enrichment analysis on these genes resulted in 31 significantly enriched gene ontology terms of which 23 had false discovery rate ≤ 0.01 (Supplementary Table 20). These 23 gene ontology terms showed that the genes near our hits were broadly related to two categories of significantly enriched biological processes, namely: (1) cellular signaling and cellular responses, including G-protein coupled receptor protein signaling and the responses to glucagon and peptides, and (2) metabolic processes in the cell (e.g. fucosylation and glycosylation) (Supplementary Fig. 10).

Tissue and gene-set enrichment analyses
In silico tissue enrichment analysis using DEPICT 6 highlights a role for hormones in HRV regulation, with enrichment for the adrenal cortex, endocrine glands, gonads, gastrointestinal tract and female reproductive organs ( Supplementary Fig. 8a,c,e). Gene-set enrichment analysis using DEPICT highlights the importance of cardiac development ( Supplementary Fig. 8b,d,f).

Known biological functions of the genes closest to the HRV SNPs
Here we describe the known biological functions of the genes closest to the eight HRV lead SNPs identified in the two-stage meta-analysis. NDUFA11 (Chr 19, full name: NADH dehydrogenase (ubiquinone) 1 alpha sub-complex, 11, also known as B14.7 and Cl-B14.7) encodes a subunit of the membrane-bound mitochondrial complex I involved in mitochondrial respiration and electron transport. NDUFA11's cDNA was cloned and sequenced from human heart mitochondria. The non-synonymous lead SNP in our meta-analysis induces a change in the structure of this complex with as yet unknown functional effects. Clinically, five phenotypes are frequently seen with homozygous mutations of this gene: severe neonatal lactic acidosis, cardiomyopathy-encephalopathy, hepatopathytubulopathy, leukodystrophy with macrocephaly, and Leigh's and Leigh-like neurodegenerative disorder, which are the most common presentations of complex I deficiency at a young age 39,40 . LINC00477 (C12orf67) (Chr 12, full name: long intergenic non-protein coding RNSA 477) is of unknown biological function but the locus has also been associated with resting heart rate 14 and PR interval 41 . PPIL1 (Chr 6, full name: peptidylprolyl isomerase (cyclophilin)-like 1) is a member of the cyclophilin family. After being recruited by Ski interaction protein (SKIP) PPIL1 participates in the activation of spliceosomes and thereby facilitates the folding of proteins in the spliceosomes 42,43 . No relation with heart rate or HRV has been described but PPIL1 has been reported to be one of two key driver genes involved in coronary artery disease networks 44 . SYT10 (Chr 12, full name: synaptotagmin-10) has a role in the regulation of calcium dependent exocytosis, including calcium regulated release of neurotransmitter from presynaptic nerve terminals 45,46 . SYT10 has been previously identified as being associated with heart rate 14 . Synaptotagmin-1 plays a known role in calcium-triggered acetylcholine release from the neuromuscular junction 47 . Our meta-analyses results hint at a possible role for synaptotagmin-10 in sinoatrial acetylcholine release but this remains to be tested. GNG11 (Chr 7, full name: guanine nucleotide binding protein (G protein, gamma 11) which is also a previously identified heart rate gene 14 encodes the γ11 subunit of Gαβγ heterotrimers. GNG11 is one of 12 genes encoding the Gγ subunits that all undergo post-translational isoprenylation of their C termini, in case of γ11 by a farnesyl 48 . Because the γ subunits have lower amino acid sequence homology then the β subunits, they are thought to be the main determinant of signaling diversity and fidelity of the Gβγ components 49,50 . γ11 shows a unique pattern of expression compared to the other 11 γ subunits in that it is abundantly expressed in the heart but poorly in the brain 50 . RGS6 (Chr 14, full name: regulator of G-protein signaling 6) is a large gene coding one of the RGS superfamily that act as GTPase-activating proteins (GAPs) for the α subunit of Gαβγ heterotrimers by accelerating their intrinsic GTPase activity. This ends both Gα and Gβγ mediated cellular signaling. RGS6 exhibits a uniquely robust expression in the heart 10 . RGS6 is known to reduce parasympathetic signaling through the M 2 R in the sinoatrial node as well as Gα o coupled adenosine A1 receptors 10,51,52 . HCN4 (Chr 15, full name: hyperpolarization activated cyclic nucleotide gated potassium channel 4) codes for a non-GIRK potassium channel in the sinoatrial node that is a core part of the 'clock circuit' generating the pacemaker potential. HCN4 is expressed abundantly in the heart from the early embryonic phase onward 12 . HCN4 was previously found to be associated with heart rate 14 . Depending on total load and severity, loss-of-function mutations in HCN4 can lead to asymptomatic bradycardia, Brugada syndrome and Sick Sinus Syndrome, atrial fibrillation and possibly left ventricular noncompaction cardiomyopathy 53,54 . NEO1 (Chr 15, full name: neogenin 1, aka HsT16534, IGDCC2, NGN or NTN1R2) is a multifunctional transmembrane receptor closely related to the immunoglobulin (Ig) superfamily and is a netrin receptor 55 . Neogenin plays a role in early cerebellar neurite outgrowth and projection in chicken and quail 56,57 and is known to be a regulator of axonal guidance in the nervous system 58 . No role in HRV or heart rate has been described previously. KIAA1755 (Chr 20, no full name in Ensembl or at NCBI) which has also been associated with heart rate 14 . Effects of the non-synonymous lead SNP (rs6123471) is unclear as the function of this gene remains unknown.