Genetic insights into resting heart rate and its role in cardiovascular disease

Higher resting heart rate (RHR) is associated with cardiovascular diseases and all-cause mortality in traditional epidemiological studies^1,2,3,4,5. However, RHR is influenced by disease status and a plethora of potential confounders, which could affect these associations.

A Mendelian randomization (MR) approach, in which genetic variants associated with RHR are used as a proxy for RHR, can also be used to study the association of RHR with cardiovascular diseases and all-cause mortality. Since genetic variants are fixed from conception and then randomly assigned from parents to offspring, they are more immune to reverse causation and confounders⁶. In our previous study, we found evidence for a positive association between genetically predicted RHR and all-cause mortaliy⁷. However, a higher genetically predicted RHR was not found to increase the risk of cardiovascular diseases^7,8,9 and appeared to decrease the risks of atrial fibrillation and cardio-embolic stroke⁹. Multiple genome-wide association studies (GWAS) have identified RHR-associated genetic variants which could be used as genetic instruments in MR studies to entangle the relationship between RHR and cardiovascular disease^7,8,10. The two largest GWAS to date have been performed in the UK Biobank^7,10 but were either performed in a subcohort of individuals with available genetic data during its time of publication⁷ or lacked a replication cohort¹⁰. We set out to perform a genome-wide meta-analysis of RHR in the largest sample size to date with internal replication of the associated genetic variants, in order to broaden our knowledge of the biological mechanisms underlying interindividual differences in RHR and identify robustly associated RHR-associated genetic variants to assess the genetic association of RHR with cardiovascular disease and all-cause mortality.

In this work, a genome-wide meta-analysis of 100 studies in up to 835,465 individuals reveals 493 independent RHR-associated genetic variants in 352 loci, including 68 genetic variants outside previously identified RHR-associated loci (Fig. 1a). In addition, in silico analysis pinpoint that identified candidate genes are mainly enriched in cardiomyocyte tissue (Fig. 1b), and Mendelian randomization analyses indicate that higher genetically predicted resting heart rate is not associated with all-cause mortality, increases the risk of dilated cardiomyopathy, but decreases the risk of developing atrial fibrillation, ischemic stroke, and cardio-embolic stroke (Fig. 1c).

Genome-wide meta-analysis of resting heart rate

We performed a meta-analysis of RHR GWAS using 99 cohorts consisting of up to 351,158 individuals, which from here on will be referred to as the International Consortium of Resting Heart Rate (IC-RHR). Second, we performed a GWAS on 484,307 individuals from the UK Biobank. These large cohorts were meta-analyzed to include up to 835,464 individuals, in whom 30,458,884 directly genotyped and imputed autosomal variants were analyzed (Supplementary Data 1, Fig. 1a). The meta-analysis revealed 493 independent genetic variants in 352 loci. Out of these 493 independent genetic variants, 68 were outside previously identified RHR-associated loci and 67 of those were internally replicated (Fig. 2a, Supplementary Data 2)^7,8,10. Out of the 425 genetic variants inside previously identified RHR-associated loci, a total of 376 were internally replicated. In addition, 332 out of the 352 loci were considered internally replicated as they showed the concordant direction of effects and nominal associations (P < 0.01) in the UK Biobank GWAS and IC-RHR meta-analysis (Fig. 2b). The RHR-associated genetic variants identified previous studies from Eppinga et al. and Den hoed et al. were all replicated in the current study (Supplementary Data 3). A total of 74 loci identified in the study by Guo et al. were not replicated in the current study, of which 40 would not have been identified as loci using the current GWAS clumping criteria. The remaining 34 loci did not reach genome-wide significance in the current meta-analysis with generally high P-values in the IC-RHR consortium, probably therefore failing replication (Supplementary Data 3). The linkage disequilibrium (LD) score regression intercept of the meta-analysis was 1.051 ± 0.002, suggesting little evidence of genomic inflation (Fig. 2c). The QQ plots for the UK Biobank GWAS and IC-RHR meta-analysis are shown in Supplementary Fig. 1. The genomic control lambda’s, LD-Score intercepts and the attenuation ratio statistics suggested no inflation due to non-polygenic signals^11,12. Single nucleotide polymorphism (SNP) heritability of RHR, as calculated by LD score regression, was estimated to be 10%. A polygenic score weighted by the effect sizes of the IC-RHR explained 5.33% of the variation in RHR in the UK Biobank. A Chow test indicated the absence of strong differences between participants with a history of any cardiovascular disease or use of RHR-altering medication versus participants without such a medical history (Supplementary Data 4). Genetic correlation analyses were performed and we observed significant correlations with anthropometric measurements, pulse wave reflection index, and physical activity measurements (Supplementary Data 5). A query of the GWAS Catalog showed that the 493 genetic variants associated with RHR were most commonly in high LD (LD > 0.8) with anthropometric measurements and blood pressure traits (Supplementary Data 6).

Candidate genes and insights into biology

We explored the potential biology of the 352 RHR-associated loci by prioritizing candidate genes in these loci (Supplementary Data 2). A total of 407 unique genes were in close proximity to the lead variant, defined as the nearest gene and any additional gene within 10 kb (Supplementary Data 2). There were 52 genes that contained coding genetic variants in LD (R² > 0.80) with RHR lead variants. Functional annotation of these coding variants is provided in Supplementary Data 7. Using summary data-based Mendelian randomization analysis (SMR) and heterogeneity in dependent instruments (HEIDI) tests, we found that the RHR-associated loci and eQTLs colocalized at 88 genes (Supplementary Data 8)¹³. Lastly, 381 unique genes were taken forward by data-driven expression-prioritized integration for complex traits (DEPICT) analyses (Supplementary Data 9). Of the 670 unique candidate causal genes identified, 33 genes were prioritized by at least three out of four established methods, which may be used to prioritize candidate genes (Fig. 2d). Of these genes, PHACTR4, ENO3, and SENP2 were prioritized by all four methods. Full annotations of all identified genes are in Supplementary Data 10.

Pathway analyses and tissue enrichment

Pathway analysis performed by DEPICT showed that RHR revolved around mainly cardiac biology, including cardiac tissue development, muscle cell differentiation, and pro-arrhythmogenic pathways. A total of 1471 reconstituted gene sets within 155 gene clusters were significantly associated with RHR (FDR < 0.05). The newly discovered gene clusters consisted of mostly protein–protein interaction pathways and were commonly located in the periphery of the network (Supplementary Data 11, Supplementary Fig. 2). The tissue enrichment analysis by DEPICT showed 28 tissues at FDR < 0.05 and implicated the cardiovascular system as the most important tissue type, with 8 of the 10 most significantly enriched tissues located within the cardiovascular system (Fig. 3, Supplementary Data 12). Non-cardiovascular tissues with enrichment included muscle and fat tissues, the adrenal glands, the esophagus, and urogenital structures. Conditional analyses showed that associations with non-cardiovascular tissues were rather due to co-expression of RHR genes in cardiovascular tissue than independent enrichment of RHR-associated genes in non-cardiovascular tissues (Fig. 3).

ECG morphology

The ECGenetics browser, which contains genome-wide summary statistics of every time-point of the complete cardiac cycle at a resolution of 500 Hz, was used to gain insights into the electrophysiological effect of the RHR genetic variants¹⁴. A total of 86 genetic variants were strongly associated with at least one ECG time point on the non-normalized and normalized association patterns across the full RR interval at a stringent Bonferroni-corrected P-value < 1 × 10⁻⁷. The associations represented a plethora of ECG morphologies (Supplementary Data 13, Supplementary Fig. 3). The ACHE, ANKRD1, and SCN5A genes exhibited their largest electrical effects on atrial depolarization, BAG3 and TTN on ventricular depolarization, and RGS6 and SYT10 on ventricular repolarization. The ECG-wide MR highlighted several loci that had not been associated with resting heart rate or cardiac rhythm and structure previously. The CCLN1 gene exhibited strong effects on atrial depolarization, RAP1A and ZBTB38 exerted strong effects on early and late ventricular repolarization, respectively. The ECG-wide MR showed that RHR variants exert the largest effect on ventricular repolarization on both the non-normalized and normalized association patterns (Supplementary Fig. 4).

Single-nucleus RNA expression

Single-nucleus RNA sequencing data obtained from a healthy human heart revealed that RHR gene expression is highest in ventricular cardiomyocytes, followed by atrial cardiomyocytes (Supplementary Data 14)¹⁵. The candidate genes of genetic variants involved in non-isoelectric parts of the ECG showed stronger expression patterns than the isoelectric parts, for example, those involved in left atrial depolarization (ANKRD1), ventricular depolarization (FOHD3, RBM20, MYO18B, TTN) and ventricular repolarization (CACNA1C) (Supplementary Fig. 3).

Mendelian randomization analyses

Series of two-sample MR analyses were performed to test whether genetically predicted RHR is associated with all-cause mortality and cardiovascular diseases (definitions provided in Supplementary Data 15, 16). We initially used the inverse variance weighted multiplicative random-effects (IVW-MRE) model, which provides a consistent estimate under the assumption of balanced pleiotropy. If we found evidence for a genetic association using the IVW-MRE model, we further interrogated these findings using several sensitivity analyses that are more robust to different sources of bias in MR analyses.

First, we assessed the association between genetically predicted RHR and all cause-mortality in the UK Biobank participants over a median follow-up of 8.9 years (interquartile range 8.2–9.5) (Supplementary Data 17–19). Genetically predicted RHR was not associated with the risk of all-cause mortality (HR 1.024, 95% CI 0.993–1.057, P = 0.13), as shown in Fig. 4a. We did not find evidence for an association between genetically predicted RHR and parental longevity. Neither did we find evidence for an association between RHR and the 35 leading causes of mortality in the UK Biobank. Systematic alteration of key differences between the current and previous Mendelian randomization study indicated that the most likely cause of the discrepancy between these studies arises from false positive findings in previous one-sample MR analyses caused by weak-instrument bias at lower P-value thresholds (Supplementary Fig. 5, Supplementary Data 20)⁷. Non-linear MR analysis showed an always-increasing dose–response relation between genetically predicted RHR and all-cause mortality that was compatible with a null effect (Fig. 4b, Supplementary Data 21, 22), providing no evidence for a U-shaped pattern that has been previously described³.

We then explored the association between genetically predicted RHR and several prevalent cardiovascular diseases. We did not find evidence for an association between genetically predicted RHR and coronary artery disease in the UK Biobank (OR 0.977, 95% CI 0.946–1.009, P = 0.16) or in the CARDIoGRAMplusC4D cohort (OR 0.976, 95% CI 0.944–1.010, P = 0.17), in line with previous analyses^8,9. Similarly, there was no evidence for an association between genetically predicted RHR and myocardial infarction in the UK Biobank or in the CARDIoGRAMplusC4D cohort (Fig. 5, Supplementary Data 23–26). We found no evidence for non-linear dose-response relations of genetically predicted RHR with coronary artery disease or myocardial infarction (Supplementary Data 21, 22).

Higher genetically predicted RHR was suggestively associated with a lower risk of atrial fibrillation development in the UK Biobank (OR 0.946, 95% CI 0.897–0.998, P = 0.04) and in the AFgen consortium (OR 0.942, 95% CI 0.897–0.989, P = 0.02), but these results were not significant after correction for multiple testing (P < 4.17 × 10⁻³). MR-Lasso, which can provide evidence for potential causal associations when there is a small number of genetic variants with heterogeneous ratio estimates, indicated that genetically predicted RHR was significantly inversely associated with atrial fibrillation (Fig. 5, Supplementary Data 23). The contamination mixture model, which provides evidence for potential causal associations if the plurality of the genetic instruments is valid, provided evidence for a negative association between genetically predicted RHR and atrial fibrillation in the UK Biobank cohort, but this was not replicated in the AFgen consortium (Fig. 5). Non-linear MR analyses showed a significant negative exponential growth pattern in the dose-response relation between genetically predicted RHR and atrial fibrillation (Supplementary Data 21, 22, Supplementary Fig. 8). Specifically, individuals at the extreme right tail of the distribution of genetically predicted RHR had a lower risk of atrial fibrillation. For example, compared with the population mean RHR of approximately 70 bpm, individuals with a genetically predicted RHR of 89 and 98 bpm had a significantly lower risk of atrial fibrillation (OR 0.969, 95% CI 0.941–0.998, P = 0.04; OR 0.922, 95% CI 0.897–0.948, P = 6.36 × 10⁻⁹), while this was not true for a genetically predicted RHR of 80 bpm (OR 1.000, 95% CI 0.968–1.034, P = 0.99).

We found that higher genetically predicted RHR is associated with the risk of any stroke (OR 0.951, 95% CI 0.926–0.976, P = 1.59 × 10⁻⁴), ischemic stroke (OR 0.940, 95% CI 0.915–0.967, P = 1.08 × 10⁻⁵) and cardio-embolic stroke (OR 0.875, 95% CI 0.828–0.925, P = 2.11 × 10⁻⁶), suggestively associated with large artery stroke (OR 0.939, 95% CI 0.884–0.998, P = 0.04) and not with small vessel stroke (OR 1.001, 95% CI 0.950–1.055, P = 0.97) in the MEGASTROKE consortium. The results were consistent across MR methods for any, ischemic and cardioembolic stroke (Fig. 5, Supplementary Data 23). The associations between genetically determined RHR and any stroke or ischemic stroke could not be replicated in the UK Biobank using a univariable MR IVW-MRE approach (OR 0.987, 95% CI 0.953–1.023, P = 0.49; OR 0.970, 95% CI 0.928–1.015, P = 0.19). We found no evidence for a non-linear association between genetically determined RHR and any ischemic stroke (Supplementary Data 21, 22). We used a multivariable MR approach to gain insights into potential mediating factors or pleiotropic pathways in the association between genetically predicted RHR and stroke. First, we found that the direct effects of RHR on cardio-embolic stroke to be attenuated by the effects of atrial fibrillation and estimated the attenuation through atrial fibrillation to be 18.4% (Fig. 6, Supplementary Data 27, 28). The direct effects of genetically predicted RHR on any stroke and ischemic stroke were most strongly attenuated by pulse pressure, with estimated attenuation of 28.1% and 31.5%, respectively (Fig. 6a, b, Supplementary Data 27, 28). There was a strong association between genetically predicted RHR and pulse pressure (β = −0.192, SE = 0.019, P = 1.81 × 10⁻²⁴), but MR-Steiger sensitivity analysis filtered a large part of the genetic variants and repeating the MR on the remaining subset did not show a significant association between RHR and pulse pressure (β = −0.005, SE = 0.008, P = 0.57, Supplementary Data 29, 30).

Lastly, we found that genetically predicted RHR is associated with an increased risk of dilated cardiomyopathy in the UK Biobank (OR 1.391, 95% CI 1.205–1.605, P = 6.27 × 10⁻⁶). The results were robust to MR-Lasso (OR 1.411, 95% CI 1.228–1.622, P = 1.20 × 10⁻⁶) and MR-contamination mixture models (OR 1.697, 95% CI 1.318–2.184, P = 4.03 × 10⁻⁵). We excluded 96 variants associated with the Q-R upslope at −18 ms of the R peak (P < 0.05), which has been established as a biomarker for dilated cardiomyopathy¹⁴, to investigate whether reversed causation contributed to the association with dilated cardiomyopathy. The results were similar to the main analyses (Supplementary Table 1). We did not find evidence for an association between genetically predicted RHR and heart failure, heart failure excluding cardiomyopathies, and hypertrophic cardiomyopathy (Fig. 5, Supplementary Data 23). We did not find evidence for a non-linear association between genetically determined RHR and any type of heart failure (Supplementary Data 21, 22). Scatterplots and dose-response curves of the association between RHR and all assessed outcomes can be found in Supplementary Figs. 6–10.

We assessed whether the Wald estimates between the RHR-associated genetic variants and the cardiovascular diseases could identify risk loci not anticipated to be associated with these outcomes in the outcome GWASs. The locus FOXC1 for coronary artery disease, USP39 for myocardial infarction, and SLC35F1 and SSPN for atrial fibrillation were significantly (P < 1.01 × 10⁻⁴) and concordantly associated with their respective outcomes in both cohorts while not reaching genome-wide significance in either one of the outcome cohorts (Supplementary Data 31).

We report 493 genetic variants in 352 loci associated with RHR, discovered in the largest GWAS meta-analysis of RHR to date in up to 835,465 individuals^7,8,16,17. This increase in sample size allowed us to report 68 novel RHR-associated genetic variants and, importantly, provide internal replication for 376 genetic variants previously associated with RHR. A total of 670 candidate genes were prioritized, providing a comprehensive data catalog for future studies on RHR and offering potential new insights into its biology. Four strategies were employed to prioritize candidate genes and PHACTR4, ENO3, and SENP2 were highlighted by all four strategies. The PHACTR4 gene regulates protein phosphatase 1 which interacts with actin and is involved in processes ranging from angiogenesis to cell cycle regulation¹⁸. It has been associated with pulse pressure and systolic blood pressure in a previous GWAS analysis¹⁹. ENO3 encodes beta-enolase, which plays an important role in glycolysis and striated muscle development²⁰. It has been implicated in cardiac myocyte development through its function in energy metabolism in both humans and rats^21,22. SENP2 encodes sentrin-specific protease 2, which deconjugates small ubiquitin-related modifiers 1 and 2 that are involved in regulating posttranslational modification of a wide variety of proteins that affect a multitude of different cellular processes^23,24. Several of these affected proteins are critical in cardiac development and mouse models have shown that alterations in SENP2 activity lead to congenital heart defects^25,26. Involvement of SENP2 in a multitude of cellular processes is reflected by its implication in GWAS of various conditions, including systolic and diastolic blood pressure²⁷, type 2 diabetes²⁸, the conduction system of the heart^29,30 and estimated glomerular filtration rate³¹. The loci we associated with coronary artery disease (FOXC1), myocardial infarction (USP39), and atrial fibrillation (SLC35F1 and SSPN) through their effects on RHR have been associated with these cardiovascular diseases in recent studies^32,33,34,35.

To obtain further biological insights into RHR, we performed pathway analyses using DEPICT and found numerous newly associated pathways. The strongest associated clusters were identified previously and their importance to RHR biology was therefore validated in the current study⁷. Conditional analyses on the tissue enrichment demonstrate that genes influencing RHR are more likely co-expressed than primarily or solely located within non-cardiovascular tissues. However, it should be noted that conditional analyses inherently attenuate tissue enrichment considering DEPICT is based on the co-regulation of gene expression³⁶. Using cardiac single-nucleus RNA data, we demonstrate that RHR genes are mostly expressed in cardiomyocytes. We provide electrophysiological insights into the biology of the RHR-associated variants and show that they exert diverse effects on ECG morphology with the largest effect on ventricular repolarization.

In-depth analyses were performed to assess genetic associations of RHR with clinical outcomes. In contrast to previous observational^1,2 and MR studies⁷, we do not find evidence for an association between genetically predicted RHR and all-cause mortality. Moreover, genetically predicted RHR was not associated with parental longevity nor with any of the 35 leading causes of mortality. Lack of such associations suggests that follow-up length or large heterogenetic effects of RHR on different causes of mortality are unlikely causes for the absence of an association between genetically predicted RHR and all-cause mortality³⁷. We demonstrate that the most likely cause of the discrepancy between current and previous results arises from false positive findings in previous one-sample MR analyses that were caused by weak-instrument bias at lower P-value thresholds³⁸. We hypothesize that RHR is not on the causal pathway to mortality itself and that previous observational studies are more likely to reflect confounders, such as stress and socio-economic status, or reversed causation, in which an individual’s disease status increases both RHR and mortality risk^39,40.

The linear MR between RHR and atrial fibrillation provided suggestive evidence for an inverse relationship between RHR and atrial fibrillation, in line with a previous linear MR study⁹. We do find a significant negative exponential dose-response curve between RHR and atrial fibrillation in support of an inverse relationship, and take the non-linear MR forward as the main result considering the fractional polynomial test indicated that a non-linear model fitted the localized average causal effect estimates better than the linear model. Previous observational studies on the relationship between RHR and atrial fibrillation have shown conflicting results and have described various relationships including inverse linear^41,42,43, U-shaped⁴³, and J-shaped⁴⁴ associations. All these association patterns support the hypothesis that individuals with a low RHR might exhibit a higher risk of atrial fibrillation development compared to those with an average RHR. A recent stratified Mendelian randomization showed an inverse genetic relationship between RHR and atrial fibrillation in individuals with an RHR below 65 bpm as well⁴⁵. Possible mechanisms that could underly an increased risk of atrial fibrillation in individuals with a low RHR include increased left atrial stroke volume and consequent atrial remodeling due to myocyte stretching⁴⁶, or an increased vagal tone promoting global disorganization in the left atrium due to increased heterogeneity of the refractory period⁴⁷. In contrast to the often hypothesized U-shaped or J-shaped association^43,44, we find a decreasing risk of atrial fibrillation development in those with a high RHR. One potential explanation is that previous observational studies were affected by collider bias through confounding factors which increase atrial fibrillation risk and typically occur in tandem with a high rather than a low RHR, such as hypertension⁴⁸ and obesity⁴⁹. We advocate for a cautious interpretation of current results due to the diverse biological mechanisms through which the RHR-associated genetic variants alter the risk of atrial fibrillation development⁸.

We found that genetically predicted RHR was inversely associated with risk of any, ischemic and cardio-embolic stroke. The results were not replicated in the UK Biobank, possibly due to the substantially lower amount of cases. The inverse association is in contrast to many observational studies and we therefore performed multivariable MR analyses to pinpoint biological mechanisms that could underly the discrepancy^4,5. We showed that atrial fibrillation attenuates the protective effect of higher genetically predicted RHR on developing cardio-embolic stroke. This indicates either biological or mediated pleiotropic effects of atrial fibrillation in the association between genetically determined RHR and cardio-embolic stroke, which cannot be distinguished based on the current results. The relationship between a low RHR and cardioembolic stroke was attenuated by only 18.4% when corrected for atrial fibrillation, which might underestimated due to atrial fibrillation being a commonly missed diagnosis. Correction for atrial fibrillation only minimally affected the association between RHR and any or ischemic stroke, despite cardio-embolic stroke accounting for a substantial amount of ischemic stroke cases⁵⁰. Although hypertension is another important risk factor for stroke, it commonly occurs in tandem with a higher and not a lower RHR^4,51,52 and we found that neither systolic nor diastolic blood pressure affects the association between RHR and stroke. We did find that pulse pressure attenuates the association between RHR and any, ischemic and large-artery stroke. Lower RHR has previously been demonstrated to increase pulse pressure due to a higher likelihood of pressure wave reflections during prolonged systole⁵³ and increased pulse pressure has been established as a risk factor of stroke^53,54,55. Moreover, the Conduit Artery Functional Endpoint Study (CAFE) study postulated that pulse pressure underlies the inferiority of β-blocker based treatment (which lowers RHR) to amlodipine-based treatment in the prevention of stroke, despite equal effects on peripheral blood pressure^53,56. Our results could be considered as support for this mechanism in the scenario that the RHR-associated genetic variants only affect pulse pressure through RHR and the association with stroke is primarily driven by RHR. However, MR-Steiger sensitivity analysis indicated that the association between the RHR-associated genetic variants and pulse pressure is unlikely mediated through RHR entirely. Biological pleiotropic effects are therefore more likely to cause the attenuation of the association between RHR and stroke when correcting for pulse pressure.

Finally, our study provides evidence that higher genetically predicted RHR increases the risk of developing dilated cardiomyopathy. The importance of decreasing RHR in the treatment of heart failure with reduced ejection fraction, the clinical phenotype of dilated cardiomyopathy, has been thoroughly studied. Beta-blockers have been shown to reduce mortality in individuals with heart failure with reduced ejection fraction and form the cornerstone of pharmacological treatment^57,58,59. There is also evidence that ivabradine lowers cardiovascular mortality in heart failure with a reduced ejection fraction⁶⁰. This protective effect is more likely due to its effect on RHR than heart rate variability as it has a larger effect on RHR⁶¹. The fact that the MR results were robust to the exclusion of SNPs associated with the −18 ms point of the R-peak, an established biomarker of dilated cardiomyopathy, supports the interpretation that current findings are driven by RHR differences that mimic pharmacological rate control¹⁴. Our MR on the compound definition of heart failure could be hampered by its phenotypical heterogeneity, as we were unable to differentiate between heart failure with reduced and preserved ejection fraction. It would be interesting to repeat the current MR analyses if more in-depth phenotyping on left ventricular ejection fraction and function becomes available, especially considering the different effects of RHR on familial dilated versus hypertrophic cardiomyopathy in the current study.

Several limitations should be considered. Although the current 493 RHR-associated genetic variants explained more than double the RHR variance compared to the 64 loci from our previous study⁷, there is still a large gap with heritability estimates from twin studies that range between 23% and 70%^62,63,64. Future studies could include whole exome sequencing data to further increase our insights into the genetic architecture of RHR⁶⁵. Second, individuals with cardiovascular diseases were included in the GWAS, which could potentially affect exposure-outcome associations. However, post-hoc analysis showed that UK Biobank participants with a history of cardiovascular disease or who used RHR-altering medication can be jointly analyzed with participants without such a medical profile. In addition, a two-sample MR strategy was adopted, reducing the risk that potential weak-instrument bias increases type 1 error rates through the reintroduction of confounding, population stratification, or correlated pleiotropy³⁸. We note the broad biological nature of RHR genetic variants as illustrated by the diverse ECG patterns the genetic variants elicit on the full cardiac cycle. These broad effects should be taken into consideration for the correct interpretation of the MR results, as pleiotropy and reversed causation might be introduced in the MR. For example, some genetic variants were included in the MR analyses which could be more specific for another trait (i.e. rs2234962 near BAG3 for dilated cardiomyopathy). We believe that the influence of reversed causation on current results is minimal because we excluded variants more strongly associated with the outcome. The MR results were generally consistent across a multitude of sensitivity analyses, strengthening the interpretation of a true relationship. However, our study is not interventional in design, and a conservative interpretation of the results as generally unconfounded rather than causal estimates should be preferred. We stress that any causal claims can only be made if interventions or drugs alter RHR equal to the biological mechanisms in which RHR-associated genetic variants affect RHR.

In conclusion, our GWAS meta-analysis identified 493 RHR-associated genetic variants within 352 loci, to which we prioritized 670 candidate causal genes. We demonstrated cardiovascular tissues as the primary enrichment sites for RHR gene effects and showed that their gene expression is highest in cardiomyocytes. ECG signatures showed that RHR-associated genetic variants exert the largest effect on RHR through ventricular repolarization. We found no evidence for linear and non-linear associations between genetically predicted RHR and all-cause mortality across several analyses, suggesting that the well-known link between higher RHR and all-cause mortality reflects confounding factors and reversed causation. The results point towards an inverse association between genetically predicted RHR and the development of atrial fibrillation and any stroke, ischemic stroke, and cardio-embolic stroke, whereas it is positively associated with dilated cardiomyopathy development. Multivariable MR analyses showed that atrial fibrillation attenuates the protective effect of higher RHR on the development of cardio-embolic stroke. Pulse pressure attenuates the protective effects on any stroke, ischemic and large artery stroke, but this likely reflects biological pleiotropy rather than true mediation.

Method details

Populations

The full RHR meta-analysis included 100 studies with data on RHR in up to 835,465 individuals. RHR was obtained from ECG in 54 studies, from pulse rate in 31 studies (of which seven were self-measured by the participants), from blood pressure monitor in nine studies, from electronic medical records in three studies, from manual measurement in one study and through a combination of multiple of the before mentioned methods in two studies. Further information on cohort characteristics is provided in Supplementary Data 1 and statistical details are provided in the “Genome-wide association studies” section.

Imputation and quality control

Genotyping and quality control before imputation were performed using different genome-wide genotyping arrays and methods, as further detailed in Supplementary Data 1.

The UK Biobank was imputed to the Haplotype Reference v1.1 panel (HRC) by the Wellcome Trust Centre for Human Genetics. Analysis has been restricted to variants that are in the HRC v1.1. Quality control of samples and variants and imputation was performed by the Wellcome Trust Centre for Human Genetics, as described in more detail elsewhere⁶⁶.

The 99 cohorts of the IC-RHR were imputed to 1000 Genomes Phase 1 and 3. For further information, please see Supplementary Data 1.

On the cohort level, we performed quality control by (1) re-formatting and SNP-name harmonization; (2) checking the used reference panel by plotting effect allele frequency plots using 1000G as a reference; (3) checking for genomic inflation by plotting QQ plots; (4) checking the betas by plotting histograms of the beta, frequency and info; (5) comparison of the expected P-value based on beta and standard error versus reported P-values.

Association with other traits

Genetic correlation analyses with GWAS of previously investigated traits were performed using LD Hub platform⁶⁷. Genetic correlations were considered significant if they achieved a Bonferroni-corrected significance threshold of P < 0 .05/855 = 5.85 × 10⁻⁵. The GWAS Catalog was queried to find previously established genetic variants (P < 1 × 10⁻⁵) in LD (R² > 0.8) with all 493 RHR variants⁶⁸. Summary statistics were downloaded from the NHGRI-EBI GWAS Catalog on 27/04/2020.

Functional annotation of genes

For all independent genetic variants that were genome-wide significantly associated in the final meta-analysis, candidate causal genes were prioritized as followed: (1) by proximity, the nearest gene and any other gene within 10 kb; (2) protein-coding genes containing variants in LD with RHR-associated variants at R² > 0.8; (3) eQTL genes in LD with RHR-associated variants at R² > 0.8; and (4) DEPICT gene mapping using variants that achieved P < 1 × 10⁻⁸ (further information described below). Annotation of all identified genes was performed by querying GeneALacart⁶⁹.

Query of dbNSFP

The dbNSFP database (version 3a) was queried to obtain functional prediction and annotation of all potential non-synonymous genetic variants⁷⁰. The dbNSFP database contains information on multiple prediction algorithms and conservation scores further detailed elsewhere⁷⁰.

eQTL analyses

Colocalization of multiple expression quantitative trait loci (eQTL) was performed using SMR and HEIDI analyses (version 0.710)¹³ in data repositories from GTEx V7⁷¹, GTEx brain⁷¹, Brain-eMeta eQTL⁷², and blood eQTL from Westra⁷³ and CAGE⁷⁴. Colocalization analyses were performed to test whether the effect size of the RHR-associated variants on the phenotype is most likely mediated by gene expression¹³. eQTL genes were considered as candidate causal genes if they achieved a significance after Bonferroni correction for the amount of eQTLs tested (P < 0.05/188,737 = 2.65 × 10⁻⁷), passed the HEIDI test at P > 0.05, and if the lead variants of the eQTL genes were in LD (R² > 0.8) with the RHR-associated genetic variants.

DEPICT analyses

DEPICT was used to find genes associated with identified variants, enriched gene sets, and tissues in which these genes are highly expressed. DEPICT.v1.beta version rel137 (obtained from https://data.broadinstitute.org/mpg/depict/) was used to perform integrated gene function analyses as stated above³⁶. DEPICT was run using all genetic variants that achieved P < 1 × 10⁻⁸. Genes were considered possible candite genes if FDR < 0.05, taking into account multiple hypothesis testing.

Pathway analyses

DEPICT was used to find enriched gene sets using the settings described above³⁶. Enriched genesets were further clustered on the basis of the correlation between scores for all genes using an Affinity Propagation method as provided by DEPICT³⁶. Each cluster was named according to the name of the most central gene set as identified using the Affinity Propagation method. Identified meta-clusters were compared to the clusters found in the study of Eppinga et al. and were determined to be new if not a single cluster within the meta-cluster had been identified before. Clustering was performed using DEPICT software in python 2.7⁷⁵ and visualization using Cytoscape 3.8.0⁷⁶.

Tissue enrichment

DEPICT was used to find enriched tissues using the settings described above³⁶. Enriched tissues were further investigated by performing conditional analyses to provide evidence for an independent association with RHR. The following formula was used:

Here, Z_t is the maximum Z-value of all tissue Z-values, Z_s a vector of all tissue Z-values, and ρ_ts the correlation between the tissue of Z_t and the tissue of Z_s⁷⁷. The maximum Z-value (Z_t) was determined for every new iteration. Conditional analyses were performed up until the highest Z-value reached 2.585, which corresponds to the lowest Z-value with an FDR < 0.05.

ECG morphology

The ECGenetics browser was used to gain insights into the electrophysiological effect of the RHR-associated genetic variants¹⁴. Detailed information on the methodology can be found in the study of Verweij et al. and is briefly discussed below. The ECGenitics browser contains genome-wide summary statistics of the complete cardiac cycle. The complete cardiac cycle was defined using two methods, including (a) the signal-averaged electrocardiographic beat surrounding the R wave at a resolution of 500 Hz resulting in 500 averaged data points and (b) R–R intervals corrected signal (made of equal length of 500 data points).

All RHR-associated genetic variants were tested for their association with both the non-normalized and normalized association patterns. A heatmap was constructed containing all associated genetic variants associated with at least one point on the ECG at a Bonferonni-corrected P-value of 0.05/493/(500 × 2) = 1 × 10⁻⁷. Effects were aligned to the most positively associated allele across all time points. The heatmap shows a hue ranging from red (positive effect) to a blue color (negative effect) color scale, with yellow indicating no effect. Secondly, the total effect of the 493 RHR-associated genetic variants on ECG morphology was assessed using an ECG-wide MR approach (inverse variance weighted fixed-effects model) on the non-normalized and normalized association pattern.

Single-nucleus RNA expression

All genes prioritized in the current study were queried in the single-cell data from the study of Tucker et al. through the Broad Institute’s Single Cell Portal (available at: https://singlecell.broadinstitute.org/single_cell/study/SCP498/transcriptional-and-cellular-diversity-of-the-human-heart under study ID SCP49) to gain insights in their transcriptional and cellular diversity¹⁵. We selected the 86 genetic variants strongly associated with at least one ECG time point. We took forward the most likely candidate gene per genetic variant, based on the amount of gene identification strategies by which the gene was identified. When a genetic variant highlighted multiple genes identified by the same amount of gene identification strategies, we took forward the gene with the highest biological plausibility of involvement in RHR biology. A dotted heatmap was constructed for this subset of genes.

UK Biobank definitions

In the UK Biobank, we captured the prevalence and incidence of functional outcomes through data collected at the Assessment Centre in-patient Health Episode Statistics (HES) and data on the cause of death from the National Health Service (NHS) Information Centre. Prevalence of disease was also based on an interview with a trained nurse at the baseline visit (self-reported). HES data were available up to 31-03-2017 for English participants, 29-02-2016 for Welsh participants, and 31-10-2016 for Scottish. Information on cause of death was available for participants from England and Wales until 31-01-2018, and from the NHS Central Register Scotland for participants from Scotland until 30-11-2016. Definitions of all-cause mortality, 35 leading causes of mortality (defined as any cause of mortality with a prevalence >0.1%), coronary artery disease, myocardial infarction, atrial fibrillation, stroke (any stroke, any ischemic stroke), heart failure and subtypes (hypertrophic cardiomyopathy, dilated cardiomyopathy) are provided in Supplementary Data 15. Longevity was obtained through questionnaires in which participants were asked to provide the age of death of both parents. Individuals were excluded in case the answer was older than 115 years, if they reported themselves as adopted, if their parent was still alive but not yet long-lived, or if their parent died prematurely (fathers <46 years or mothers <57 years), in line with previously established methods³⁷. Combined parental longevity was assessed by summing Z-scores of the age of death from both parents if the information on both parents was provided³⁷. Systolic and diastolic blood pressure values were obtained through two automated and/or two manual blood pressure measurements. The average value of all available blood pressure measurements was used per phenotype. Automated measurements were corrected according to the previously described methodology⁷⁸. In addition, we corrected systolic and diastolic blood pressure for medication use, by adding respectively 15 and 10 mmHg to the blood pressure trait⁷⁹. Pulse pressure was calculated by subtracting diastolic from systolic blood pressure.

Statistical details of the analyses on functional outcomes are provided in the “Genetics and regression analyses on functional outcomes in the UK Biobank” and “Mendelian randomization” sections.

External cohort definitions

External cohorts included the CARDIoGRAMplusC4D⁸⁰, AFGen⁸¹, MEGASTROKE⁵⁰, and ICBP-plus¹⁹ consortia and descriptions have been detailed previously. Effect sizes from the ICBP-plus consortium were obtained from the meta-analysis of the UK Biobank and ICBP-plus, after subtracting the effects from the UK Biobank (for further details, see the section “Meta-subtract of blood pressure traits”). An overview of these studies is provided in Supplementary Data 16. We searched for proxies (LD > 0.8) in case variants could not be found within the outcome datasets.

Quantification and statistical analysis

Genome-wide association studies

All included cohorts performed genetic variant association analyses on RHR using linear regression analyses assuming an additive genetic model (Supplementary Data 1). No transformation of heart rate was performed and extreme (>4 SD) phenotypic outliers were excluded analogous to previous GWAS on RHR and as per predefined criteria⁷. The GWAS model was adjusted for age, age², body mass index, sex, and study-specific covariates (e.g. principal components, genotyping array, and RHR measuring method in case multiple RHR methods were used within a study).

The UK Biobank GWAS was performed using BOLT-LMM v2.3beta2, employing a mixed linear model that corrects for population structure and cryptic relatedness⁸². A total of 484,307 participants from the UK Biobank remained available for the GWAS after the exclusion of 14,242 individuals for whom no genetic data was available, 1341 individuals who failed genetic quality control, 1058 individuals who were outside the 4SD range for RHR, and 1587 individuals due to missing covariates (Supplementary Data 1).

Study-specific details and methodology of the 99 cohorts of the IC-RHR are provided in Supplementary Data 1. Several software programs were used for the analysis, including mach2qtl⁸³, minimac⁸⁴, minimac2⁸⁵, IMPUTE2⁸⁶, GEEPACK⁸⁷, ProbaABEL⁸⁸, and MMAP⁸⁹, for which version details per cohort are provided in Supplementary Data 1. A fixed effect meta-analysis using the inverse variance method in METAL was performed on all 99 cohorts, including up to a total of 351,158 individuals⁹⁰. Genomic control was applied at the study level by correcting for the study-specific lambda.

All genetic variants were excluded if they had poor imputation quality score (Info < 0.3) and effective sample size (N_eff) < 25 for the genetic variants computed as sample size × Info × 2 × minor allele frequency (MAF)×(1−MAF). After these exclusions, a total of 19,400,415 and 27,082,649 variants remained available for the UK Biobank and IC-RHR GWAS respectively.

Again, a fixed-effects meta-analysis using the inverse variance method in METAL was performed to pool the data from the UK Biobank and IC-RHR up to 835,465 participants using ~30M genetic variants⁹⁰. LD score regression software (v1.0.0) was used to calculate linkage disequilibrium score regression intercepts and attenuation ratios^11,12. We corrected for genomic inflation prior to the meta-analysis by multiplying the standard errors with the square root of linkage disequilibrium score regression intercepts in the UK Biobank (1.132 ± 0.017) and the IC-RHR (1.020 ± 0.010)^11,12.

PLINK (version 1.9) was used to prune genetic variants in a set of independently associated variants⁹¹. An independent genetic variant was defined as a genome-wide significant genetic variant in low LD (R² < 0.005) with another genome-wide significant variant within a five-megabase window. A genetic locus was defined as the most significant variant in a megabase region at either side of the independent genetic variant.

A single one-stage replication analysis was performed next. The following criteria had to be satisfied for a signal to be reported as a replicated signal for RHR:

1. 1.

   the sentinel genetic variant has P < 1 × 10⁻⁸ in the discovery (UKB + IC-RHR) meta-analysis;

2. 2.

   the sentinel genetic variant shows support (P < 0.01) in the UKB GWAS alone;

3. 3.

   the sentinel genetic variant shows support (P < 0.01) in the IC-RHR meta-analysis alone;

4. 4.

   the sentinel genetic variant has the concordant direction of effect between UKB and IC-RHR datasets;

The sentinel genetic variants were compared with previous loci from previous GWAS of RHR and were determined novel if located outside a 1-megabyte distance of previously RHR-associated loci^7,8,10. We selected the P-value thresholds to be an order of magnitude more stringent than a genome-wide significance P-value to ensure robust results and to minimize false positive findings.

Post-hoc quality control

We performed additional analyses to investigate whether individuals with a history of cardiovascular disease or those who took RHR-altering medication could influence the results of the GWAS. The UK Biobank population was stratified by a medical history of any cardiovascular disease or reported intake of RHR-altering medication. A history of any cardiovascular disease was defined according to the definition in Supplementary Data 15. RHR altering medication was defined as intake of beta-blockers, calcium antagonists, sotalol, amiodarone, flecainide, anti-depressants, atropine, other anti-cholinergic medication, cardiac glycosides, diuretics, ACE-inhibitors or angiotensin II receptor blockers, analogous to previous methods⁹². Linear regressions on RHR were performed in both populations, using cluster-robust standard errors with genetic family IDs as clusters to account for the relatedness among participants. Exclusions and covariates were similar to those used for the GWAS. Individuals belonging together based on 3rd-degree or closer as indicated by the kinship matrix (kinship coefficient > 0.0442) provided by UK Biobank received a family ID. A Chow test was used to investigate whether there were significant differences in beta estimates in participants with and without cardiovascular disease or RHR-altering medication⁹³. The post-hoc quality control was performed using the statistical software STATA 15 (StataCorp LP).

Genetics and regression analyses

All of the outcomes assessed in the UK Biobank that are reported in this manuscript have been adjusted for age, sex, the first 30 principal components (PCs) to account for population stratification, and genotyping array (Affymetrix UK Biobank Axiom® array or Affymetrix UK BiLEVE Axiom array). The exclusions were performed according to the above-mentioned methods for the GWAS of RHR in the UK Biobank. In addition, we excluded 74,471 individuals based on familial relatedness, after which 412,481 individuals remained available for further analyses.

SNP-outcome associations for all outcomes (see the section “UK Biobank definitions”) were obtained for all 493 variants with a P < 1 × 10⁻⁸ in the RHR GWAS (see section “Mendelian randomization analyses”). The associations with all-cause mortality and 35 leading causes of mortality within the UK Biobank (defined as a prevalence higher than 0.1%) were obtained using a Cox proportional hazard model during a median (interquartile range) follow-up of 8.9 (8.2–9.5) years. The associations with parental longevity was assessed using linear regression analyses. The associations with both prevalent and incidence of cardiovascular diseases were assessed using logistic regression analyses. Cox and linear regressions were corrected for age at baseline, while the logistic regression analysis was corrected for age until the last date of follow-up to correctly account for both prevalent and incident disease.

We performed an in-depth assessment of the association between RHR and all-cause mortality in the UK Biobank by systematically altering the differences between the current study and the previous study from Eppinga et al., which included (a) the set of SNPs, (b) the P-value threshold for SNP inclusion, (c) the assessment of the outcome in an independent cohort, and (d) the follow-up length⁷. Genetic risk scores for RHR were created following an additive model by summing the number of alleles (0, 1, or 2) for each individual after multiplication with the effect size for RHR. Genetic risk scores were constructed using the 493 discovered variants within the full meta-analyses using the effect sizes of the IC-RHR, the effect sizes of the UK Biobank, and using the 73 previously discovered variants at five P-value thresholds (1 × 10⁻⁸, 5 × 10⁻⁸, 1 × 10⁻⁷, 1 × 10⁻⁶, 1 × 10⁻⁵). These were transformed to translate to a change of 5 bpm. The association with all-cause mortality was tested using Cox regression analyses in different populations of the UK Biobank. One population included all individuals (ncases = 16,289, ncontrols = 396,183). Another population included a subset of individuals which were genotyped for the UK Biobank interim release from May 2015, which included in the GWAS by Eppinga et al. (ncases = 4953, ncontrols = 133,102). The final population consisted of a subset of individuals without genetic information at the time of the UK Biobank interim release, which was therefore not included in the GWAS by Eppinga et al., ncases = 11,336, ncontrols = 283,081). Please note that sample sizes might slightly differ from those in the previous GWAS due to updated exclusions. Lastly, point d) was taken into account by re-performing above-mentioned steps using mortality data up until the previously available follow-up (All individuals, ncases = 7099, ncontrols =405,373; individuals not included in the GWAS by Eppinga et al., ncases = 5000, ncontrols = 289,417; and those included in the GWAS by Eppinga et al., ncases = 2099, ncontrols = 115,956). All regression analyses were performed using the statistical software STATA 15 (StataCorp LP).

Mendelian randomization analyses

All 493 independent genetic variants at P < 1 × 10⁻⁸ in the final meta-analysis were taken forward in the MR. To minimize overlap between exposure and outcome cohorts, effect sizes were taken from the IC-RHR data to test the associations with outcomes within the UK Biobank, whereas effect sizes were taken from the UK Biobank to test the association within other independent cohorts. Proxies (LD > 0.8) were searched in case genetic variants could not be found within the UK Biobank or IC-RHR. All effect sizes were transformed to translate to a change in RHR of 5 bpm.

Potential weak instrument bias was assessed by calculating the F-statistic using the following equation:

In this formula, n is the sample size of the exposure and R² is the amount of variance of the exposure explained by the SNP⁹⁴. R² was calculated based on summary statistics using a previously established formula⁹⁵. Genetic variants were not excluded from further analyses if the F-statistic was <10 as this can exacerbate bias by increasing the chance of winner’s curse³⁸. Exposure and outcome summary statistics were then harmonized using the TwoSample MR package⁹⁶. Forward strand alleles were inferred using allele frequency information and palindromic SNPs were removed if the MAF was above the recommended setting of 0.42⁹⁶. MR-Steiger filtering was applied to explore pleiotropic effects through the assessment of potential reversed causation. R² for both the exposure and outcome were calculated and variants were removed from further analyses if the R² of the exposure is significantly lower (P-value < 0.05) than the R² of the tested outcome⁹⁷. R² for linear traits was calculated as described above⁹⁵, R² for binary outcomes was calculated on the liability scale⁹⁸. A true causal direction was assumed if the R² for binary outcomes was too small to be correctly estimated. Variants were excluded from further analyses in case a false causal direction was indicated.

The linear association between genetically determined RHR on all outcomes was initially assessed using the IVW multiplicative random-effects method, which provides a consistent estimate under the assumption of balanced pleiotropy. The Rücker framework was applied to assess heterogeneity and thus potential pleiotropy within the MR effect estimates⁹⁹. A Cochran’s Q P-value of <0.05 was considered as proof of heterogeneity within the IVW estimate and, as a consequence, balanced horizontal pleiotropy. An I² index >25% supports this conclusion¹⁰⁰. The MR-Egger test was performed to allow SNPs to exert unbalanced horizontal pleiotropy¹⁰¹. The Rücker framework assesses heterogeneity within the MR-Egger regression (Rucker’s Q) and calculates the difference between heterogeneity within the IVW effect estimate (Q–Q’)⁹⁹. A significant Q–Q’ (P < 0.05) in combination with a significant non-zero intercept of the MR-Egger regression (P < 0.05) was considered an indication of unbalanced horizontal pleiotropy. We then moved from an IVW model to the MR-Egger model as initial analysis, as the MR-Egger can provide causal estimates if SNPs exert unbalanced horizontal pleiotropy under the assumption that the Instrument Strength Independent of Direct Effect (InSIDE) assumption holds. Weak instrument bias in the MR-Egger regression analysis was assessed by I²_GX and was considered to indicate a low risk of measurement error if larger than 95%¹⁰². The MR-Lasso method was used to find consistent estimates under the same assumptions as the IVW method, but only for the set of genetic variants not identified as outlier¹⁰³. This method is most valuable in the scenario that a small proportion of the genetic variants is invalid and shows heterogeneous ratio estimates¹⁰³. The weighed median approach was used to provide a consistent estimate if up to half of the variants are invalid. Finally, we performed the MR contamination mixture method to provide a consistent estimate if no larger subset of invalid genetic variants estimates the same causal association than the subset of valid genetic variants¹⁰⁴.

The non-linear associations of genetically predicted RHR with all-cause mortality and cardiovascular diseases were assessed using a fractional polynomial method^105,106. These associations were assessed using the UK Biobank as an outcome cohort, considering this was the largest cohort with individual-level data available to us. Consequently, we used the independent weights of the IC-RHR meta-analysis to construct a weighted polygenetic risk score of RHR by summing the number of alleles (0, 1, or 2) for each individual after multiplication with the effect size between the genetic variant and RHR. We first calculate residual RHR by subtracting the results of the regression of RHR from the polygenetic score of RHR from RHR itself. Covariates of the regression included age, age², sex, BMI, genotyping array, and PC1-PC30, analogous to the GWAS. Residual RHR, which characterizes the predicted RHR for an individual if their polygenetic score took the value zero, was then divided into 30 quantiles. Stratifying on residual RHR rather than total RHR avoids overadjustment and collider bias as residual RHR is not downstream of the effect of the genetic variants on the outcome in a causal diagram. We then calculated the genetic associations with the exposure in each stratum of residual RHR using linear regression analyses, correcting for the same covariates as described above. Two tests for non-linearity in the genetic association with the exposure (trend and Cochran’s Q tests) were performed to investigate heterogeneity in the polygenetic score of RHR on residual RHR in different strata. We then calculated the genetic associations with the outcome in each stratum. The same methodology was used as described in the section “Genetics and regression analyses”, including the same covariate model (age, sex, genotyping array, and PC1-PC30) and regression type (Cox regression for all-cause mortality and logistic regression for cardiovascular diseases). The outcome regression coefficient was then divided by the exposure regression coefficient as a ratio of coefficients to obtain local average causal effects (LACE) in each stratum. These localized average causal effects were meta-regressed against the mean of the exposure in each stratum in a flexible semiparametric framework, using the derivative of fractional polynomial models of degrees 1 and 2. All possible fractional polynomials of degrees 1 and 2 were fitted using the powers −2, −1, 0, 0.5, 1, 2, and 3¹⁰⁶. The fractional polynomial of degree 1 is fit to the data if the fractional polynomial of degree 1 was as good of a fit (P > 0.05) as the degree 2 as indicated by the likelihood ratio test. Three tests for non-linearity of the association between genetically predicted RHR and the outcomes are reported: a trend test, which assesses for a linear trend among the localized average causal effect estimates, a Cochran’s Q test, and a fractional polynomial test, which assesses whether a non-linear model fits the localized average causal effect estimates better than a linear model. Please note that before fitting the fractional polynomials, we subtracted 45 from the values of RHR as the most flexible fit is achieved when the exposure is close to 0 but still positive. A reference of RHR of 70 bpm was taken as this was close to the mean RHR of 69.3 bpm. An additional 1390 individuals were dropped compared to the linear MR estimates obtained from the UK Biobank cohort due to missing BMI values necessary for the correction of the exposure regression coefficients.

A multivariable MR approach was used to gain additional insights into the relationship between RHR (effect sizes of the UK Biobank) and (subtypes of) stroke from the MEGASTROKE consortium. We used either atrial fibrillation (AFgen consortium⁸¹), systolic blood pressure, diastolic blood pressure, or pulse pressure (ICBP consortium, please see the section “Meta-substract of blood pressure traits” for further details¹⁹) as secondary exposures to obtain insights in the direct effect of RHR on (subtypes of) stroke that are independent of this secondary exposures¹⁰⁷. First, a multivariable MR-IVW method was used, in which for each exposure the instruments are selected and regressed together against the outcome, weighting for the inverse variance of the outcome¹⁰⁷. Weak instrument bias for any of the exposures was assessed using Q_x1 and Q_x2¹⁰⁷. When both are larger than the critical value on the χ² distribution, there is little evidence of weak instrument bias. The critical value on the χ² distribution was calculated by subtracting one degree of freedom from the amount of SNPs at a P-value of 0.05. Q_a was considered to indicate potential pleiotropy when larger than the critical value on the χ² distribution as calculated by the amount of SNPs minus two degrees of freedom at a P-value of 0.05¹⁰⁷. Multivariable MR-Egger was performed to allow for unbalanced horizontal pleiotropy¹⁰⁸. An MR-Egger intercept with a P-value < 0.05 in combination with a significant Q_a was considered proof of unbalanced horizontal pleiotropy, and consequently the MR-Egger regression to provide a robust causal estimate¹⁰⁸. Multivariable MR-Lasso analysis was performed as this method provides consistent estimates even when half of the genetic variants are invalid instruments and display unbalanced pleiotropy¹⁰⁹. We also performed multivariable weighted median analysis as this type of analysis has been shown to perform well under higher levels of pleiotropy¹⁰⁹. We did not search for proxies in the multivariable MR setting as this could introduce uncertainty through different LD patterns between the secondary exposure and outcome and we therefore re-estimate the univariable effect of RHR on the outcome with the eligible SNPs to allow for a better comparison of the results.

We assessed whether the Wald estimates between the RHR-associated genetic variants and the cardiovascular disease outcomes could identify risk loci not previously associated with these outcomes in their respective GWASs. The genetic variants were considered associated with the outcome if (a) the Wald estimates had concordant effects within the UK Biobank as well as either the CARDIoGRAMplusC4D⁸⁰, AFGen⁸¹, or MEGASTROKE⁵⁰ cohorts, (b) when the Wald estimates were significant at a Bonferonni corrected threshold of P < 1.01 × 10⁻⁴, that is, α = 0.05 with Bonferroni correction for a maximum of 493 independent tests, and (c) the genetic variant did not reach a genome-wide significant threshold of P-value < 5 × 10⁻⁸ in either one of the outcome cohorts used in the current study.

MR analyses were performed using R (version 3.6.3), the TwoSampleMR package (version 0.5.3)⁹⁶, and the MR-Lasso source code¹⁰³. The multivariable MR analyses were performed using the MVMR (version 0.3)¹⁰⁷ and MendelianRandomization (version 0.5.1) packages¹¹⁰. Non-linear MR analyses were performed based on previously described methods¹⁰⁵. For the MR on (subtypes of) mortality, we considered a liberal two-sided P-value of P < 0.05 significant for any of the outcomes using the IVW-MR random effects model. For the MR on cardiovascular diseases, we considered a Bonferroni corrected two-sided P-value for the amount of unique outcomes (P = 0.05/12 = 4.17 × 10⁻³) to be significant for the main IVW-MR random effects analyses, and a two-sided P-value between 4.17 × 10⁻³ and 0.05 to indicate suggestive evidence for an association. A two-sided P-value threshold of P < 0.05 was adopted for the sensitivity analyses.

Meta-subtract of blood pressure traits

We used the MetaSubtract (version 1.60) package in R to remove the effects sizes of the UK Biobank from the largest blood pressure GWAS’s to date in order to obtain the independent effect sizes of the ICBP consortium^19,111. Effect sizes for the UK Biobank were obtained through linear regression analyses using every RHR SNP available in the UK Biobank as exposure, and systolic, diastolic, and pulse pressure as outcomes. Covariates included age, age², sex, BMI, genchip, and PC1-PC30. Cluster-robust standard errors with genetic family IDs as clusters were used to account for the relatedness among participants. Individuals belonging together based on 3rd-degree or closer as indicated by the kinship matrix (kinship coefficient > 0.0442) provided by UK Biobank received a family ID. To keep the cohort similar to the one used in the study from Evangelou et al., we excluded those who self-reported as of non-European ancestry (n = 18,405) and pregnant women (n = 306) from the 484,307 included in the GWAS, leaving 465,659 individuals for the analysis¹⁹. We note that we did not correct the standard errors for the genomic inflation reported for the GWASs of blood pressure traits in the UK Biobank as our linear regression estimates, while resembling the GWAS data, will not be exactly equal to BOLT-LMM estimates due to different methodologies⁸².

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.