A genetic model of ivabradine recapitulates results from randomized clinical trials.

BACKGROUND: Naturally occurring human genetic variants provide a valuable tool to identify drug targets and guide drug prioritization and clinical trial design. Ivabradine is a heart rate lowering drug with protective effects on heart failure despite increasing the risk of atrial fibrillation. In patients with coronary artery disease without heart failure, the drug does not protect against major cardiovascular adverse events prompting questions about the ability of genetics to have predicted those effects. This study evaluates the effect of a variant in HCN4, ivabradine's drug target, on safety and efficacy endpoints.
METHODS: We used genetic association testing and Mendelian randomization to predict the effect of ivabradine and heart rate lowering on cardiovascular outcomes.
RESULTS: Using data from the UK Biobank and large GWAS consortia, we evaluated the effect of a heart rate-reducing genetic variant at the HCN4 locus encoding ivabradine's drug target. These genetic association analyses showed increases in risk for atrial fibrillation (OR 1.09, 95% CI: 1.06-1.13, P = 9.3 ×10-9) in the UK Biobank. In a cause-specific competing risk model to account for the increased risk of atrial fibrillation, the HCN4 variant reduced incident heart failure in participants that did not develop atrial fibrillation (HR 0.90, 95% CI: 0.83-0.98, P = 0.013). In contrast, the same heart rate reducing HCN4 variant did not prevent a composite endpoint of myocardial infarction or cardiovascular death (OR 0.99, 95% CI: 0.93-1.04, P = 0.61).
CONCLUSION: Genetic modelling of ivabradine recapitulates its benefits in heart failure, promotion of atrial fibrillation, and neutral effect on myocardial infarction.

Introduction

Human genetics can be a powerful tool to guide drug development. The identification of mutations in important coronary artery disease associated genes has led to the development of new drugs and the approach of Mendelian randomization (MR) is widely used to predict the effect of interventions on biomarkers [1], to validate drug targets and to predict the effect of drug combinations [2]. There are limitations, however, in the value of human genetics to predict the effects of drugs. The main problems are caused by pleiotropic effects of genetic variants [3], the difference between a lifelong exposure to a risk factor and interventions that are administered after disease onset [2] and the generalizability of the results to specific patient populations and to different ethnic populations [4].

Here, we investigate whether human genetics can reproduce the diverging results obtained on different clinical outcomes in randomized clinical trials of ivabradine. This heart-rate lowering drug was demonstrated to reduce the composite of cardiovascular death and hospitalization for worsening heart failure in patients with symptomatic heart failure and a heart rate above 70 bpm (beats per minute) at baseline in the SHIFT trial [5]. In this study, there was a placebo-adjusted reduction in heart rate of 10.9 (10.4, 11.4) bpm after 28 days on treatment with ivabradine and the hazard ratio (HR) for the cardiovascular composite endpoint was 0.82 (95% CI: 0.75–0.90, p<0.0001). In contrast, in the SIGNIFY trial, ivabradine did not reduce the composite of cardiovascular death or myocardial infarction in patients with stable coronary artery disease (CAD) without heart failure and with a heart rate > 70 bpm at baseline (HR 1.08, 95% confidence interval (CI): 0.96–1.20, P = 0.20) [6]. In both studies, there was an increase in the risk of atrial fibrillation in patients randomized to ivabradine. The incidence of atrial fibrillation was 9% and 8% in the ivabradine and placebo arms in SHIFT (P = 0.012) respectively, whereas it was 5.3% and 3.8% in SIGNIFY (S1 Table). The heart rate reduction induced by ivabradine is due to the inhibition of the “funny” current (If), which is important for cardiac depolarization during phase 4 of the action potential in the sino-atrial node [7]. The hyperpolarization-activated cyclic nucleotide-gated channel 4 encoded by the HCN4 gene is responsible for this current [8]. Here, we study naturally occurring variants at this locus as a genetic model of ivabradine therapy to predict the effects of the drug on heart failure, atrial fibrillation and CAD.

Methods

Data sources

The UK Biobank is a prospective population cohort of over 500,000 individuals aged between 40 and 69 at recruitment, and has been previously described [9]. We used hospitalization data between the beginning of the Health Episode Statistics (HES) linkage (April 1st 1997) and the last available date for the current data release (March 1st 2016). Codes used to define the clinical variables are presented in S2 Table. All UK Biobank participants were previously genotyped. We applied genetic quality control leaving 413,083 individuals for analysis (Methods in S1 Appendix). All reported genomic positions are reported with respect to build GRCh37. We also used the largest available meta-analysis of genome-wide association studies (GWAS) with summary statistics reporting the effect of the HCN4 rs8038766 variant on stroke, atrial fibrillation, CAD, myocardial infarction and heart failure (Methods in S1 Appendix) [10-14]. All participants of the UK Biobank gave their informed consent and the present study was approved by the institutional ethics review board of the Montreal Heart Institute.

Statistical analyses

To identify independent variants at the HCN4 locus (chr15:73,612,200–73,661,605 ± 200kb) associated with resting heart rate at baseline in the UK Biobank dataset, we used forward stepwise linear regression with additive allele coding and a genome-wide significance threshold (p≤5.0 ⨉ 10−8). Association between HCN4 variant rs8038766 and clinical endpoints was assessed using multivariable logistic regression. All models were adjusted for age, sex and the first 10 principal components. For the prospective and competing risk analyses we used Cox proportional-hazards regression. For the competing risk analyses, we estimated the cause-specific hazards where individuals are censored at the time of occurrence of the competing risk if it occurred prior or if it was reported at the same time as the event of interest [15]. We used time from the first baseline visit in years as the timescale and the censure was the date of death or end of follow-up period. For the construction of the heart rate Genetic Risk Score (GRS), we used 64 previously reported genome-wide significant heart rate associated SNPs (with r2 < 0.1) [16]. We split the participants based on the GRS quintiles with the group formed by the 5th GRS quintile (and above) corresponding to the higher heart rate group and the odds ratio for CAD, heart failure and atrial fibrillation were obtained by comparing the first 4 groups individually to the 5th group used as reference in logistic regression. All analyses were performed using the R (v.3.5.2) programming language unless otherwise specified.

Mendelian randomization

We used the inverse variance weighted (IVW) [17], MR-Egger [18], contamination mixture [19] and Mendelian Randomization Pleiotropy RESidual Sum and Outlier (MR-PRESSO) methods [20]. We present results of all four methods in order to outweigh the drawbacks of individual approaches and help guide conclusions, and we report all causal effect estimates for a standard deviation reduction in heart rate as measured in the UK Biobank (11.1 bpm). Coincidentally, this heart rate reduction is similar to the effect of ivabradine in randomized clinical trials (e.g. 10.9 bpm in the SHIFT trial) and is comparable in magnitude to the heart rate reduction by ivabradine. For the heart rate GRS, we used the two-stage method with individual level data and the effect estimates are for a 11.1 bpm reduction in heart rate as well [21]. Analyses were performed with the “MendelianRandomization” R package and MR-PRESSO [20]. Refer to Methods in S1 Appendix for additional details.

Results

Genetic model of ivabradine

To construct a genetic model of ivabradine treatment, we tested the association between variants at the HCN4 locus (defined as the gene boundaries ± 200 kb) and heart rate in the UK Biobank by stepwise forward regression analysis. Two independent signals were identified, the first was led by rs8038766 with every copy of the “G” allele reducing heart rate by 0.57 bpm (95% CI 0.51–0.64, P = 2.76 × 10−66). This variant is in high linkage disequilibrium (LD) with variants previously associated with resting heart rate, heart rate variability traits and atrial fibrillation (S3 Table). The second variant was rs3743496 with the “T” allele reducing heart rate by 0.30 bpm (95% CI 0.25–0.35, P = 3.96 × 10−30). The region spanned by the association signal led by this variant was wide and overlapped more of the neighbouring gene (NEO1), than HCN4 (S1 Fig). Furthermore, the lead variant was in linkage equilibrium with rs8038766 (D’ = 0.43 and χ2 test of independence p-value <0.0001 in 1000 Genomes phase III Europeans) suggesting that the secondary association signal could in fact not be independent of the first. For these reasons, we were not confident that rs3743496 could be used as a specific and independent genetic instrument of HCN4 activity and excluded it from the genetic model of ivabradine. Additionally, HCN4 is a short gene (49,405 bases) that is intolerant to loss-of-function mutations (probability of being loss of function intolerant of 1 in the gnomAD database) [22] which could explain the scarcity of functional variants to be used as genetic instruments.

Genetically predicted effect of ivabradine on safety endpoints

We tested the association between the heart rate reducing allele of the HCN4 variant rs8038766 and atrial fibrillation and stroke, a common and well-known consequence of atrial fibrillation. HCN4 has previously been implicated in atrial fibrillation, and we replicated these results using rs8038766 [23]. In the UK Biobank, rs8038766 was strongly associated with atrial fibrillation (OR 1.09, 95% CI 1.06–1.13; P = 9.3 × 10−9) but not with any stroke or ischemic stroke (Fig 1). The association between rs8038766 and atrial fibrillation was also observed in summary statistics from previously published GWAS of atrial fibrillation with OR = 1.11 (P = 1.8 × 10−26), and 1.12 (P = 5.4 × 10−35) for Roselli et al. and Nielsen et al. respectively (Fig 1) [12, 13]. Previous epidemiologic studies have shown that chronic atrial fibrillation leads to a five-fold increase in the risk of stroke [24]. We did not find a significant association between rs8038766 and stroke in the UK Biobank, potentially because of the low number of cases (4,158 cases for ischemic stroke). Summary results from the MEGASTROKE consortium show an association between rs8038766 and cardioembolic (OR = 1.08, 95% CI 1.03–1.13, P = 1.54 × 10−3) and ischemic stroke (OR = 1.03, 95% CI 1.01–1.05, P = 0.0152) (Fig 1) [25].

Fig 1

Association between the heart rate lowering allele (G) of the HCN4 variant rs8038766 and safety outcomes in the UK Biobank and in published GWAS summary statistics from large consortia.

For the UK Biobank, reporting results from logistic regression comparing the combined prevalent and incident cases to non-cases. References: Roselli et al. [13], Nielsen et al. [12], MEGASTROKE [25]. * rs7174098 (LD r2 = 1 in 1000 genomes Europeans) was used instead of rs8038766 as the latter was unavailable in the MEGASTROKE summary statistics for this outcome.

Genetically predicted effect of ivabradine on efficacy endpoints

We tested for association of the heart rate-reducing allele at the HCN4 variant rs8038766 with combined prevalent and incident heart failure in the UK Biobank and found no association in this naïve model (OR = 0.96, 95% CI 0.91–1.00, p = 0.071) (Fig 2). However, because atrial fibrillation is an important risk factor for heart failure [26], it is a possible that the increased risk of atrial fibrillation attenuates a possible protective association with heart failure. Indeed, after adjustment for any prevalent or incident atrial fibrillation, the association of rs8038766 with heart failure was OR = 0.91, 95% CI 0.87–0.96 (P = 6.7 × 10−4). In a model including the interaction term between rs8038766 and atrial fibrillation, the estimated coefficient for the variant was -0.136, 95% CI -0.205, -0.067, (P = 0.00011) and the interaction term coefficient was 0.110 95% CI 0.005, 0.215 (P = 0.04). These coefficients correspond to an estimated OR of the SNP on heart failure of 0.87 in individuals without atrial fibrillation and 0.97 in individuals with atrial fibrillation. However, these associations could be biased if both the SNP and the outcome increase atrial fibrillation risk resulting in a possible collider bias. To account for this, we used a cause-specific hazards model for the incidence of heart failure and atrial fibrillation separately using 404,767 UK Biobank participants that were free of both diseases at baseline. In this group, there were 3,385 incident heart failure cases and the HCN4 variant rs8038766 showed a non-significant trend for a protective effect (HR = 0.96, 95% CI: 0.89–1.02; P = 0.177) (Table 1). However, in a competing risk model accounting for incident occurrences of atrial fibrillation, the protective effect of the heart rate-reducing variant on heart failure was brought to focus with HR = 0.90, 95% CI 0.83–0.98 (P = 0.013) (Table 1). We conducted a similar analysis using incident myocardial infarction or cardiovascular death corresponding to the primary endpoint in the SIGNIFY trial, which was also potentially exposed to the opposing effects of the heart rate-reducing variant on atrial fibrillation and myocardial infarction. There was no detectable association of the heart rate-reducing variant with myocardial infarction or cardiovascular death in the simple Cox proportional-hazards model (HR = 0.99, 95% CI 0.94–1.05) or in the cause-specific competing risk model (HR = 0.99, 95% CI 0.93–1.04) (Table 1). We did see, however, an association between rs8038766 and prevalent or incident cases of unstable angina (OR 0.92 95% CI 0.86–0.98, P = 0.0056) (Fig 2).

Fig 2

Association between the heart rate lowering allele (G) of the HCN4 variant rs8038766 and efficacy outcomes in the UK Biobank and in GWAS summary statistics from large consortia.

For the UK Biobank, reporting results from logistic regression comparing the combined prevalent and incident cases to non-cases. References: CARDIoGRAMplusC4D [10], UKB SOFT CAD + CARDIoGRAMplusC4D + MI Genetics [11], HERMES [14].

Table 1

Association of the HCN4 variant with outcomes in the UK Biobank using prospective and cause-specific hazard competing risk analyses.

Outcome	Model	N total	N events	HR (95% CI) *	P-value
Genetic model for SHIFT**
Using participants without a history of atrial fibrillation or heart failure at recruitment
Heart failure	Cox proportional-hazards	404,767	3,385	0.96 (0.89, 1.02)	0.18
Atrial fibrillation	Cox proportional-hazards	404,767	8,461	1.08 (1.04, 1.13)	9.4 × 10−5
Heart failure	Competing risk (atrial fibrillation)	404,767	2,380	0.90 (0.83, 0.98)	0.013
Atrial fibrillation	Competing risk (heart failure)	404,767	7,663	1.08 (1.04, 1.13)	3.2 × 10−4
Genetic model for SIGNIFY**
Using participants without a history of atrial fibrillation or MI at recruitment
MI or CV Death	Cox proportional-hazards	397,008	4,976	0.99 (0.94, 1.05)	0.84
Atrial fibrillation	Cox proportional-hazards	397,008	7,880	1.08 (1.04, 1.13)	3.1 × 10−4
MI or CV Death	Competing risk (atrial fibrillation)	397,008	4,534	0.99 (0.93, 1.04)	0.61
Atrial fibrillation	Competing risk (MI or CV death)	397,008	7,482	1.09 (1.04, 1.13)	1.4 × 10−4

* Reporting the effect of the heart rate reducing allele of rs8038766 at the HCN4 gene. All models were adjusted for age, sex and the first 10 principal components. In the Cox proportional-hazards models, individuals were censored at the time of death or end of follow up; in the competing risk models, individuals were censored at the time of occurrence of the competing event, death, or end of follow up.

** Our model aims to match the outcome and exposure of interest from the SHIFT and SIGNIFY trials, but we did not emulate the trials in any other way such as by matching the inclusion / exclusion criteria.

CV, cardiovascular; HR, hazard ratio; MI, myocardial infarction.

In the HERMES case-control consortium, the heart rate reducing allele of rs8038766 was only weakly associated with heart failure (OR 0.98 95% CI 0.96–1.00, p = 0.079) (Fig 2), but when using the mtCOJO method to adjust for atrial fibrillation using summary statistics [27], the protective effect was increased with a conditional OR = 0.96 95% CI 0.94–0.98 (P = 9.7×10−4) [14]. In the CARDIoGRAMplusC4D consortium, there was no association between rs8038766 and CAD or myocardial infarction (Fig 2).

Bi-directional MR

Bi-directional MR supports a causal effect of atrial fibrillation on heart failure with a causal OR estimate of 1.22 and ranging up to 1.25 according to different MR models (Table 2, S7 Table), and supports a causal effect of heart failure on atrial fibrillation with OR ranging from 1.21–1.94 (excluding the contamination mixture model estimate of 6.82 which is an outlier among the other methods). These results are concordant with observational longitudinal studies that have observed an increased incidence of atrial fibrillation in new heart failure patients and vice versa and where both diseases are often diagnosed on the same day [26]. Bi-directional MR also supported a causal effect of CAD and myocardial infarction on atrial fibrillation, but not the opposite. The point estimates ranged from OR: 1.12–1.17 for the effect of CAD on atrial fibrillation and OR: 1.11–1.22 for the effect of myocardial infarction on atrial fibrillation (Table 2, S7 Table).

Table 2

Bi-directional Mendelian randomization estimates.

Exposure	Outcome	MR Causal OR (95% CI) *	P-value
Atrial fibrillation (152 variants)	Heart failure	1.23 (1.20, 1.27)	3.7 × 10−52
Atrial fibrillation (152 variants)	Coronary artery disease	1.00 (0.98, 1.03)	0.76
Atrial fibrillation (152 variants)	Myocardial infarction	0.98 (0.95, 1.02)	0.30
Heart failure (11 variants)	Atrial Fibrillation	1.45 (1.11, 1.90)	0.0067
Coronary artery disease (68 variants)	Atrial Fibrillation	1.15 (1.11, 1.21)	1.7 × 10−10
Myocardial infarction (31 variants)	Atrial Fibrillation	1.11 (1.06, 1.16)	1.3 × 10−5

Summary statistics for atrial fibrillation taken from Nielsen et al. [12] for myocardial infarction and CAD from CARDIoGRAMplusC4D and CARDIoGRAMplusC4D + UKB SOFT + MiGen [10, 11], for heart failure from HERMES [14].

* IVW MR model. For MR results using MR-Egger, the contamination mixture model and MR-PRESSO, see S7 Table. Causal ORs relate the odds of the outcome in exposed individuals vs non-exposed.

IVW, inverse-variance weighted; MR, Mendelian randomization; OR, odds ratio.

Effect of heart rate on cardiovascular outcomes

To compare our observations of heart rate reduction attributable to the HCN4 variant to that of polygenic origin, we constructed a genetic risk score (GRS) with 64 variants previously associated with heart rate (S4 Table) [16]. An increase of 1 standard deviation in the heart rate GRS was associated with a 1.76 (1.73, 1.80) bpm increase in heart rate explaining 2.5% of the variance in the UK Biobank data (S2 Fig). The two stage method causal estimate scaled for a 11.1 bpm (corresponding to 1 standard deviation of heart rate in the UK Biobank and concordant with the effect of ivabradine in clinical trials) genetic reduction in heart rate was OR = 1.25 (95% CI: 1.13–1.39) for atrial fibrillation, OR = 1.03 (95% CI: 0.88–1.21) for heart failure and 1.03 (95% CI: 0.96–1.11) for CAD. To account for the presence of pleiotropy, we also used MR-Egger, contamination mixture model and the MR-PRESSO methods (S5 Table), and saw an increase in the risk for atrial fibrillation associated with heart rate reduction (OR 1.54, P = 1.3 × 10−7 for MR-PRESSO), but not with heart failure or CAD, although these did not take into account the possible competing effect of atrial fibrillation. MR analyses using effect estimates derived from larger GWAS consortia using the same set of 64 heart rate variants supported results observed with the UK Biobank data (S6 Table).

Discussion

In the present study, we used genetics to infer the causal effect of ivabradine on safety and efficacy outcomes in an attempt to reproduce observations from randomized clinical trials, and to assess the value of genetic approaches to support drug targets and trial design issues such as target patient population and clinical outcomes.

Effect of HCN4 on atrial fibrillation

Genetically predicted heart rate reduction from the HCN4 gene variant rs8038766 was associated with an increase in risk of atrial fibrillation, recapitulating the observations from the SIGNIFY and SHIFT trials. In our MR analyses using methods robust to the inclusion of invalid instruments, we observed that a genetically predicted reduction in heart rate of approximately 11 bpm conferred an increased risk of atrial fibrillation with an OR of 1.54 in the UK Biobank and OR ranging from 1.36 to 1.56 using summary statistics from previous large GWAS. These results are also coherent with a recent MR study reporting a protective effect of increased heart rate for atrial fibrillation and cardioembolic stroke [28]. Atrial fibrillation is known to increase the risk of stroke by 3 to 5-fold [24]. In clinical trials of ivabradine, there was no treatment association with stroke, but the small number of incident atrial fibrillation events would have made such an observation unlikely [29]. The estimated effects of heart rate on atrial fibrillation are smaller than the effect predicted from the HCN4 variant alone, whose scaled OR for a comparable 11 bpm reduction would be greater than 5. This may be explained partly by the inaccuracy of extrapolation of the OR estimate derived from a single genetic variant, and also possibly by an effect of HCN4 on atrial fibrillation that may be specific to modulation of the If current or other structural consequences of HCN4 variants. For example, genetic mutations in HCN4 have been associated to Brugada syndrome and sick sinus syndrome as well as left ventricular noncompaction and it is possible that common polymorphisms in the HCN4 gene have more subtle effects on myocardium structure or conduction parameters that may be independent of heart rate [30, 31]. Additionally, altering HCN4 function or levels during embryogenesis in other species have been shown to structurally alter heart development which could explain effects beyond heart rate modulation alone [32]. The MR estimates from the UK Biobank are also based on mostly healthy individuals with a low heart rate (mean of 69 bpm) possibly limiting clinical interpretation [33]. However, the effect estimates of heart rate reduction on atrial fibrillation are similar when using the GWAS results from Nielsen et al. which include both population-based and clinical cohorts [12]. The possibility that the effect is greater in healthy patients is also supported by the previously described association between low heart rate during physical activity and the increased incidence of atrial fibrillation [34].

Effect of HCN4 on ischemic endpoints

We tested the association between the HCN4 heart rate-reducing variant and various ischemic endpoints in the UK Biobank. The largest effect we observed was with unstable angina, which is coherent with the use of ivabradine to alleviate anginal symptoms. Nonetheless, the effect sizes of the association with CAD, angina and myocardial infarction were small and marginally significant in the UK Biobank and importantly they were not supported by results from larger GWAS consortia. This suggests that the effect of HCN4 on CAD may be null or of a very small effect size so as to not be detectable in the context of a clinical trial such as in the SIGNIFY study [6]. We also investigated the possibility that the increased risk of atrial fibrillation offsets the beneficial effects on the SIGNIFY primary endpoint of myocardial infarction or cardiovascular death using a prospective competing risk analysis in individuals that did not develop atrial fibrillation and showed that accounting for atrial fibrillation had no impact on the risk for myocardial infarction or cardiovascular death. There was no detectable association of the HCN4 heart rate-reducing variant with myocardial infarction or cardiovascular death in the cause-specific competing risk model. This was further supported by the bi-directional MR analysis that showed that CAD caused atrial fibrillation but not the opposite. Finally, the MR study did not show a causal link between heart rate and CAD suggesting that reducing heart rate is not sufficient to prevent the disease.

Relationship with clinical trials of ivabradine

The analysis of the subgroup of participants with angina class 2 or greater at baseline in SIGNIFY showed a nominal increase in the rate of the primary endpoint of cardiovascular death or myocardial infarction with ivabradine [5, 6]. Whether this observation represented a chance finding in the context of a neutral result in the overall SIGNIFY population or a potential signal of harm in this subset of patients was a matter of discussion. The results of the current analyses support neutral effects of If current inhibition on the composite endpoint of cardiovascular death and myocardial infarction, without evidence of harm.

In the SHIFT trial, ivabradine reduced the rate of the primary composite endpoint of cardiovascular death or hospitalization for worsening heart failure in patients with heart failure with reduced ejection fraction and without atrial fibrillation. In the genetic model of ivabradine, the competing risk analysis accounting for atrial fibrillation showed that the HCN4 heart rate-reducing variant protected against heart failure (HR = 0.90, 95% CI: 0.83–0.98, P = 0.013). The results from the marginal models and the competing risk analyses do suggest opposing effects of the heart rate-reducing HCN4 variant on atrial fibrillation and heart failure. The importance of these effects is also highlighted by the bi-directional MR of atrial fibrillation and heart failure that confirmed that both diseases are mutually causal of one another.

Study limitations

As for any MR study, our analyses were subject to the assumptions of the underlying models and the possibility of unobserved horizontal pleiotropy. Additionally, our genetic model of ivabradine corresponds to a lifelong effect as opposed to an exposure after drug initiation. Generally, common variants also result in effects of smaller magnitude than ones resulting from pharmacological modulation and extrapolation is required to compare them. We also used data from individuals of predominantly European ancestry both in the UK Biobank and in summary statistics from large GWAS consortia which could limit the generalizability of our results to other populations both in terms of clinical profile and ancestry. Finally, we defined clinical variables based on combinations of hospitalization and death record codes in the UK Biobank which is likely to result in imperfect coding of disease status.

Conclusion

In conclusion, genetic modelling of ivabradine recapitulates its benefits in heart failure, promotion of atrial fibrillation, and neutral effect on myocardial infarction. This study supports the use of methods that leverage naturally occurring genetic variants to predict diverging results on different clinical outcomes and support the design of randomized clinical trials, even in a situation where more complex disease risks are at play.

Supplementary methods and references for supplementary methods, figures and tables.

(DOCX)

Click here for additional data file.

Results from the stepwise regression of HCN4 variants on heart rate in the UK Biobank.

(DOCX)

Click here for additional data file.

Effect of heart rate genetic risk score groups based on quintiles on atrial fibrillation, heart failure and coronary artery disease in the UK biobank dataset.

(DOCX)

Click here for additional data file.

Summary of ivabradine cardiovascular outcomes trials.

(DOCX)

Click here for additional data file.

Self-reported, hospitalization (ICD10) and operation (OPCS) codes used to define clinical variables based on the UK Biobank available data.

(DOCX)

Click here for additional data file.

Results from the NHGRI-EBI GWAS catalog mapped to the HCN4 gene.

(DOCX)

Click here for additional data file.

Variants and weights used for the computation of the heart rate GRS.

(DOCX)

Click here for additional data file.

MR estimates based on 64 heart-rate associated variants and their effect on outcomes in the UK Biobank.

(DOCX)

Click here for additional data file.

MR estimates based on the effect of 64 heart-rate associated variants in external summary statistics from large GWAS consortia.

(DOCX)

Click here for additional data file.

Participant overlap between GWAS meta-analysis studies used for observational and Mendelian randomization analyses and the UK Biobank.

(DOCX)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

13 May 2020

PONE-D-20-09954

A genetic model of ivabradine recapitulates results from randomized clinical trials

PLOS ONE

Dear Dr. Dubé,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

We would appreciate receiving your revised manuscript by Jun 27 2020 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Ify Mordi

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for stating the following in the Competing Interests section:

I have read the journal's policy and the authors of this manuscript have the following competing interests: A patent pertaining to pharmacogenomics-guided CETP inhibition was granted and J.-C.Tardif and M.-P. Dubé are mentioned as authors. J.-C.Tardif and M.-P. Dubé have a minor equity interest in DalCor. J.-C.Tardif has received research support from Amarin, AstraZeneca, DalCor, Eli-Lilly, Ionis, Pfizer, RegenexBio, Sanofi and Servier, and honoraria from DalCor, Pfizer, Sanofi and Servier. M.-P. Dubé has received honoraria from Dalcor and research support (access to samples and data) from AstraZeneca, Pfizer, Servier, Sanofi and GlaxoSmithKline. B. Tyl is an employee of Laboratoires Servier. Simon de Denus has received grants from Pfizer, AstraZeneca, Roche Molecular Science, DalCor and Novartis. R. Thomas Lumbers has received research grants from Pfizer. The remaining authors have nothing to disclose.

We note that one or more of the authors are employed by a commercial company: Laboratoires Servier

a. Please provide an amended Funding Statement declaring this commercial affiliation, as well as a statement regarding the Role of Funders in your study. If the funding organization did not play a role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript and only provided financial support in the form of authors' salaries and/or research materials, please review your statements relating to the author contributions, and ensure you have specifically and accurately indicated the role(s) that these authors had in your study. You can update author roles in the Author Contributions section of the online submission form.

Please also include the following statement within your amended Funding Statement.

“The funder provided support in the form of salaries for authors [insert relevant initials], but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the ‘author contributions’ section.”

If your commercial affiliation did play a role in your study, please state and explain this role within your updated Funding Statement.

b. Please also provide an updated Competing Interests Statement declaring this commercial affiliation along with any other relevant declarations relating to employment, consultancy, patents, products in development, or marketed products, etc.

We also note that you have a patent relating to material pertinent to this article. Please ensure your amended statement of Competing Interests declares this patent (with details including name and number), along with any other relevant declarations relating to employment, consultancy, patents, products in development or modified products etc.

Within your Competing Interests Statement, please confirm that this commercial affiliation and patent does not alter your adherence to all PLOS ONE policies on sharing data and materials by including the following statement: "This does not alter our adherence to  PLOS ONE policies on sharing data and materials.” (as detailed online in our guide for authors http://journals.plos.org/plosone/s/competing-interests) . If this adherence statement is not accurate and  there are restrictions on sharing of data and/or materials, please state these. Please note that we cannot proceed with consideration of your article until this information has been declared.

c. Please include both an updated Funding Statement and Competing Interests Statement in your cover letter. We will change the online submission form on your behalf.

4. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors studied genetic variants in the HCN4 gene for association with heart failure, atrial fibrillation and coronary artery disease. They reconfirm that genetic variants in this region, which primarily associate with heart rate, have also significant association with atrial fibrillation and mildly with cardioembolic stroke. There were also some signals for association with heart failure, particularly after adjustment for atrial fibrillation. Moreover, the authors studied a genetic score for heart rate and cardiovascular outcome and confirm that a reduction in heart rate is positively associated with atrial fibrillation but no other obvious cardiovascular conditions. The manuscript largely recapitulates findings that have been made for genetic variants at the HCN4 locus published in genomewide association studies. The strength of the paper is the focus on predicting Ivabradine related effects (i.e. pharmacological effects similar to those of genetic variants) as observed in previous clinical trials.

Minor comments

It would be interesting to see whether the genetic variants tested affect expression of the HCN4 gene as can be studied in publically available databases.

In the context of the reported variants, the authors should not use the term mutation.

Reviewer #2: Sorry, I don't have a lot of time to re-read. I reviewed this manuscript at a previous journal and liked it, but it wasn't prioritized for publication there. I've copy-pasted my comments from the previous review - if any comments have already been addressed in the editing process between journals, please indicate that the comment is no longer relevant. Thanks and best wishes, Steve

---

Please note that I am a statistician, so my comments relate to the statistical aspects of the analysis - I am not best placed to comment on other aspects.

Generally speaking, this was a persuasive and thorough exposition of the interrelations between various cardiovascular diseases. I have only minor comments:

1) I'd appreciate if someone with more knowledge than me was able to comment on the use of the rs8038766 variant only, and the fact that the rs3743496 variant was ignored. Having a high LD-score isn't necessarily a bad thing (depends on what its neighbours are), and you would expect a variant that is selected in a region-wide search to have a higher LD score than average (and 87th %ile isn't especially high!). The argument that the variant overlaps another gene region is a more persuasive reason for excluding it from the analysis, but many variants have LD regions that overlap different genes. Anyway, I trust that this was a decision made for principled reasons rather than for convenience (the results looked better when looking at one variant only), but the reasons given seem to me a little thin.

2) Estimates in the tables in units. In Table 2, what are the units for the causal ORs? (are they per unit increase in the log odds of the exposure?). In Figures 1/2, could write "Odds ratio per heart rate lowering allele" on the figure (this info is in the legend). Is there any systematic difference between these figures? In the Take-home figure, I'm not fully clear what "Genetic model of SHIFT" means. I presume 11 bpm is 1 SD for heart rate, but it could be the mean effect of taking ivabradine.

3) I'd appreciate if someone with more knowledge than me was able to comment on the plausibility of the bidirectional analyses. Is it plausible to suggest that there could be bidirectional effects of AF on heart failure? Does AF sometimes precede HF, and sometimes HF precedes AF? If not, this doesn't invalidate the analyses (could be that the genetic variants are picking up an effect of subclinical HF on AF risk or vice versa), but it changes the interpretation - in my (limited) experience, true bidirectional effects are rare.

4) This is really picky, but generally the term "2SLS method" implies a continuous outcome and a linear regression model. If you used linear regression and treated the binary outcome as continuous, then fair enough (it doesn't make much difference). But if you used logistic/Cox PH regression, then "two-stage method" is my preferred term (others have been suggested).

Otherwise, I don't have much else to say - it was an interesting read!

---

Stephen Burgess

Reviewer #3: The authors present their findings on using genetic variation in HCN4, the gene which encodes the target of ivabradine, to replicate the findings of previously published drug trials of ivabradine. Major strengths are the various complementary analyses using large-scale data sources and the detailed documentation. Overall the manuscript is well-written. Some observations and questions:

- It seems that the UK Biobank was a large contributor to many of the sources of summary statistics. It would be good to provide the reader insight in the % sample overlap (always with respect to the larger study) for the various data sources which are combined across the different analyses.

- Given that the genetic risk score for heart rate was derived from a GWAS meta-analysis where the UK Biobank formed the discovery stage, any MR analyses performed in the one-sample setting of the UK Biobank with this score might suffer from the Winner’s Curse. How would this influence the reported results?

- On page 10 the authors describe how adjusting for atrial fibrillation would be problematic if both the SNP and heart failure collide on AF, which would introduce collider bias. However, isn’t it equally likely and problematic that, if the SNP has an effect on AF, and AF has an effect on heart failure, that both the SNP and the confounders of the AF-heart failure association would collide on AF (leading to collider bias when you adjust for AF)?

- The interpretation of main effect estimates are less straightforward when interaction effects have been added to the model. Therefore, please be explicit how the reader should interpret the sentence providing both the main and interaction (with AF) estimates between rs8038766 and heart failure.

- Figure 1: Why was rs7174098 used for just one outcome in MEGASTROKE?

- Supposedly you choose the transethnic data of MEGASTROKE for its large number of cases. However, rs8038766 need not necessarily be a strong genetic proxy for HCN4 in non-Europeans. Does using METASTROKE's European dataset show comparable results?

- Please report (perhaps in supplemental material) whether the various GWAS meta-analyses were based on incident or prevalent cases and whether recurrent events were included.

- There exist MR methods which can incorporate correlated genetic variants to boost power. Did you consider these methods for HCN4-variants?

- Table 1 describes the ‘genetic model for SHIFT/SIGNIFY’. Please be explicit this only refers to the outcome definition (and intervention), i.e., not also the inclusion criteria for participants.

- The selected population of the UK Biobank may give rise to issues like selection bias, also for MR studies. Could this have influenced your results?

Minor:

- Please mention that the heart rate GRS includes variants which are just relatively independent (r2<0.1). Lower r2 thresholds are now typically advised for MR (e.g., 0.001). In extension, please be more explicit regarding the independence of the various sets of instruments used in the bidirectional MR analyses.

- Perhaps of interest, FINNGEN could serve as an additional source of publicly available summary statistics for (ICD-code based) heart failure

- For the analyses with the GRS it seems two-sample methodology was applied in the one-sample MR setting (e.g., Supplemental Table 5). Perhaps of interest: this (very) recent preprint (https://www.biorxiv.org/content/10.1101/2020.05.07.082206v1) suggests that particularly the MR-Egger method may be easily biased in this setting. Calculating the I2 may give insight.

- Not a fan of ‘non-significant trend’ – perhaps rephrase to ‘provided (very) weak evidence’?

- With regard to the kinship threshold of >0.0884 – wouldn’t this reflect 2nd degree relationships or more (rather than ‘or less’)?

- Please note that quintiles are not the group themselves but rather the cut-offs to define these groups

- The abbreviation GRS is introduced twice in the methods section

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Stephen Burgess

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

17 Jun 2020

Editor comment 1: Style requirements

We have now modified the manuscript as needed to comply with the PLOS One style requirements.

Editor comment 2: competing interests and funding statements

We have updated the competing interests and funding statements as required to highlight the possible role of funders in our study. The new statements were enclosed within the cover letter. Note that the patent on pharmacogenomics-guided CETP inhibition is only distantly related to the topic of the manuscript as it concerns a drug of a different class (CETP-inhibition) and used for a different indication than ivabradine. We will clarify the information in the individual declarations of Competing Interests.

The patent title is “Methods for Treating or Preventing Cardiovascular Disorders and Lowering Risk of Cardiovascular Events”, the patent number is US20190070178A1, the patent is owned by Dalcor and authors receive no royalties.

Editor comment 4: Supporting information

We have added captions for the S1 Appendix, S1-S2 Figures and S1-S7 Tables at the end of the manuscript. We also updated the in-text references.

Reviewer’s comments to the Authors :

Reviewer #1:

The authors studied genetic variants in the HCN4 gene for association with heart failure, atrial fibrillation and coronary artery disease. They reconfirm that genetic variants in this region, which primarily associate with heart rate, have also significant association with atrial fibrillation and mildly with cardioembolic stroke. There were also some signals for association with heart failure, particularly after adjustment for atrial fibrillation. Moreover, the authors studied a genetic score for heart rate and cardiovascular outcome and confirm that a reduction in heart rate is positively associated with atrial fibrillation but no other obvious cardiovascular conditions. The manuscript largely recapitulates findings that have been made for genetic variants at the HCN4 locus published in genomewide association studies. The strength of the paper is the focus on predicting Ivabradine related effects (i.e. pharmacological effects similar to those of genetic variants) as observed in previous clinical trials.

Author response: Thank you for your comments and agreeing to review our manuscript.

Minor comments

Comment #1: It would be interesting to see whether the genetic variants tested affect expression of the HCN4 gene as can be studied in publically available databases.

Author response: Thank you for this suggestion. We agree that using eQTLs of HCN4 would have been ideal, however the tissue of interest in this study is the sinoatrial node which is not available in gene expression resources. As such, the identification of strong eQTLs valuable for MR in our study is not possible.

Comment #2: In the context of the reported variants, the authors should not use the term mutation.

Author response: As suggested by the reviewer, we have adapted the wording where we refer to common HCN4 polymorphisms. Specifically, we replaced “mutation” by “variant” when referring to common variants everywhere in the manuscript.

Reviewer #2:

Sorry, I don't have a lot of time to re-read. I reviewed this manuscript at a previous journal and liked it, but it wasn't prioritized for publication there. I've copy-pasted my comments from the previous review - if any comments have already been addressed in the editing process between journals, please indicate that the comment is no longer relevant. Thanks and best wishes, Steve

Please note that I am a statistician, so my comments relate to the statistical aspects of the analysis - I am not best placed to comment on other aspects.

Generally speaking, this was a persuasive and thorough exposition of the interrelations between various cardiovascular diseases. I have only minor comments:

Author response: Thank you for your comments and accepting to review our manuscript. We had included many of your previous recommendations into the current version of the manuscript. We will highlight the changes in response to the specific points below.

Comment #1: I'd appreciate if someone with more knowledge than me was able to comment on the use of the rs8038766 variant only, and the fact that the rs3743496 variant was ignored. Having a high LD-score isn't necessarily a bad thing (depends on what its neighbours are), and you would expect a variant that is selected in a region-wide search to have a higher LD score than average (and 87th %ile isn't especially high!). The argument that the variant overlaps another gene region is a more persuasive reason for excluding it from the analysis, but many variants have LD regions that overlap different genes. Anyway, I trust that this was a decision made for principled reasons rather than for convenience (the results looked better when looking at one variant only), but the reasons given seem to me a little thin.

Author response: We agree with the reviewer’s assessment that the selection of a single variant as genetic instrument is not ideal. This is our only course of action, unfortunately, as HCN4 is a short and highly constrained gene with very few variants known to affect function (pLI = 1 on gnomAD v2.1.1).

We have opted to not use the second variant identified in the stepwise analysis, rs3743496, because it spans an association signal that overlaps more with the NEO1 gene than the HCN4 gene (Online Figure 1). Furthermore, the Online Figure 1 suggests residual LD between both signals questioning the true independence of both signals. In the 1000 genomes phase III Europeans, the D’ between both SNPs was 0.43 suggesting a strong deviation from independence even though the r2 was 0.02. This discrepancy between metrics of linkage disequilibrium can be explained by the large difference in MAF (0.16 for rs8038766 and 0.38 for rs3743496). In light of these findings, we believe that the use of rs8038766 alone is likely to be a good proxy for the true HCN4 modulating variants and that the inclusion of the second variant would not lead to added benefit. Moreover, there are more independently identified GWAS associations with heart rate parameters and heart rate variability traits in SNPs in LD with rs8038766 than with rs3743496 (Supplementary Table 3) suggesting that it is a better tag for HCN4 functional variants.

We have adapted the text to add these justifications. The “results / genetic model of ivabradine” section now reads:

“This variant is in high linkage disequilibrium (LD) with variants previously associated with resting heart rate, heart rate variability traits and atrial fibrillation (Supplementary Table 3). The second variant was rs3743496 with the “T” allele reducing heart rate by 0.30 bpm (95% CI 0.25-0.35, P=3.96 × 10-30). The region spanned by the association signal led by this variant was wide and overlapped more of the neighbouring gene (NEO1), than HCN4 (Online Figure 1). Furthermore, the lead variant was in linkage equilibrium with rs8038766 (D’ = 0.43 and �2 test of independence p-value <0.0001 in 1000 Genomes phase III Europeans) suggesting that the secondary association signal could in fact not be independent of the first. For these reasons, we were not confident that rs3743496 could be used as a specific and independent genetic instrument of HCN4 activity and excluded it from the genetic model of ivabradine. Additionally, HCN4 is a short gene (49,405 bases) that is intolerant to loss-of-function mutations (probability of being loss of function intolerant of 1 in the gnomAD database) 22 which could explain the scarcity of functional variants to be used as genetic instruments.”

Comment #2: Estimates in the tables in units. In Table 2, what are the units for the causal ORs? (are they per unit increase in the log odds of the exposure?).

Author response: We have added that the IVW OR estimates from Table 2 are to be interpreted as the increase in odds of the outcome in exposed vs non-exposed.

In Figures 1/2, could write "Odds ratio per heart rate lowering allele" on the figure (this info is in the legend). Is there any systematic difference between these figures?

Author response: We have adapted Figures 1 and 2 to read “Odds ratio per heart rate lowering allele” as recommended.

Figures 1 and 2 are very similar but Figure 1 shows results for safety outcomes whereas Figure 2 shows results for efficacy outcomes (with respect to ivabradine treatment). Apart from that distinction there are no systematic differences.

In the Take-home figure, I'm not fully clear what "Genetic model of SHIFT" means. I presume 11 bpm is 1 SD for heart rate, but it could be the mean effect of taking ivabradine.

Author response: We have removed the “take home” figure for this journal as it was not required. To clarify, the 11 bpm does correspond to 1 s.d. of heart rate in the UK Biobank, but the effect of taking ivabradine in RCTs was comparable in magnitude (e.g. placebo adjusted net reduction of 10.9 bpm in ivabradine arm of the SHIFT trial at 28 days). To clarify this important point, we have added the following text in the “Mendelian randomization” subsection of the “Methods” section:

“We present results of all four methods in order to outweigh the drawbacks of individual approaches and help guide conclusions, and we report all causal effect estimates for a standard deviation reduction in heart rate as measured in the UK Biobank (11.1 bpm). Coincidentally, this heart rate reduction is similar to the effect of ivabradine in randomized clinical trials (e.g. 10.9 bpm in the SHIFT trial) and is comparable in magnitude to the heart rate reduction by ivabradine.”

Comment #3: I'd appreciate if someone with more knowledge than me was able to comment on the plausibility of the bidirectional analyses. Is it plausible to suggest that there could be bidirectional effects of AF on heart failure? Does AF sometimes precede HF, and sometimes HF precedes AF? If not, this doesn't invalidate the analyses (could be that the genetic variants are picking up an effect of subclinical HF on AF risk or vice versa), but it changes the interpretation - in my (limited) experience, true bidirectional effects are rare.

Author response: Thank you for this remark. In general clinicians are not surprised by this bi-directional relationship as it is consistent with clinical observations. More rigorously, these results are consistent with observations from the Framingham longitudinal study as mentioned in the “Results / Bi-directional MR” section:

“These results are concordant with observational longitudinal studies that have observed an increased incidence of atrial fibrillation in new heart failure patients and vice versa and where both diseases are often diagnosed on the same day 26.”

Furthermore, the shared risk factors and etiology of AF and HF is well established in the literature. For a recent review, see:

Carlisle, Matthew A., Marat Fudim, Adam D. DeVore, and Jonathan P. Piccini. 2019. “Heart Failure and Atrial Fibrillation, Like Fire and Fury.” JACC. Heart Failure 7 (6): 447–56.

Comment #4: This is really picky, but generally the term "2SLS method" implies a continuous outcome and a linear regression model. If you used linear regression and treated the binary outcome as continuous, then fair enough (it doesn't make much difference). But if you used logistic/Cox PH regression, then "two-stage method" is my preferred term (others have been suggested).

Author response: We agree with this comment and the importance of avoiding misleading wording of statistical methods. We have adapted the text to use two-stage method when the 2nd stage involves a logistic regression. Concretely, the text in “Methods / Mendelian randomization” now reads:

“For the heart rate Genetic Risk Score (GRS), we used the two-stage method with individual level data and the effect estimates are for a 11.1 bpm reduction in heart rate as well”

Otherwise, I don't have much else to say - it was an interesting read!

---

Stephen Burgess

Reviewer #3:

The authors present their findings on using genetic variation in HCN4, the gene which encodes the target of ivabradine, to replicate the findings of previously published drug trials of ivabradine. Major strengths are the various complementary analyses using large-scale data sources and the detailed documentation. Overall the manuscript is well-written. Some observations and questions:

Author response: We thank the reviewer for the comments on the manuscript.

Comment #1: It seems that the UK Biobank was a large contributor to many of the sources of summary statistics. It would be good to provide the reader insight in the % sample overlap (always with respect to the larger study) for the various data sources which are combined across the different analyses.

Author response: We now provide in Supplementary Table 7 an estimate on the overlap of UK Biobank samples from the large consortiums based on the published numbers of cases and controls.

The “Supplementary Methods / External summary statistics” section now reads:

“To provide insight into the overlap between samples from the external summary statistics and the UK Biobank, we summarized the shared number of cases and controls in Supplementary Table 7.”

And the table shows an overlap in cases ranging from 13.7% (HERMES) to 24.5% (Roselli et al.) and an overlap in controls ranging from 39.3% (Nielsen et al.) to 64.1% (Roselli et al.).

Comment #2: Given that the genetic risk score for heart rate was derived from a GWAS meta-analysis where the UK Biobank formed the discovery stage, any MR analyses performed in the one-sample setting of the UK Biobank with this score might suffer from the Winner’s Curse. How would this influence the reported results?

Author response: We thank the reviewer for highlighting this important consideration. In our analyses based on the heart rate score, we have used 5,000 bootstrap resamples to empirically estimate the variance of the causal estimate. This procedure is briefly described in the “Supplementary Methods / Mendelian randomization” section. We believe that this procedure accounts for the possible inflation of the weights used to build the GRS. Furthermore, our MR analyses of the effect of heart rate is consistent with other methods that have been estimated using individual variant summary statistics (Supplementary Table 6).

Comment #3: On page 10 the authors describe how adjusting for atrial fibrillation would be problematic if both the SNP and heart failure collide on AF, which would introduce collider bias. However, isn’t it equally likely and problematic that, if the SNP has an effect on AF, and AF has an effect on heart failure, that both the SNP and the confounders of the AF-heart failure association would collide on AF (leading to collider bias when you adjust for AF)?

Author response: We believe that the reviewer’s assessment that adjusting for AF in estimating the SNP-HF relationship would be susceptible to collider bias.

In the manuscript, we propose the cause-specific competing risk model of incident HF as a means of estimating the effect of the SNP without the mediation by AF and without problematic statistical adjustment (Table 1). This result from the competing risk model is also supported by the mtCOJO analysis reported by the HERMES case-control consortium as mentioned in the “Results / Genetically predicted effect on efficacy endpoints” section (see ref. 14). We believe that the evidence from this method that is robust to collider bias as well as our results from the competing risk model reduce the risk of effect distortion due to collider bias.

Comment #4: The interpretation of main effect estimates are less straightforward when interaction effects have been added to the model. Therefore, please be explicit how the reader should interpret the sentence providing both the main and interaction (with AF) estimates between rs8038766 and heart failure.

Author response: We thank the reviewer for this comment which helped present the results of this analysis more clearly. We have adapted the text to report the raw beta coefficients as well as the OR for the SNP in individuals with and without atrial fibrillation. Specifically, the text now reads:

“In a model including the interaction term between rs8038766 and atrial fibrillation, the estimated coefficient for the variant was -0.136, 95% CI -0.205, -0.067, (P=0.00011) and the interaction term coefficient was 0.110 95% CI 0.005, 0.215 (P=0.04). These coefficients correspond to an estimated OR of the SNP on heart failure of 0.87 in individuals without atrial fibrillation and 0.97 in individuals with atrial fibrillation.”

Comment #5: Figure 1: Why was rs7174098 used for just one outcome in MEGASTROKE?

Author response: The selected variant, rs8038766, was not in the MEGASTROKE results file for large artery stroke which we downloaded from http://www.megastroke.org/. It was available for the other outcomes. We selected a very highly correlated variant to be used instead. We added this precision to the manuscript. The Figure legend now reads:

“* rs7174098 (LD r2=1 in 1000 genomes Europeans) was used instead of rs8038766 as the latter was unavailable in the MEGASTROKE summary statistics for this outcome.”

Comment #6: Supposedly you choose the transethnic data of MEGASTROKE for its large number of cases. However, rs8038766 need not necessarily be a strong genetic proxy for HCN4 in non-Europeans. Does using METASTROKE's European dataset show comparable results?

Author response: The European dataset did show similar results and we indeed decided to only present the trans-ethnic results for brevity and given the larger sample size.

Specifically, the European OR for cardioembolic stroke was 1.07 (1.01, 1.12) p = 0.013.

Comment #7: Please report (perhaps in supplemental material) whether the various GWAS meta-analyses were based on incident or prevalent cases and whether recurrent events were included.

Author response: Large GWAS meta-analyses typically use a combination of prevalent, incident and recurrent events in a case/control association model to increase the number of cases. For example, in the HERMES consortium analyses, 23 of the 51 included studies (45%) were RCTs suggesting that some events were incident and other were recurrent (or “worsening” events) with respect to the inclusion criteria of the underlying studies. Among the remaining studies are 16 (31%) cohort studies which included both prevalent and incident events. To summarize the data we used for large consortia was unspecific with respect to the event type and represented case/control analyses based on liberal case definitions.

We have not updated the manuscript because we feel that referring to the original publication is the best way to gain insight into the individual components of the meta-analysis on a case-by-case basis.

Comment #8: There exist MR methods which can incorporate correlated genetic variants to boost power. Did you consider these methods for HCN4-variants?

Author response: Because HCN4 is a short gene with a large constraint score (pLI = 1 in gnomAD), there are few functional variants that make good genetic instruments. In other words, we believe that there are few variants of functional relevance (refer to Comment #1 from reviewer #2 above). Additionally, Figures 1 and 2 report observational effects for a single variant based on various studies. We believe that for these non-MR estimates using a single variant adds clarity and avoids difficulties that could arise from using weighted scores like variability in variant availability in summary statistics, differences in LD patterns between studies. The other MR analyses we conducted were not “cis-MR” analyses and included independent variants from across the genome.

Comment #9: Table 1 describes the ‘genetic model for SHIFT/SIGNIFY’. Please be explicit this only refers to the outcome definition (and intervention), i.e., not also the inclusion criteria for participants.

Author response: We thank the reviewer for this suggestion that clarifies an important aspect of our manuscript. We have added the following text to the Table 1 legend:

Comment #10: The selected population of the UK Biobank may give rise to issues like selection bias, also for MR studies. Could this have influenced your results?

Author response: We agree with the reviewer’s comment. In the discussion, the manuscript reads “The MR estimates from the UK Biobank are also based on mostly healthy individuals with a low heart rate (mean of 69 bpm) possibly limiting clinical interpretation 33.” Which aims at highlighting the enrichment of healthy volunteers in the UK Biobank. In the “Study limitations” section, we also mention the possible differences in ethnicity and clinical profile: “We also used data from individuals of predominantly European ancestry both in the UK Biobank and in summary statistics from large GWAS consortia which could limit the generalizability of our results to other populations both in terms of clinical profile and ancestry.”.

We believe that these aspects are to be considered in the interpretation of our findings. Nonetheless, we combined evidence from many independent data sources and analytical models and our results are in line with the observations from ivabradine RCTs. Hence, in this case, the risk of a selection bias that would lead to incorrect conclusions is low.

Minor:

Comment #11: Please mention that the heart rate GRS includes variants which are just relatively independent (r2<0.1). Lower r2 thresholds are now typically advised for MR (e.g., 0.001). In extension, please be more explicit regarding the independence of the various sets of instruments used in the bidirectional MR analyses.

Author response: We thank the reviewer for this comment and have adapted to text to highlight the r2 threshold. The text in “methods > statistical analyses” now reads:

“For the construction of the heart rate Genetic Risk Score (GRS), we used 64 previously reported genome-wide significant heart rate associated SNPs (with r2 < 0.1) 16.”

Additionally, we respectfully disagree with the advice of using LD threshold as low as 0.001. To support our claim, we simulated genotypes for two independent SNPs with MAF 0.2 and corresponding to 1000 individuals. We then repeated this simulation 10,000 times and estimated that the probability of obtaining an r2 less than 0.001 by chance in this setting is around 32%. The R code for this short simulation is as follows:

> set.seed(5)

> n <- 1000

> n.sim <- 10000

> mafs <- c(0.2, 0.2)

> r2 <- sapply(1:n.sim, function(x) cor(rbinom(n, 2, mafs[1]), rbinom(n, 2, mafs[2])) ** 2)

> sum(r2 > 0.1) / length(r2)

[1] 0

> sum(r2 > 0.01) / length(r2)

[1] 0.001

> sum(r2 > 0.001) / length(r2)

[1] 0.3159

In light of this result, we believe that a r2 threshold of 0.001 may be too strict and lead to the exclusion of independent variants.

For the bi-directional analyses, the variants were selected using an r2 threshold of 0.15 as described in the “Supplementary methods > bi-directional MR” section.

Comment #12: Perhaps of interest, FINNGEN could serve as an additional source of publicly available summary statistics for (ICD-code based) heart failure

Author response: We thank the reviewer for this suggestion which may be useful for future work. For the current project, we feel that the HERMES case/control and UK Biobank were sufficient to support our claims.

Comment #13: For the analyses with the GRS it seems two-sample methodology was applied in the one-sample MR setting (e.g., Supplemental Table 5). Perhaps of interest: this (very) recent preprint (https://www.biorxiv.org/content/10.1101/2020.05.07.082206v1) suggests that particularly the MR-Egger method may be easily biased in this setting. Calculating the I2 may give insight.

Author response: We thank the reviewer for pointing out this recent preprint. As the methods other than MR-Egger seem to be less affected by this problem, we feel that our results are unlikely to be strongly distorted as we report estimates based on various MR techniques (IVW, MR-Egger, contamination mixture model, MR-PRESSO) and the 2-stage estimate based on the GRS. As such, we are in favour of keeping the MR-Egger results presented alongside with the other models in Supplemental Table 5.

Comment #14: Not a fan of ‘non-significant trend’ – perhaps rephrase to ‘provided (very) weak evidence’?

Author response: We agree with the reviewer that this wording has been abused. In this context, we used it because we show, later in the manuscript, that there is indeed a significant association between the HCN4 variant and heart failure when the increased risk in atrial fibrillation is properly accounted for. To avoid the stigma associated with this term, we have reworded as follows:

“We tested for association of the heart rate-reducing allele at the HCN4 variant rs8038766 with combined prevalent and incident heart failure in the UK Biobank and found week evidence of association (OR = 0.96, 95% CI 0.91-1.00, p=0.071) (Figure 2).”

Comment #15: With regard to the kinship threshold of >0.0884 – wouldn’t this reflect 2nd degree relationships or more (rather than ‘or less’)?

Author response: We thank the reviewer for spotting this mistake. The text now reads: “To avoid including related individuals, we used the kinship estimates from the UK Biobank and randomly selected an individual for pairs with a kinship coefficient > 0.0884 consequently keeping only individuals with 2nd degree relationships or less.”.

Comment #16: Please note that quintiles are not the group themselves but rather the cut-offs to define these groups

Author response: We thank the reviewer for this comment.

The “Methods / Statistical analyses” section now reads: “We split the participants based on the GRS quintiles with the group formed by the 5th GRS quintile (and above) corresponding to the higher heart rate group and the odds ratio for CAD, heart failure and atrial fibrillation were obtained by comparing the first 4 groups individually to the 5th group used as reference in logistic regression.”.

We have also adapted the wording in the Supplementary Figure 2 legend which now reads: “Effect of heart rate genetic risk score groups based on quintiles on atrial fibrillation, heart failure and coronary artery disease in the UK biobank dataset. For every outcome, the highest heart rate group (determined by the 5th quintile) is used as the reference group and the reported odds ratios are adjusted for age, sex and the first 10 principal components.”.

Comment #17: The abbreviation GRS is introduced twice in the methods section

Author response: We thank the reviewer for noticing this. The mention of “GRS” in the “Methods / Mendelian randomization” section now uses the acronym.

Submitted filename: response_to_reviewers.docx

Click here for additional data file.

1 Jul 2020

A genetic model of ivabradine recapitulates results from randomized clinical trials

PONE-D-20-09954R1

Dear Dr. Dubé,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Ify Mordi

Academic Editor

PLOS ONE

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

Reviewer #2: Thanks for responding to the comments in a comprehensive way. I particularly liked the justification of the R^2<0.01 threshold rather than R^2<0.001.

Best wishes, Stephen Burgess

Reviewer #3: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Stephen Burgess

Reviewer #3: No

7 Jul 2020

PONE-D-20-09954R1

A genetic model of ivabradine recapitulates results from randomized clinical trials

Dear Dr. Dubé:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Ify Mordi

Academic Editor

PLOS ONE