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.

Heart rate variability (HRV) is a physiological variation in cardiac cycle duration. When measured under supine or sitting conditions, resting HRV is most prominently centred around the frequency of respiration (∼0.25 Hz) and the intrinsic blood pressure rhythm (∼0.1 Hz). This reflects modulation of tonic activity in the cardiac vagal nerves originating in cortical and subcortical nuclei1 by oscillatory input at the brainstem level from cardiorespiratory coupling, lung stretch-reflexes, and arterial chemo- and baroreceptors12. This vagal gating gives rise to oscillatory vagal effects on the pacemaker potentials in the sinoatrial node that scales with the tonic activity in the vagal nerves and provides a source of beat-to-beat variation in heart rate. Due to its good reproducibility3 and ease of measurement, HRV is a widely used non-invasive research and clinical tool to quantify the degree of vagal control of heart rate4.

Loss of cardiac vagal control as indexed by low HRV is associated with mortality in patients with cardiovascular disease5. Animal research further supports a role for cardiac vagal activity in preventing sudden death and ventricular fibrillation6. In addition, hypertension7, end-stage renal disease8 and diabetes9 are all associated with low HRV. Although the above associations may partly reflect impaired cardiac vagal control caused by these diseases, lowered HRV does not simply indicate disease severity as it also predicts all-cause mortality10 and cardiac morbidity and mortality1112 in apparently healthy individuals.

Large inter-individual differences in HRV exist in the basal resting state. Family and twin studies have uniformly confirmed a substantial genetic contribution to resting HRV with heritability estimates between 25 and 71% (ref. 13). Candidate gene studies based on current knowledge of parasympathetic nervous system biology have not yielded results that hold up in replication14. To improve our understanding of the genetic basis of HRV, we performed a two-stage meta-analysis of genome-wide association studies (GWAS) in up to 53,174 individuals of European ancestry on three HRV traits (the s.d. of the normal-to-normal inter beat intervals (SDNN), the root mean square of the successive differences of inter beat intervals (RMSSD) and the peak-valley respiratory sinus arrhythmia or high frequency power (pvRSA/HF)). These HRV traits were measured during resting, basal recordings ranging in length from ultrashort 10-s electrocardiograms (ECGs) to up to 90 min of sitting or from 2 to 12 h of daytime recording. Relevance of the identified loci for other ethnicities was examined in data from 11,234 Hispanic/Latino and 6,899 African-American individuals. In silico post-GWAS analyses were performed to test for association with cardiac disease risk factors and disease outcomes and to provide insights into the biological mechanisms by which the identified loci influence cardiac vagal control and its effect on HRV.

We detect 17 SNPs in eight loci harbouring several genes preferentially expressed in the sinoatrial node and significant negative genetic correlations of HRV with heart rate and blood pressure. These findings provide clinically relevant biological insight into heritable variation in vagal heart rhythm regulation, with a key role for genetic variants in proteins (RGS6, GNG11) known to influence G-protein heterotrimer action in GIRK-channel induced pacemaker membrane hyperpolarization.

Results

New loci associated with HRV

We meta-analysed results from GWAS on three HRV traits (see Methods section for details) performed by 20 cohorts of European ancestry in up to 28,700 individuals (Fig. 1; Supplementary Figs 1–3; Supplementary Tables 1–4). Using a significance threshold of 1 × 10−6, 23 single-nucleotide polymorphism (SNPs) in 14 loci that were associated with one or more of these HRV traits were taken forward for wet-lab genotyping or in silico replication in 11 cohorts including up to 24,474 additional individuals of European ancestry, followed by a second stage meta-analysis (Supplementary Data 1).

Figure 1

Manhattan plots of the meta-analyses of stage 1 GWAS results.

(a) SDNN, (b) RMSSD and (c) pvRSA/HF in up to 28,700 individuals of European ancestry. Only SNPs with a minor allele frequency >1% and that were present in at least 1/3 of the sample are plotted. Significant loci are shown in blue, suggestive ones in red. The blue horizontal line represents the genome-wide significance threshold. Genes closest to the lead SNPs are indicated for the loci that were genome-wide significantly associated with the trait after the stage 1+2 combined meta-analysis.

After stage 2, we identified 17 lead SNPs (11 independent) in eight loci (Table 1) that reached genome-wide significance (P<5 × 10−8). The loci on chromosomes 14 and 15 contained three and two independent signals, respectively, (Supplementary Fig. 3). Conditional analysis confirmed the presence of independently associated variants in these loci (Supplementary Table 5). In total, nine independently associated SNPs in seven loci were detected for SDNN, nine independently associated SNPs in eight loci for RMSSD, and five independently associated SNPs in five loci for pvRSA/HF. Many of the SNPs were associated with at least two of the HRV traits (Supplementary Data 1). In four loci, the lead SNPs differed between traits but were in linkage disequilibrium (LD) with each other (0.24

Table 1

Stage 1+2 combined meta-analysis results for SDNN, RMSSD and pvRSA/HF of loci that were genome-wide significant (P<(5 × 10−8)/3) in the analysis of individuals of European ancestry.

Variance explained

Weighted genetic risk scores based on the independent SNPs that reached genome-wide significance after the second stage meta-analysis were computed for the three HRV traits and used to predict RMSSD, SDNN and pvRSA/HF in adults from the Lifelines (n=12,101) and NESDA (n=2,218) cohorts, adolescents from the TRAILS-Pop cohort (n=1,191), and children from the ABCD cohort (n=1,094) (Table 2). The multi-SNP genetic risk scores were all significantly associated with HRV and the percentages of variance explained for the corresponding traits were 1.0–1.4% for SDNN, 1.1–2.4% for RMSSD, and 0.9–2.6% for pvRSA/HF. Cross-trait explained variances of genetic risk scores were close to those for the corresponding trait.

Table 2

Explained variance in HRV traits in the Lifelines (n=12,101), NESDA (n=2,118), TRAILS-Pop (n=1,191) and ABCD (n=1,094) cohorts by the weighted multi-SNP genetic risk score based on the independent genome-wide significant SNPs in the stage 1+2 meta-analysis.

To test the contribution of SNPs that did not reach genome-wide significance, we performed polygenic risk score analyses using increasingly more lenient significance thresholds and determined the percentages of explained HRV in the same four cohorts (Supplementary Fig. 5; Table 3). Maximal variance explained by the polygenic risk score was 0.8–1.4% for SDNN, 0.9–2.3% for RMSSD and 0.9–2.3% for pvRSA/HF. This was reached at relatively small numbers of SNPs (≤71) with additional SNPs adding more noise than signal.

Table 3

Explained variance in HRV traits in the Lifelines (n=12,101), NESDA (n=2,118), TRAILS-Pop (n=1,191) and ABCD (n=1,094) cohorts by the optimal polygenic risk scores computed at the P value threshold that explained the largest percentage of phenotypic variance.

Trait	Risk score	Cohort	P cutoff	No. SNPs	P value	ΔR2
SDNN	SDNN	Lifelines	<5E-7	13	6.8E-27	0.82%
		NESDA	<5E-8	6	2.6E-08	1.16%
		TRAILS-Pop	<5E-5	64	1.1E-04	1.23%
		ABCD	<5E-5	71	9.4E-05	1.39%
SDNN	RMSSD	Lifelines	<5E-8	8	2.4E-23	0.71%
		NESDA	<5E-6	23	1.2E-07	1.05%
		TRAILS-Pop	<5E-8	8	1.2E-04	1.23%
		ABCD	<5E-7	13	2.8E-06	2.00%
SDNN	pvRSA/HF	Lifelines	<5E-8	7	3.1E-19	0.58%
		NESDA	<5E-8	4	3.5E-05	0.64%
		TRAILS-Pop	<5E-7	6	9.7E-06	1.61%
		ABCD	<5E-5	67	9.2E-04	1.01%
RMSSD	SDNN	Lifelines	<5E-7	13	8.9E-31	0.95%
		NESDA	<5E-8	6	1.6E-10	1.46%
		TRAILS-Pop	<5E-8	7	8.3E-06	1.63%
		ABCD	<5E-5	71	1.6E-04	1.30%
RMSSD	RMSSD	Lifelines	<5E-7	12	2.8E-30	0.94%
		NESDA	<5E-7	10	2.7E-10	1.43%
		TRAILS-Pop	<5E-7	11	3.4E-07	2.13%
		ABCD	<5E-7	13	3.8E-07	2.34%
RMSSD	pvRSA/HF	Lifelines	<5E-8	7	1.4E-25	0.78%
		NESDA	<5E-8	4	3.6E-09	1.25%
		TRAILS-Pop	<5E-7	6	3.7E-08	2.47%
		ABCD	<5E-8	67	8.4E-04	1.02%
pvRSA/HF	SDNN	NESDA	<5E-8	6	1.1E-12	1.52%
		TRAILS-Pop	<5E-8	7	5.0E-05	1.36%
		ABCD	<5E-5	71	5.4E-04	1.09%
pvRSA/HF	RMSSD	NESDA	<5E-7	10	5.6E-14	1.69%
		TRAILS-Pop	<5E-7	11	3.3E-06	1.78%
		ABCD	<5E-7	13	1.9E-06	2.06%
pvRSA/HF	pvRSA/HF	NESDA	<5E-8	4	4.4E-13	1.58%
		TRAILS-Pop	<5E-7	6	1.6E-07	2.25%
		ABCD	<5E-5	67	1.6E-03	0.90%

NA, not available.

NOTE: Weighted polygenic risk score was determined based on independent SNPs in the stage 1 meta-analysis. For NESDA and TRAILS-Pop the weights (that is, effects sizes) and P values were adjusted by analytically extracting the cohort’s effect size and s.e. from the meta effect size and s.e., respectively, and recalculating the P value based on these adjusted effect size and s.e., since these cohorts were included in stage 1.

The total variance explained by common SNPs (SNP-based heritability) estimated by Genomic Restricted Maximum Likelihood or LD score regression analysis varied between 10.8 and 13.2%, with only small differences in estimates across methods and HRV traits (Supplementary Note 2).

Generalization to other ethnicities

In data from up to 11,234 Hispanic/Latino individuals, five SNPs in five of the eight loci identified for RMSSD, seven SNPs in six of the seven loci for SDNN and three SNPs in three of the five loci for pvRSA/HF showed a statistically significant association that was consistent in direction with the association in individuals of European ancestry (Table 4). In data from 6,899 African–Americans, four SNPs from four of the eight loci were associated with RMSSD, three SNPs in three of the seven loci with SDNN and none with pvRSA/HF. In the combined meta-analysis in a maximum of 71,675 participants from all ethnicities, one SNP (rs6123471 on chromosome 20) was no longer significant (Table 4).

Table 4

Meta-analysis results for the identified loci in other ethnicities and combined meta-analysis results with European ancestry.

Correcting HRV for heart rate

The strong inverse association between HRV and heart rate reflects the well-established simultaneous biological effect of cardiac vagal activity on heart rate15 and HRV16, but it also expresses a mathematical dependency of the variance in inter beat interval (IBI) on the mean IBI that is unrelated to the underlying biology. We conducted three analyses to test whether the association of the HRV SNPs was robust to correction of the HRV traits for heart rate (Supplementary Table 8). First, we used a recently developed analytical technique17 to obtain the meta-analysis for the coefficient of variation of SDNN and RMMSD from the summary statistics of the HRV and resting heart rate meta-analyses18. The coefficient of variation detects the amount of IBI variability relative to the mean IBI of each subject, and deals with the proportionality-based dependence of HRV on heart rate19. Second, we established the effect of the 17 HRV SNPs on the coefficients of variation for SDNN and RMSSD in the Lifelines, NESDA and TRAILS-Pop cohorts, and meta-analysed the results. Third, we use a mediation analysis in these same cohorts to see how much of the SNP effects on the three HRV measures was mediated by heart rate. In all three analyses, we find some attenuation of the HRV SNP associations. The average mediation of the association by heart rate was ∼28%. However, the correction for heart rate left most of the HRV SNP associations intact, particularly in the first analysis that used the full discovery sample.

Association of the HRV SNPs with resting heart rate

Because the HRV traits reflect cardiac vagal activity, we expected the HRV SNPs to have an effect on resting heart rate. We performed a lookup of the 17 HRV SNPs in a GWAS meta-analysis on resting heart rate in 85,787 individuals18. Out of the 17 HRV lead SNPs, 11 were associated with heart rate after correcting for multiple testing (Supplementary Table 9, panel a). All effects were in the expected direction such that the HRV decreasing allele was associated with higher heart rate (Supplementary Fig. 6). Six of the HRV SNPs were not significantly associated with heart rate, including our top hit on chromosome 19 (rs12974991 in NDUFA11: p RMSSD=4.6 × 10−46; p heart rate=0.18). Analysis of summary statistics of the HRV and heart rate meta-analyses as implemented in the gtx R package showed that multi-SNP genetic risk scores for HRV were significantly associated with heart rate (Supplementary Table 9, panel b).

In addition, genetic risk scores based on the independent genome-wide significant HRV SNPs from the combined stage 1 and 2 meta-analysis were tested for association with heart rate in the Lifelines, NESDA, TRAILS-Pop and ABCD cohorts (Supplementary Table 9, panel b). The three multi-SNP risk scores of the HRV traits explained a small, but mostly significant percentage of variance in heart rate (0.09–1.13%). Polygenic risk score analysis showed that adding HRV SNPs below the genome-wide significance threshold did not further increase the variance explained in heart rate (Supplementary Fig. 5; Supplementary Table 9, panel c).

The reverse question, whether SNPs with effects on heart rate are associated with HRV, was also investigated. The 21 heart rate SNPs identified by the GWAS meta-analysis on heart rate18 explained between 0.2 and 0.9% of the variance in the three HRV traits (Supplementary Table 10).

Association with cardiometabolic traits and diseases

In addition to heart rate, we examined the association of the 17 HRV-associated SNPs with other confirmed risk factors for cardiac, metabolic and renal disease traits and endpoints using data from large-scale GWAS meta-analyses (Supplementary Table 11). Multi-SNP risk scores were computed based on our 17 top SNPs and we tested their association with the outcomes.

No effects of risk scores using the 17 HRV SNPs were observed for systolic or diastolic blood pressure, body mass index, renal function, heart failure, sudden cardiac death, coronary artery disease, atrial fibrillation or type 2 diabetes. Only for atrial fibrillation we observed individually significant SNPs. These two highly significant SNPs (rs10842383 near LINC00477, P=3.45 × 10−7 and rs2680344 in HCN4, P=4.34 × 10−7 (Supplementary Table 12) had large opposite effects on atrial fibrillation, while both decreased HRV. In addition to these lookups that were restricted to genome-wide significant SNPs, we employed bivariate LD score regression20 that uses the full GWAS summary statistics of the HRV and cardiometabolic traits and diseases to compute genetic correlations. The genetic correlations systematically pointed to an overlap in the genetic variants causing low HRV and increased risk for disease (that is, negative correlations with systolic and diastolic blood pressure, coronary artery disease, heart failure, sudden cardiac death, BMI and type 2 diabetes) compatible with clinical relevance of the HRV SNPs identified, although significance was reached only for systolic and diastolic blood pressure after correction for number of outcomes tested (Supplementary Table 11).

Potential functional impact of the HRV variants

To identify functional variants tagged by the 17 HRV SNPs, we performed various post-GWAS annotation (Supplementary Fig. 7). In silico annotation (Supplementary Data 2) showed that the lead SNP for SDNN on chromosome 19 was a non-synonymous SNP (rs12980262 in NDUFA11) and that the lead SNPs for RMSSD (rs12974991) and pvRSA/HF (rs12974440) were in perfect LD with this SNP (Table 1; Supplementary Data 2). SNP rs129080262 was characterized as deleterious, with a sorting intolerant from tolerant (SIFT) score of 0.01 and a polymorphism phenotyping (PolyPhen) score of 0.753 indicating a possibly damaging effect. Functional variant analyses using RegulomeDB confirmed that rs12980262 and rs12974440 in NDUFA11 on chromosome 19 likely have functional consequences (Supplementary Table 13) by binding to transcription factors or influencing the chromatin state. SNP rs6123471 in the locus on chromosome 20 was in high LD with two non-synonymous SNPs in the KIAA1755 gene (rs3746471 [r2=0.94] and rs760998 [r2=0.55]) that are predicted to yield tolerated, benign amino acid changes (Supplementary Data 2).

We examined if the 17 HRV SNPs were eQTLs in a large whole blood database. Four of the HRV SNPs were significantly (false discovery rate <5%) associated with gene expression in blood (Supplementary Table 14): rs1812835 with expression of NEO1, rs4899412 with expression of RGS6, and rs180238 and rs4262 with expression of GNG11. These four SNPs were all in strong LD with the top eQTL SNPs for these genes (r2>0.70) and lost significance after conditioning on the corresponding top eQTL. The eQTLs for NEO1 and RGS6 were replicated in at least one other whole blood eQTL study (Supplementary Table 14). The eQTL for GNG11 was replicated in the medulla (P=2.8 × 10−4) and the anterior tibialis artery (P=8.1 × 10−9). None of the 17 SNPs reached significance in a smaller heart eQTL database.

Nine of the 17 HRV SNPs were in high LD (r2>0.70) with SNPs associated with methylation level of one or multiple CpG sites (methylation quantitative trait loci [mQTLs]) in whole blood (Supplementary Table 15). Two of the HRV SNPs that were eQTLs also influenced methylation of the same gene in whole blood, strongly suggestive of a regulatory function for those SNPs. eQTL rs1812835 in NEO1 was associated with methylation level of cg11357013, cg19281068, cg11552023 and cg17150474. eQTL rs4262 was associated with methylation level of cg08038054 and cg06439941 in GNG11. The other two eQTL SNPs did not achieve genome-wide significance level for an association with methylation, but eQTL rs4899412 in RGS6 was in high LD with a proxy SNP (rs2238280) that was associated with methylation level of cg19493789, which is located in a CpG island shelf near RGS6.

Five other HRV SNPs were (in high LD with) mQTLs but were not themselves eQTLs. For example, HRV SNPs rs12974991, rs12974440 and rs12980262 (chromosome 19) were associated with methylation level of multiple CpG sites (cg22854549, cg03715305 and cg19211619) located in or nearby NDUFA11, but were not associated with expression level of NDUFA11 in whole blood. Such mQTLs may well exert a regulatory effect on NDUFA11 in other tissues. DEPICT tissue enrichment analysis (Supplementary Data 3; Supplementary Table 16, Supplementary Fig. 8) showed NDUFA11 expression was weak in blood, but enriched in heart, sensory and endocrine tissues.

Discussion

This meta-analysis of GWAS for HRV yielded 17 lead SNPs (11 independent) in eight loci that were genome-wide significantly associated, six of which generalized to individuals of African-American and Hispanic/Latino ethnicity. Various ways that correct HRV for its mathematical dependency on resting heart rate attenuated the SNP effects, but largely left the associations intact. Together, the hits in the eight loci explained 0.9–2.6% of the variance in resting HRV in four independent cohorts of European ancestry. Details of known biological functions of the genes closest to these loci are given in the Supplementary Note 6.

We noted a strong enrichment of our HRV loci in a previously conducted meta-analysis of GWAS for resting heart rate18, a known risk factor for cardiac morbidity and mortality2122. SNPs in five of the 21 resting heart rate loci (that is, LINC00477 (C12orf67), SYT10, GNG11, HCN4 and KIAA1755) were associated with HRV at genome-wide significance level and six more attained nominal significance, with associations always in the expected direction. Genetic risk scores for HRV traits were also significantly associated with heart rate and LD score regression confirmed that the allelic variants that decrease HRV in parallel increase heart rate. This suggests to us that part of the HRV SNPs exert their effect on heart rate through oscillatory modulation of pacemaker activity by the vagal nerves.

Supplementary Fig. 9 depicts the two routes by which acetylcholine released by the vagal nerves in the sinoatrial node is known to influence heart rate, both of which are supported by our results in GNG11, RGS6 and HCN4. By binding the muscarinic type 2 receptor (M2R) and dissociating the G-protein heterotrimer (Gαβγ) into a Gαi/o subunit and a Gβγ component, acetylcholine inhibits the ongoing depolarization of the pacemaker cells by β1/β2-adenylatecyclase activation of funny (If) channels and calcium channels23. In parallel, it acts to actively hyperpolarize the pacemaker cells by activation of the GIRK1/4 channel. Each route accounts for about half of the tonic decrease in heart rate upon vagal stimulation23, but the response time for M2R-GIRK effects on the sinus rate is much shorter than for the M2R-HCN2/4 or the β1/β2-adenylatecyclase signalling pathways. Only signalling through the Gβγ component is fast enough (∼0.3 s) to rapidly track changes in vagal outflow to the sinoatrial node, for example, as they occur within the duration of a single respiration (∼4.5 s), whereas signalling through the α subunit is too slow (>3 s) to track such phasic changes in acetylcholine release124. GIRK signalling, therefore, accounts for most of HRV due to the phasic oscillation in vagal activity24, but it accounts for only half of the tonic vagal effects on heart rate.

The above leads to HRV only partially capturing the vagal effects on heart rate. Additional reasons for the imperfect relation between HRV and vagal effects on heart rate12 are individual differences in: (i) resting respiration rate and depth; (ii) the amplitude of the intrinsic 0.1 Hz oscillations related to both vagal and sympathetic blood pressure regulation through the baroreflex loops; (iii) mechanotransduction or intracellular pathways stimulated by sinoatrial stretch or (iv) the efficiency of the actual vagal gating process. These processes can have a strong impact on HRV, but less so on mean heart rate. We found six SNPs in four loci, including our top hit (rs12974991 in NDUFA11), that may act on the individual differences in these processes as they had no discernible effect on heart rate, in spite of their significant impact on HRV.

The genome-wide significant SNPs in GNG11, RGS6 and NEO1 were eQTLs and in strong LD with the top mQTLs and eQTLs for the corresponding genes. Two of these (GNG11, RGS6) readily provide a biological hypothesis to account for the associations detected in the meta-analysis. The C alleles of rs4262 and rs180238 of GNG11 coding for the γ11 subunit of the heterotrimeric G-protein complex Gαβγ cause decreased expression of this subunit and were associated with lower HRV. The effects of the GNG11 eQTLs associated with lower HRV are likely to lower the availability of the γ11 subunit, thereby reducing Gβγ component-induced GIRK activation. This potentially blunts the heart rate change in response to the oscillatory changes in cardiac vagal activity.

The regulator of heterotrimeric G-protein complex signalling, type 6 (RGS6) gene on chromosome 14 was found to be linked to three independent signals for SDNN and RMSSD. RGS6 acts as a critical negative regulator of M2R signalling in the sinoatrial node of the heart rapidly terminating Gβγ signalling and thus curtailing vagal lowering of the heart rate2526. The results of our meta-analysis are consistent with a role for RGS6 in decreasing HRV previously hinted at by animal experimentation2327 and a human case report2728. The T allele of our eQTL RGS6 SNP (rs4899412) causes increased expression of RGS6. By increasing RGS6 expression, the T allele acts as a gain-of-function mutation that gives rise to a decrease in GIRK-channel signalling and the observed decrease in HRV. Of note, Rgs6−/− mice, that show the expected increase in HRV, are characterized by a strong bradycardia and an increased susceptibility to AV block and atrial fibrillation which is attributed to an enhancement of GIRK-induced sinoatrial and atrioventricular node hyperpolarization by removing the negative regulation of Gβγ by RGS6 (refs 23, 26, 28).

The association of the rs2680344 SNP in HCN4 is puzzling because HCN signalling does not involve the fast M2R-GIRK channels and cannot translate rapid vagal fluctuation into beat-to-beat variation in IBI length, that is, HRV. The effect of the HCN4 SNP on HRV may be secondary to its effects on the average slope of the diastolic depolarization29. The HCN4 protein is a key component of the If channel303132 that generates the pacemaker potential by a gradual depolarization of the sinoatrial myocyte cell membrane during diastole. This ‘pacemaker depolarization’ phase is known to be slowed by loss-of-function mutations in the HCN4 that lead to lower heart rate31 and the If is the known site of action for ivabradine and other therapeutic agents used to slow heart rate in angina patients32. Of note, both ivabradine treatment33 and loss-of-function mutations increase the risk for atrial fibrillation34. In contrast, gain-of-function mutations in the sensitivity of HCN4 for cAMP lead to higher heart rate30. This leads us to hypothesize that the A allele of rs2680344 in HCN4 either is itself a gain-of-function mutation or tags such a mutation because it increases heart rate18.

High HRV is associated with lower morbidity and mortality in patients with cardiovascular disease5, hypertension7, end-stage renal disease8 and diabetes9, but also in apparently healthy individuals1112. Using LD score regression on meta-GWAS summary statistics from various risk factors and endpoints we find some evidence for overlap in the genetic variants causing low HRV and increased risk for disease, but significance was reached only for systolic and diastolic blood pressure after correction for multiple outcomes tested. These genetic correlations are compatible with causal effects of cardiac vagal control in the aetiology of disease, but they could also be ascribed to reversed causality, where the disease process leads to lower cardiac vagal control. A strength of this study in this regard is that analyses were confined to individuals in good cardiac health, that is, cohorts excluded patients with existing cardiovascular diseases or medication potentially impacting HRV. Because we selected individuals in good cardiac health reverse effects of disease on HRV seem less likely, although some latent pathology could have been present. However, an alternative explanation that is harder to rule out is that the genetic correlation derives from pleiotropic effects of genetic variants common to both outcomes.

Further strengths of this study were the consistency of results across the different HRV traits used to capture cardiac vagal control and the generalization of the HRV SNP effects to different ancestries, in spite of known ethnic differences in absolute resting HRV35. Results also held in men and women separately and across a very large range of mean cohort ages spanning from early childhood to the late middle ages; in spite of a strong reduction in HRV values with aging36. Although effects of age and sex on HRV were taken into account in the analyses, many other factors were not. The ideal design would have corrected for the known effects of respiration depth and rate on HRV, which are independent of vagal activity37. These could not be added as covariates because they were not available in most cohorts. We were liberal in excluding other covariates like BMI, smoking and exercise in the GWAS analyses. These traits are substantially heritable themselves and adjusting for heritable covariates can bias the genome-wide association effects38 or even induce non-existing associations through collider bias39. Finally, instructions on pre-ECG recording behaviours like physical activity, and caffeine, alcohol or nicotine use were not rigorously standardized across cohorts.

Direct clinical relevance of most current GWAS findings is still low and our study is no exception. Potential future clinical use of our findings hinges on the ability of our genetic variants to capture (sub)cortical, brainstem and medullary transmission of tonic vagal activity to the sinoatrial node, not just the impact of that activity on heart rate. Subcortical generation of tonic vagal activity is an important biomarker for cardiovascular health and potentially modifiable by interventions on psychosocial stress40 and lifestyle habits41. It can even be a transdiagnostic biomarker for psychopathology and executive cognitive functioning possibly by reflecting the integrity of prefrontal cortex functioning42. Genetic markers for HRV may prove useful as instrumental variables in Mendelian randomization43 to test causal hypotheses on the effects of centrally generated vagal activity on behavioural and health outcomes.

In conclusion, this meta-analysis detects a critical role for genetic variation in Gβγ and HCN signalling in explaining individual differences in HRV. The HRV variants detected can help guide further investigations of the functional consequences and potential therapeutic implications of individual differences in sinoatrial Gβγ signalling.

Methods

Study cohorts

Appropriate IRB approval and informed consent from participants in all participating cohorts was obtained. Full information on consent procedures and details of the IRB boards are provided in the Supplementary Note 8.

HRV measurement

In this study, we investigated three HRV traits: SDNN, the RMSSD and pvRSA or HF. SDNN and RMSSD were derived from the IBI time series obtained from the R waves in the ECG4. HF was calculated from Wavelet or Fourier decomposition with power obtained from a high frequency band of either 0.15–0.40 Hz or 0.15–0.50 Hz. A time domain measure of RSA was derived by pvRSA using a respiratory signal co-registered with the ECG. Estimates of pvRSA are obtained by subtracting the shortest IBI during heart rate acceleration in the inspiration phase from the longest IBI during heart rate deceleration in the expiration phase.

HRV traits were extracted from the IBI time series preferably based on 2–10 min periods of ECG in a standardized setting, at rest and in a sitting/supine position. If ambulatory data were available, we advised cohorts to extract a period of sitting still in the evening, when this proved feasible. Supplementary Table 2 lists the actual way HRV was assessed by the participating cohorts. For the cohorts analysed in stage 2, we extended our HRV measurements to include cohorts with 10 s and/or 20 s ECG recordings, as RMSSD and SDNN based on these ultra-short recordings have shown a good agreement with 4–5 min recordings3. Furthermore, since IBI time series require reliable detection of the R-wave only, a three-lead ECG was considered sufficient while the use of more leads was encouraged. For pvRSA, an additional respiration signal of sufficient quality to detect beginning and end of inspiration and expiration was needed.

SDNN and RMSSD have prevailed in epidemiological studies because they are more easily assessed in large cohorts and, as noted above, can be obtained even from short ECG recordings. HF and pvRSA were available in fewer cohorts, but they better reflect the cardiorespiratory coupling that drives the oscillatory modulation of vagal effects in the sinoatrial node. In the typical resting respiratory frequency range, these time- and frequency-domain measures of RSA are much less contaminated by oscillations in cardiac sympathetic control than SDNN (and other measures of HRV that span a broader frequency range). This is due to the temporal dynamics of the sinoatrial node signalling pathway that acts as a low pass filter allowing only oscillations in vagal effects to translate into HRV, whereas for sympathetic effects or vagal effects at progressively higher respiratory frequencies the node acts as a leaky integrator causing more tonic changes in heart rate1. Phasic modulation of vagal effects is therefore captured most purely by pvRSA or HF. Because pvRSA and HF are conceptually similar and highly correlated with each other (r>0.80) across a wide range of values for respiration and heart rates44, we grouped the analyses on pvRSA and HF under the label pvRSA/HF.

Study population

Cohorts that had data on at least one of the three HRV traits and genome-wide data were invited to participate in the first (discovery) stage of the Genetic Variance in Heart Rate Variability (VgHRV) consortium. The stage 1 discovery analysis was performed in up to 28,700 individuals of European ancestry from a maximum of 20 cohorts. Independent cohorts with either genome-wide or gene-centred array data or with the ability to perform wet-lab genotyping on the single-nucleotide polymorphism (SNPs) taken forward from the first stage were included in the second (replication) stage. This stage included additional data from up to 24,474 individuals from 11 cohorts of European ancestry (see Supplementary Tables 1–4 for cohort descriptions and details).

Association analysis: stage 1 (discovery)

The following exclusion criteria were applied a priori: (1) individuals with heart disease (for example, angina, past myocardial infarction, left ventricular failure) and (2) individuals known to use antidepressants (particularly tricyclic antidepressants) and all anticholinergic agents (for example, digoxin, atropine and acetylcholinesterase inhibitors) because of the strong effects that these drugs have on HRV. Individuals reporting over the counter use of anticholinergic agents were not excluded.

Imputation of SNPs was done to extend and create similar SNP databases between cohorts using different genotyping platforms. Most of the cohorts used the HapMap Phase II release 22 CEU panel as reference, but later releases (for example, release 24) or other reference data sets (for example, 1000 Genomes) were also used (Supplementary Table 4).

Each cohort performed linear regression analyses on all available HRV traits using an additive SNP model adjusting for age at the time of ECG recording, sex, principal components—to adjust for population stratification—and other study-specific parameters; all HRV traits were log-transformed because of the skewness of their distributions. Only autosomal associations were examined. Analyses were performed for all individuals as well as for men and women separately.

Stage 1 meta-analysis

Before meta-analysis, quality control of all uploaded cohort files was performed using the QCGWAS package45. In case of issues the cohorts were notified and problems were solved. Using the QCGWAS results, specific imputation quality and allele frequency thresholds were set for each cohort.

An inverse-variance, fixed-effects meta-analysis was performed for RMSSD and SDNN for which SNPs of the different cohorts were merged based on rs-id. For pvRSA/HF, we performed a sample size weighted meta-analysis using z-scores with METAL46, since we combined results of two HRV phenotypes (pvRSA and HF) that have different units and ranges, and therefore incomparable SNP effect sizes. To get an idea of the size of the SNP effect on pvRSA/HF, we obtained effect sizes and s.e.’s from an additional fixed-effect meta-analysis on the GWAS results of the (majority of) cohorts that measured HF. Results of the meta-analyses were double genomic control corrected47 to control for potential inflation as a result of population stratification within and between cohorts. The results included all SNPs that met the following selection criteria: (a) a minor allele frequency in the meta-analysis of >1%, and (b) present in at least one third of the cohorts. This resulted in 2,555,913 SNPs being analysed for SDNN, 2,534,714 SNPs for RMSSD and 2,628,894 SNPs for pvRSA/HF. For each trait separately, SNPs with a P<1 × 10−6 were clumped for LD using pairwise LD checking in SNAP48 to ascertain independent primary and secondary signals (r2<0.1). A total of 23 lead SNPs in 14 loci were selected for follow-up in the second (replication) stage.

Stage 2 meta-analysis

Stage 2 cohorts applied the same exclusion criteria and performed the same association analysis as in the discovery stage, but analyses were restricted to the 23 lead SNPs. If a SNP was not available in a cohort, the best available proxy was used instead based on strongest LD according to the 1000 Genomes database. To verify homogeneity of the results in the stage 2 cohorts with those in the stage 1 cohorts, the stage 1 meta-analysis effect sizes of the 23 SNPs were correlated to the effect sizes obtained in each cohort for each of the HRV traits. If a negative correlation (r<0) was found, the cohort/trait pair was excluded from stage 2 analysis. For this reason results from one cohort for SDNN were excluded. The replication results were then meta-analysed per trait using an inverse-variance fixed-effects meta-analysis for RMSSD and SDNN and a sample size P weighted meta-analysis using z-scores in METAL46 for pvRSA/HF. SNPs were matched based on rs-id. Next the association results from both stages were combined in the same way. A SNP was only considered to be significantly associated to HRV if it satisfied the following criteria: (1) it had P<1 × 10−6 in stage 1, (2) it had a one-sided P<0.05 in the stage 2 meta-analysis congruent with the direction of effect in the stage 1 meta-analysis and (3) it had a genome-wide significant P<5 × 10−8/3 (two-sided) in the combined meta-analysis of stage 1 and 2 results, correcting for the testing of three separate traits.

Conditional analysis

In the discovery stage independent SNPs were selected for follow-up based on LD clumping (r2<0.1). To confirm independence between these SNPs within the loci on chromosome 14 and 15, we applied the conditional-and-joint analysis as implemented in the Genome-wide Complex Trait Analysis software package49 to the stage 1 summary statistics of RMSSD and SDNN with the genotype data of the NESDA cohort50 of 1,925 individuals as the LD reference data set. In addition, cohort-level individual data on log-transformed RMSSD and SDNN of 12,101 individuals from the Dutch Lifelines cohort51 were analysed using linear regression analysis with age and sex as covariates conditioned on the other associated SNP(s) within the locus.

Gene-based association analysis (VEGAS)

We performed gene-based testing with the full set of ∼2.5 M HapMap SNPs from GWAS results of all three phenotypes, using VEGAS (Supplementary Table 17). This software has the advantage of accounting for LD structure and the possibility to define a range beyond the gene bounds to include promoter, 5′UTR, intronic and 3′UTR regions into the analysis. We defined a 50 kb extra window beyond the genes, considered every SNP in this window for the gene-based analysis, and ran the analyses per chromosome with up to 106 permutations. A P<2.5 × 10−6 (=0.05/∼20,000 genes) was considered as the threshold for significance.

The Lifelines and NESDA cohorts were used for genetic risk score and polygenic risk score analyses to determine the percentage of variance explained by independent HRV SNPs that were genome-wide significant, and by SNPs meeting increasingly lenient significance thresholds, respectively. Lifelines and NESDA represent examples of a population-based cohort and a cohort ascertained on case-control status (for major depressive disorder). Both recruited adult participants. To test the stability of explained variance across the life span, we repeated this analysis in two other Dutch cohorts, the adolescent TRAILS-Pop cohort52 (age 10–18) and the ABCD cohort consisting of young children (age 5–7) (ref. 53).

For the genetic risk score, stage 1+2 summary statistics were used for the selection of HRV SNPs. No correction was needed for ABCD as genotyping in this cohort had finished only after completion of the meta-analyses. However, the NESDA cohort had been included in both stage 1 and 2, TRAILS-Pop in stage 1, and Lifelines in stage 2, so the effect sizes and s.e.’s of the HRV SNPs were corrected to subtract the effects of those cohorts to obtain independent validation cohorts54. Also, only SNPs were used in the genetic risk score if they remained genome-wide significant after analytically subtracting these cohort’s effects from the meta-analysis. Genetic risk scores of the remaining SNPs (Lifelines: SDNN(9), RMSSD(7), pvRSA/HF(5); NESDA: SDNN(9), RMSSD(11), pvRSA/HF(5); TRAILS-Pop: SDNN(10), RMSSD(11), pvRSA/HF(5)) weighted by the adjusted effect size were calculated for the participants of all four cohorts and regressed on the three HRV traits (pvRSA/HF was not available in Lifelines). Explained variance was computed as the change in R2 from a model with and without the genetic risk score, while adjusting both for age, sex and principal components.

To compute the polygenic risk scores, the imputed genotypes were first converted to best-guess genotypes. This was done regardless of the imputation quality, since it was previously shown that even low-quality SNPs might contribute to the variance explained by SNPs (ref. 54). The SNP set was further pruned for LD using PriorityPruner (http://prioritypruner.sourceforge.net/) to select independent SNPs, taking the significance of the SNP in the discovery meta-analysis of each of the HRV traits into account. This provided three LD-pruned SNP sets. Polygenic risk scores were then calculated in PLINK55 using significance thresholds of 5 × 10−8, 5 × 10−7, 5 × 10−6, 5 × 10−5, 5 × 10−4, 0.005, 0.05, 0.5 and 1 and associated with the three HRV traits and resting heart rate in the Lifelines, NESDA, TRAILS-Pop and ABCD cohorts. For NESDA and TRAILS-Pop pruning and polygenic risk score analysis was based on analytically corrected results, since these cohorts were part of stage 1 of our study54.

Heritabilities and genetic correlations

We applied genomic restricted maximum likelihood analysis implemented in the Genomic Complex Trait Analysis software package56 in the Lifelines cohort (Supplementary Table 18) to estimate the percentages of additive phenotypic variance that can be explained by common SNPs (that is, common SNP heritability denoted as h2SNP). For this analysis, SNPs from the HapMap Phase 3 project were selected to obtain a set of independent SNPs. We further used LD score regression to estimate the heritabilities of the three HRV traits and the genetic correlation among HRV traits and with heart rate20. The GWAS meta-analysis summary statistics for RMSSD, SDNN and pvRSA/HF were obtained from stage 1 of the current study, and the GWAS meta-analysis summary statistics for heart rate from the discovery stage of a recent GWAS meta-analysis for heart rate18. The LD scores required by the method were computed using 1000 Genomes data of Europeans. The heritabilities of the three HRV measurements were estimated using the univariate model of this method. Cross-phenotype LD score regression analysis was performed using the LDSC tool (LD SCore) to estimate genetic correlations between pairs of phenotypes20.

In addition, we used the Oman Family Study57 to perform univariate and bivariate analyses in five multigenerational highly inbred pedigrees to estimate the heritabilities for and the genetic correlations between log-transformed RMSSD, SDNN, HF and heart rate using SOLAR (v7.2.5) (ref. 58).

We further examined the generalization of loci identified after meta-analysis of stage 1 and 2 results to other ethnicities using data from 11,234 individuals of two Hispanic/Latino cohorts, and 6,899 individuals from five African-American cohorts (Supplementary Tables 1–4). Stage 3 meta-analyses were performed in the same way as in stage 2 of this study to assess the effect of the HRV-associated SNPs in individuals of Hispanic/Latino and African-American ancestry, in the combined set of European and Hispanic/Latino ancestry, in the combined set of European and African-American ancestry, and in all three ethnicities combined. Here, we applied the same criteria for significance as in stage 2 described above, that is, a SNP was only considered to be significantly associated to HRV if: (1) it had P<1 × 10−6 in stage 1 meta-analysis in European individuals, (2) it had a one-sided P<0.05 in the new ethnicity specific meta-analysis congruent with the direction of effect in the stage 1 meta-analysis in European individuals and (3) it had a genome-wide significant P<5 × 10−8/3 (two-sided) in the combined meta-analysis.

The well-known inverse association between HRV and heart rate in part reflects a dependency of the variance in IBI on the mean IBI that is unrelated to cardiac vagal activity59. That is, the slower the heart rate, the longer the IBI, and therefore, any proportionally minor beat-to-beat differences in IBI are more pronounced at slower heart rates. This occurs on top of the well-established dual effect of cardiac vagal activity that lowers heart rate and increases HRV1516. Although these two mechanisms (biological, mean-variance dependency) are impossible to completely separate, we conducted three analyses to test whether the HRV SNPs were robust to correction for the mean IBI.

First, we corrected SDNN and RMSSD for their dependency on mean IBI by using the coefficient of variation, which is a more parsimonious solution19 than the logarithmic approach suggested by Monfredi et al.29. We obtained the summary statistics for the resting heart rate GWAS meta-analysis18 from: https://walker05.u.hpc.mssm.edu/ and used the GWIS procedure17 to infer a GWA analysis of the coefficient of variation of the SDNN and the RMSSD. We approximated the coefficients of variation by (SDNN/X) × 100% and (RMSSD/X) × 100%, respectively, where X equals 60,000 per heart rate. Transformation from heart rate to IBI is required as both terms in the coefficient of variation (HRV and IBI) are in milliseconds, whereas the heart rate GWAS meta-analysis used heart rate in beats per minute. As the coefficients of variation were skewed, we used a log-transformation. As an example of the linear approximation by GWIS we assume that the increaser effect of one allele for an SDNN SNP is +0.2 with the same SNP reducing heart rate by −0.1. Given a mean SDNN of 100 and mean heart rate of 60, we can then approximate (omitting some nuances adequately explained in Nieuwboer et al.17) the effect of the SNP on the coefficient of variation of the SDNN as:

We used the delta method to approximate a s.e. for the effect of the SNP given that we know the s.d.’s for the SNP effects on SDNN and HR, and their dependence. We obtain the dependence from analysis with LD score regression20.

Second, we performed association analyses for our 17 top SNPs on the actual log-transformed coefficients of variation of SDNN and RMSSD computed in the Lifelines, NESDA and TRAILS-Pop cohorts and then meta-analysed these results. Because pvRSA and HF are expressed on different scales, such a meta-analysis was not feasible for pvRSA/HF.

Third, we repeated the association analysis for our 17 top SNPs on SDNN, RMSSD and pvRSA/HF in the Lifelines, NESDA and TRAILS-Pop cohorts with and without adjusting for heart rate as a covariate and performed mediation tests with the Sobel test to assess the mediation effect of heart rate on the HRV SNP association. Significance of the Sobel t-value was determined using a bootstrap procedure (n=10,000 permutations). The mediation P values of the three cohorts for SDNN and RMSSD and two for pvRSA/HF (as this was not available in Lifelines) were next meta-analysed to determine the significance of mediation and to compute the percentage of the SNP effect on HRV that was mediated through its effects on heart rate. We note that this is likely an overcorrection because the HRV SNPs are expected to influence heart rate through a common biological mechanism, that is, changes in cardiac vagal activity.

Association of the HRV SNPs with heart rate

We conducted a lookup of the 17 (11 independent) HRV lead SNPs identified in this study using the results of a recent GWAS meta-analysis for heart rate18. A HRV-associated SNP was considered to be significantly associated with resting heart rate if the GWAS meta-analysis result for heart rate was <0.05/11=0.0045. Three separate HRV weighted multi-SNP genetic risk scores were calculated from 10 (SDNN), 11 (RMSSD) and five (pvRSA/HF) HRV SNPs, respectively, (based on all genome-wide significant SNPs for the respective HRV trait in the stage 1+2 meta-analysis). These were tested for their effect α on resting heart rate using the gtx package in R (https://cran.r-project.org/web/packages/gtx), which approximated α by (Σω × β × se−2)/(Σω2 × se−2) with se≅✓(1/Σω2 × se−2), where ω is the effect of the SNP on HRV, β is the effect of the SNP on heart rate and se is the s.e. of β. This approximation requires only single SNP association summary statistics extracted from GWAS results60. The effects of the multi-SNP genetic risk scores were considered as statistically significant when the P was less than 0.0045 (correcting for 11 traits; heart rate and the 10 cardiometabolic traits described below).

In addition to these lookups that were restricted to genome-wide significant SNPs, we employed LD score regression20 that uses the full summary statistics of the HRV and heart rate GWAS meta-analyses to compute genetic correlations.

We further examined the variance in resting heart rate explained by multi-SNP genetic risk scores (based on the lead SNPs only) and of the full polygenic risk scores for HRV in the four Dutch cohorts Lifelines, NESDA, TRAILS-Pop and ABCD. The identical approach was used as done previously for the computation of variance explained in the HRV traits themselves.

Association of heart rate SNPs with HRV

We also performed reverse analyses to detect the effects of heart rate SNPs on the HRV traits. In our GWAS meta-analysis results for SDNN, RMSSD and pvRSA/HF, we performed a lookup for the 21 previously identified heart rate SNPs by den Hoed et al.18. A heart rate associated SNP was considered to be significantly associated with HRV if the P was <0.05/21=0.0024. The 21 heart rate SNPs were tested in a multi-SNP risk score for their effect on the HRV traits using the gtx approach as described above.

To examine the variance explained in the HRV traits by the 21 heart rate SNPs, multi-SNP genetic risk scores and polygenetic risk scores based on the heart rate SNPs were computed in the Lifelines, NESDA, TRAILS-Pop and ABCD cohorts and these were tested for association with the available HRV traits. For the multi-SNP genetic risk scores weights were either the original SNP effect sizes on heart rate (for NESDA, TRAILS-Pop and ABCD) or corrected because of participation of the cohort in the GWAS meta-analysis (Lifelines). Only 15 of the 21 SNPs were used in the Lifelines cohort because five SNPs lost genome-wide significance after subtracting the SNP effects of the Lifelines cohort. One other SNP (rs826838) was removed because it was in LD (r2=0.15) in Lifelines with a more significant heart rate SNP (rs7980799).

We estimated the joint effect of the HRV SNPs on cardiometabolic and renal disease traits and endpoints. The traits included were systolic and diastolic blood pressure, body mass index and urinary albumin excretion as well as estimated glomerular filtration rate based on creatinine. The clinical outcomes used were heart failure, coronary artery disease, atrial fibrillation, sudden cardiac death and type 2 diabetes. The relevant consortia (Supplementary Table 11) and/or corresponding authors of the studies were contacted with the request to perform a lookup and provide summary GWAS meta-analysis results for our list of 17 SNPs.

The association analyses consisted of the same three steps as used for heart rate. First, we checked the P of our HRV SNPs (or their proxies) in the cardiometabolic trait or disease GWAS meta-analysis results. Second, three separate HRV weighted genetic risk scores were calculated from 11 (RMSSD), 10 (SDNN) and five (pvRSA/HF) HRV SNPs, respectively (based on all genome-wide significant SNPs for the respective HRV trait in the stage 1+2 meta-analysis). These were tested for their effect on the clinical outcomes using a regression model in the gtx package in R as described above for the association of the HRV SNPs with heart rate. The effects of the genetic risk scores were considered as statistically significant when the P was less than 0.0045 (0.05/11, correcting for heart rate and the 10 traits and diseases).

In addition to these lookups that were restricted to genome-wide significant SNPs, we employed LD score regression20 that uses the full GWAS summary statistics of the HRV and cardiometabolic traits and diseases to compute genetic correlations.

Search for known functional SNPs (in silico annotation)

We followed an in silico bioinformatics-based approach61 to search and annotate SNPs in the regions surrounding the 17 identified HRV SNPs. For this purpose SNP positions were converted from National Center for Biotechnology Information (NCBI) build 36, Human Genome 18, to NCBI build 37, Human Genome 19, (GRCh37/hg19) using the NCBI Genome Remapping service tool (http://www.ncbi.nlm.nih.gov/genome/tools/remap). For ±1 Mb regions surrounding the SNPs, we downloaded the according variance call format file from the 1000 Genomes Project. We used data of 503 European ancestry individuals from 1000 Genomes Project Phase 3 (version 5.a.) to calculate LD between the HRV SNP and all other SNPs within the area. SNPs in moderate to high LD (r2≥0.5) were subsequently selected and annotated by ANNOVAR software62 for functionality. For all non-synonymous SNPs loss-of-function and gain-of-function was determined by using the SIFT and PolyPhen prediction scores. A SNP was categorized as deleterious if the SIFT score was ≤0.05 or the PolyPhen score was between 0.957 and 1 (probably damaging).

We used RegulomeDB to integrate results from the RoadMap Epigenomics and ENCODE projects to identify variants that are likely to have functional consequences using the lead SNPs identified for the three HRV traits. We distilled information on transcription factor binding and chromatin states for SNPs that showed most evidence of being functional, that is, for SNPs with a RegulomeDB score <4.

Finally, all the HRV SNPs and those that were in high LD (r2≥0.8) with them were looked up in the National Human Genome Research Institute GWAS catalogue to check for association with other complex traits or diseases identified in previous GWAS studies63.

eQTL analyses

We performed expression quantitative trait locus (eQTL) analysis in whole blood to identify regulatory variants that were associated with the HRV SNPs using the gene-expression database from NESDA50 and NTR64 cohorts. The sample used for this analysis consisted of 4,896 individuals of European ancestry. For complete details on the sample and the procedures, see65.

eQTL effects were tested with a linear model approach using MatrixeQTL66 with expression level as dependent variable and SNP genotype values as independent variable. In this study we only tested cis effects for our HRV SNPs, meaning that the probe was at a distance <1 Mb from the SNP on the genome according to GRCh37/hg19. For each probe set that displayed a statistically significant association with at least one SNP in the cis region, we identified the most significantly associated SNP (top eQTL). Conditional eQTL analysis was carried out by first residualizing probe set expression using the corresponding top eQTL and then repeating the eQTL analysis using the residualized data.

All HRV SNPs with significant results in the NESDA/NTR eQTL data were looked up in two other independent whole blood eQTL databases, eQTLs in lymphoblastoid cell lines, eQTLs in 10 different brain regions, and a heart eQTL database.

mQTL analyses

We obtained mQTL results from a previously published study67. In short, genome-wide DNA methylation data were generated using Illumina 450 k arrays for 3,841 whole blood samples. Corresponding genotype data were imputed using the Genome of The Netherlands68 reference panel. To determine the effect of nearby genetic variation on methylation levels (cis-mQTLs) a Spearman rank correlation and corresponding P value was computed for each CpG-SNP pair, in which the CpG and SNP location were no further than 250 kb apart. To control for multiple testing, we used a permutation procedure to empirically control the false discovery rate at 5%. The distribution of observed P values was compared to the P value distribution obtained from the analyses on permuted data. For a permutation the sample identifiers of the genotype data set were shuffled, breaking the link between the genotype data set and the methylation data set. This was repeated 10 times to obtain a stable distribution of P values under the null hypothesis. To determine the false discovery rate only the strongest effect per CpG in both the real analysis and in the permutations were selected.

Gene prioritization using four bioinformatics approaches

Potentially causal genes for the associations identified by GWAS were identified using four previously described bioinformatics tools: ToppGene, Endeavour, MetaRanker and DEPICT (Supplementary Table 19). To this end, we first retrieved positional coordinates for all lead SNPs according to GRCh37/hg19. These coordinates were used to extract all genes located within ±40 kb of lead SNPs using the UCSC genome browser. The identified genes subsequently served as input for ToppGene and Endeavour, together with two genes with established roles in sinus node function (HCN4) and synaptic signal transmission (ACHE) that served as training genes. For MetaRanker, we first combined results of the stage 1+2 meta-analyses of GWAS for the three HRV traits, retained the association with the lowest P for lead SNPs that were identified for multiple traits, and subsequently provided SNPs, P values, and the same two test genes (HCN4 and ACHE) as input. For DEPICT—arguably the most powerful and informative of the four methods—we used results from the stage 1 meta-analysis for all SNPs that reached a P for association <10−5 as input, for each of the three HRV outcomes separately. In order for genes to be prioritized by the combined four approaches, they needed to be either: (1) selected by DEPICT for at least one of the three HRV outcomes; or (2) identified by at least two of the three remaining tools (ToppGene, Endeavour and/or MetaRanker).

Network and functional enrichment analyses

We performed gene network and enrichment analysis using the GeneMANIA algorithm, which uses data resources on genetic interactions, protein–protein, co-expression, shared protein domains and co-localization networks. To build a functional interaction network, we selected genes as input for this analysis using the following criteria: (a) genes implicated by gene prioritization using the four bioinformatics approaches described above, (b) the genes closest to our 17 HRV SNPs, (c) genes to which linked (r2>0.50) non-synonymous SNPs mapped, (d) genes to which other linked (r2>0.80) SNPs mapped, (e) genes identified by VEGAS, and (f) expression probe gene names significantly associated with HRV eQTLs (false discovery rate <0.01). The input gene list was extended to 100 by their most strongly interacting genes and a weighted composite functional association network was constructed61. Subsequently, functional enrichment analysis of all genes of the constructed interaction network against Gene Ontology (GO) terms was performed to find the most enriched GO terms (Supplementary Table 20). Significantly enriched GO terms (false discovery rate <0.10) were visualized as highlighted boxes within their corresponding GO tree depicted by the RamiGO R package69 (Supplementary Fig. 10).

Tissue and gene-set enrichment analyses

We used DEPICT for a tissue enrichment analysis to tabulate tissues that are enriched for expression of genes located within ±40 kb of SNPs with a P<10−5 association with the HRV traits. DEPICT calculates the likelihood of every known gene to be a member of, amongst others, KEGG, GEO or REACTOME-based gene sets (N=14,461) to create reconstituted gene sets. It then determines which of these reconstituted gene sets are enriched for the HRV genes. A graphical representation of DEPICT’s reconstituted gene-set enrichment analysis (P<0.05 after Bonferroni correction for examining three HRV traits) was generated using a script that is based on an affinity propagation clustering algorithm by Frey et al.70. Interactions between gene sets are considered significant if the Pearson coefficient, which is based on the number of genes that are shared between gene sets, is >0.3.

Data availability

Summary statistics of the meta-analyses are available on request from the corresponding authors after a formal data access request procedure and approval by the VgHRV consortium.

Additional information

How to cite this article: Nolte, I. M. et al. Genetic loci associated with heart rate variability and their effects on cardiac disease risk. Nat. Commun. 8, 15805 doi: 10.1038/ncomms15805 (2017).

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.