Therapeutic potential of IL6R blockade for the treatment of sepsis and sepsis-related death: A mendelian randomisation study

In this study, we aimed to perform a two-sample MR study in order to proxy the effect of IL6R blockade on sepsis and other infections. Fig 1 shows the overall study design, comparing our analysis with a randomised controlled trial of IL-6 receptor antagonist therapy.

For our main outcomes, we used UK Biobank, a large UK adult volunteer cohort described in detail elsewhere [24]. UK Biobank has linked genetic and physical data with direct links to national healthcare datasets. For secondary outcomes, we utilised FinnGen (Round 6), a large prospective cohort study in Finland, linked to electronic health record data [25]. For measurement of COVID-19 outcomes, we included data from the COVID-Host Genetics Initiative (HGI) (Round 7), a large meta-analysis of multiple studies including participants with COVID-19 [26].

Finally, for additional data on sepsis survival, we utilised summary statistics from a previous genome-wide association study (GWAS) on survival from sepsis, which included data from the GaINS and GenOSept consortium [27]. Full details of each cohort, inclusion criteria, and genetic quality control are available in S1 Text.

For secondary outcomes, we included (a) a set of 9 other common infections that present to primary or secondary healthcare (S2 Table) and (b) COVID-19 infection, as a comparator.

Admissions with sepsis were identified in UK Biobank in ICD-coded linked secondary care data. ICD-10 codes A02, A39, A40, and A41 were used to identify sepsis, in line with recent literature [28]. Cases were included if the code was in the primary or secondary diagnostic position in Hospital Episode Statistics (HES) data (or similar datasets in the devolved nations), provided by UK Biobank. We did not include self-reported cases, or cases only occurring in primary care. To ensure no contamination with COVID-19-related codes, we excluded codes for sepsis that occurred after 1 February 2020.

Other infections were defined similarly and included by the presence of ICD-10 codes, derived by two authors (FH and DA), with a code list available in the S2 Table. Codes were derived from recent publications but altered to reflect 3-digit ICD-10 coding available in HES. Controls were defined by the absence of the ICD code.

To generate summary statistics for downstream MR, we performed a case–control GWAS on each infection outcome using regenie v2.2.4 on all UK Biobank participants of European ancestry (see critical care exclusions above) [29]. We used the in-house MRC-IEU GWAS pipeline (version 2) to quality control our data [30]. Full details of this are published elsewhere [30], with further details in S1 Text.

Our IL6R instruments were selected based on a recent MR study that performed a meta-analysis of a high-sensitivity CRP (hsCRP) GWAS of 522,681 European individuals from the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium and UK Biobank [31].

Conceptually, our genetic instrument is similar to the action of anti-IL6R monoclonal antibodies (e.g., tocilizumab) that lead to inhibition of IL6R signalling by blocking both IL6 classical and trans-signalling, although the effect of inhibition is much reduced with our genetic instrument compared to therapeutic IL-6 inhibition [16–19,32,33].

Specifically, we aimed to assess the effect of decreased activity of the IL6/IL6R pathway, modelled using independent (r² < 0.1) variants within 300 kb of IL6R in order to proxy IL6RA effect. We call this our cisIL6R instrument. We weighted this variant on the effect on hsCRP, as justified in S1 Text, and in line with previous analyses [17–19]. For simplicity, we henceforth refer to this exposure as “IL6R blockade.”

It is recognised that this is an oversimplification of the IL6 pathway, with evidence that effects on health outcomes are mediated by classical and trans-signalling in differing ways and that the effect of our instrument may not act in the same way as IL6RAs [8]. In our sensitivity analyses (below), we explore other ways of defining our exposure, including alternative weighting strategies. In total, 26 variants were included. The minimum F-statistic of an included variant was 31.1 [34]. All SNPs were available in the GWAS performed above. The included SNP list is available in S1 Table.

For our main analysis, we identified and extracted SNPs (or proxies, for secondary cohorts) from our outcome GWAS and performed two-sample MR using harmonised SNPs on each outcome in term. MR estimates from each SNP were generated and then meta-analysed by inverse variance weighting (IVW) using first-order weights [21]. Analyses were performed using the TwoSampleMR package (version 0.5.6) in R (version 4.1.3) [35].

Our primary instrument includes only variants around IL6R, weighted by their effect on CRP, as a readout of IL-6 function. However, this does not imply that CRP itself is causal. In order to understand whether CRP might be the causal target, we undertook a subsequent analysis using four well-understood cisCRP-related variants from a recent study (S1 Table) [36]. These variants are highly likely to alter CRP through a pathway independent of IL-6 down-regulation, and, therefore, evidence of any MR effect would support CRP as a potential causal mechanism and also provide additional genetic support for targeting of this pathway. These variants were weighted using the same UK Biobank-CHARGE meta-analysis of hsCRP [37].

gp130, also known as IL6ST, is the other component of the IL-6 receptor, although it is also present in other cytokine receptors. Candidate gene analyses have suggested variation in gp130 has phenotypic consequences [38], although this is less established than the association with variation at IL6R. To explore this further, we generated instruments for gp130 plasma protein levels using a recent GWAS of plasma protein levels from the DECODE (n = 35,287) consortium. [39] We extracted independent (r² < 0.1) SNPs that were associated with levels of gp130 and that were within 300 kb of gp130 for use in downstream MR analyses on the infection outcomes using the TwoSampleMR package. Again, MR was performed in turn on each outcome.

We performed three broad types of sensitivity analysis. Firstly, we attempted to test one of the assumptions of MR (exclusion restriction) using alternative meta-analytic approaches (MR-Egger and weighted median approaches) [21], although these methods cannot falsify this restriction but can detect violations. Secondly, we tested whether variant weighting altered the results, by weighting variants solely on their effect estimates from CHARGE rather than from the UKB-CHARGE meta-analysis to avoid overfitting of data (“winners curse”) [21], and ran the analysis weighted entirely on the beta-coefficients from the SNP-outcome association (e.g., unweighted by CRP), with these results then meta-analysed together by the inverse variance weighting of each SNP-outcome (beta coefficient) association.

Subsequently, we tested whether our choice of SNPs to include altered the results. Firstly, we used the R package RadialMR to identify SNP outliers and reran analyses with outliers excluded [40]. Secondly, we performed iterative leave-one-out analyses, where each SNP is left out of the model in turn and IVW estimates recalculated. Thirdly, we reran the analysis including only variants within 10 kb of IL6R, and, finally, we ran the canonical and well-described Asp358Ala SNP (rs2228145) as a single instrument [16].

As an additional sensitivity analysis to explore heterogeneity identified in our primary analysis (using cisIL6R SNPS), we performed clustering to identify SNPs that had a consistent effect on multiple traits [41]. In this analysis, we identified all outcomes for which we had effect estimates for all 26 SNPs (31/35 outcomes, all except the GenOSept and GaINS consortium) and performed a form of mixture modelling (noise augmented directional clustering) to identify SNPs that had correlated effects on different traits and cluster them. This was performed using the navmix package in R using default settings (K = 10) on standardised associations between traits. Subsequently, we performed MR within each cluster to identify any evidence of heterogeneity [41].

This study is reported in line with the STROBE-MR guidance, with the checklist available in the Supporting information (S1 STROBE Checklist) [42].

In UK Biobank, we identified 11,643 cases of sepsis, with 474,841 controls of European ancestry. A total of 1,896 patients died within 28 days of admission, with 484,588 controls. A total of 1,380 patients had critical care admission with sepsis, with 429,925 controls (including only UK Biobank participants in England). Approximately 347 patients died within 28 days of critical care admission, leaving 431,018 controls. We subsequently performed a case–control GWAS for each of these outcomes, with links to quantile–quantile plots and Manhattan plots available in S2 Table, as well as details of other included infections in Table 1.

For our primary outcome—sepsis—we identified a severity-dependent effect, with evidence suggesting that IL6R blockade is increasingly protective with more severe disease. The OR for sepsis was 0.80 (95% confidence interval (CI) 0.66 to 0.96); for 28-day death in sepsis was 0.74 (95% CI 0.47 to 1.15), for sepsis requiring critical care admission was 0.48 (95% CI 0.30 to 0.78), and for 28-day death in sepsis requiring critical care admission was 0.37 (95% CI 0.14 to 0.98) (Fig 2A and Table 1).

Fig 2.

IVW MR estimates (ORs) with 95% CIs for IL6R blockade and each outcome (A: Sepsis, B: COVID, C: other infections) in UK Biobank. CI, confidence interval; COVID-19, Coronavirus Disease 2019; IVW, inverse variance weighted; LRTI, lower respiratory tract infection; MR, mendelian randomisation; OR, odds ratio; URTI, upper respiratory tract infection; UTI, urinary tract infection.

https://doi.org/10.1371/journal.pmed.1004174.g002

We performed IVW MR to investigate the association of IL6R blockade on the odds of infection in the absence of sepsis within UK Biobank. In line with trial and registry literature, there was evidence suggesting an increase in susceptibility to evaluated infections. The largest effect sizes observed were for upper respiratory tract infection (URTI) (OR 1.67; 95% CI 1.13 to 2.48) and for endocarditis (OR 1.66; 95% CI 0.88 to 3.15). For respiratory infections—lower respiratory tract infection (LRTI) and pneumonia—we saw no strong evidence of effect (OR approximately 1 for both) (Fig 2C and Table 1).

For both of these respiratory infections, despite the largely null effect with incidence of disease, IL6R blockade was associated with decreased odds of critical care admission (Fig 3A and Table 1), with effect estimates concordant with the estimates from sepsis requiring critical care admission: OR for pneumonia requiring critical care admission was 0.69 (95% CI 0.49 to 0.97) and the OR for LRTI requiring critical care admission was 0.51 (95% CI 0.23 to 1.21).

Fig 3.

IVW MR estimates of IL6R blockade with 95% CIs for (A) respiratory infection, (B) survival from sepsis related to critical care admission, and (C) FinnGen replication cohort. CI, confidence interval; IVW, inverse variance weighted; LRTI, lower respiratory tract infection; MR, mendelian randomisation; OR, odds ratio; URTI, upper respiratory tract infection; UTI, urinary tract infection.

https://doi.org/10.1371/journal.pmed.1004174.g003

We were able to extract parallel GWAS results from Round 6 of the FinnGen consortium (S3 Table) and to rerun all analyses. The frequency of incident cases differed across infection between studies. This was most notable for LRTI (frequency of incident cases 2.9% UK Biobank compared to 1.05% FinnGen). We identified three sepsis outcomes to compare: one combined one (all cases of ICD-coded sepsis) and two specific ones (subsets of the combined outcome). For the combined sepsis outcome, the frequency of incident cases was 54% higher in FinnGen (2.4% UK Biobank, 3.7% FinnGen), while the frequency of mortality at 5 years—the only mortality figure available within FinnGen—was 20.2% as opposed to 37.6% within UK Biobank.

For the combined FinnGen sepsis outcome, the OR was close to the null: OR 0.98 (95% CI 0.79 to 1.21). However, when focusing on the two specific outcomes of streptococcal sepsis and sepsis due to pneumonia, effect estimates were similar to those using UK Biobank: OR 0.79 (95% CI 0.48 to 1.31) and OR 0.63 (95% CI 0.29 to 1.35), respectively, although with more uncertainty due to a smaller sample size (Fig 3C).

Across UK Biobank and FinnGen, effect estimates were largely consistent in direction for all infections except pneumonia, although effect estimates were generally smaller. The meta-analysed summary estimate for sepsis, utilising UK Biobank and the main FinnGen sepsis outcome, had an OR of 0.86 (0.74 to 0.99), when meta-analysing with streptococcal sepsis in FinnGen, the OR was 0.80 (0.67 to 0.95), and when meta-analysing with respiratory sepsis in FinnGen, the OR was 0.79 (0.66 to 0.95). FinnGen-specific and meta-analysed summary effects are available in S4 Table, with a summary forest plot in S1 Fig.

In both studies, we performed IVW MR and effects were concordant with our primary data (Fig 3B), with a summary OR of 0.22 (95% CI 0.04 to 1.31) for sepsis and 0.06 (95% CI 0.01 to 0.55) for the CAP subset. These estimates were imprecise, given the number of included cases and controls.

In IVW MR analyses using the cisCRP instrument, we found evidence that CRP may be a part of the IL6R blockade effect, with evidence of reduced odds of sepsis and sepsis-related mortality: OR for sepsis outcome (0.91, 95% CI 0.82 to 1), OR for sepsis critical care admission (0.88, 95% CI 0.65 to 1.17), OR for death (0.72, 95% CI 0.59 to 0.93), OR for critical care death (0.58, 95% CI 0.32 to 1.04) (Fig 4 and S6 Table). This effect was also seen in FinnGen: OR for combined sepsis outcome (0.76, 95% CI 0.59 to 0.98), with similar results for other sepsis definitions. These effects were generally smaller in size than those seen using the cisIL6R instrument, with less evidence of a graduated effect with severity as seen with the cisIL6R instrument. Effect estimates were again broadly similar (although imprecise) for protection against critical infection (e.g., OR for critical care admission with LRTI 0.77, 95% CI 0.49 to 1.21).

Fig 4.

IVW MR effect estimates for CRP from cisCRP variants and each outcome (A: Sepsis, B: COVID, C: other infections). CI, confidence interval; COVID-19, Coronavirus Disease 2019; CRP, C-reactive protein; IVW, inverse variance weighted; LRTI, lower respiratory tract infection; MR, mendelian randomisation; OR, odds ratio; URTI, upper respiratory tract infection; UTI, urinary tract infection.

https://doi.org/10.1371/journal.pmed.1004174.g004

Subsequently, we investigated cis variants in gp130 (also known as IL6ST). Using a recently performed large-scale GWAS of gp130 plasma protein level (DECODE, n = 35,287), we identified 24 independent (r² < 0.1) SNPs (S1 Table). IVW MR was performed as described above in this analysis (reported in S7 Table and S2 Fig).

In this analysis, we identified one cluster of 14 SNPs (heatmap shown in S3 Fig) with all other SNPs excluded as noise (included SNPs listed in S8 Table). Subsequently, we performed MR using just these instruments (S9 Table and S4 Fig). These analyses showed similar results to our primary analysis, although effect estimates were generally larger, and p-values generally lower.

We performed a range of sensitivity analyses that did not alter the interpretation of primary results. Firstly, where possible, we tested some of the assumptions of MR using different meta-analytic techniques [21]. For these analyses, we performed meta-analysis using MR-Egger and weighted median approaches, which test the exclusion restriction assumption. Weighted median analyses were broadly similar to IVW results, but MR-Egger estimates were more imprecise, with nearly all confidence intervals for outcomes in UK Biobank, FinnGen, and COVID-HGI crossing the null. These results are shown in S5 Fig and S10 Table, while representative scatter plots for two representative outcomes (sepsis and critical care admission with sepsis) are available in S6 Fig.

In a second set of sensitivity analyses, we tested whether the way we designed and weighted our instrument materially affected the results. Firstly, we weighted our exposure SNPs using betas solely from the CHARGE consortium [34] (rather than using a combined CHARGE-UK Biobank meta-analysis; S11 Table). These results were similar (although more imprecise, due to the smaller sample size of the original GWAS). Secondly, we ran completely unweighted analyses (S12 Table and S7 Fig), which simply meta-analysed the effect of each SNP-outcome association, using IVW of the beta-coefficients at each SNP. These results were again concordant with our primary analysis, although were highly imprecise (given the small individual effect of each SNP).

In a third set of sensitivity analyses, we tested whether our selection of IL6R SNPs influenced results. Firstly, we identified and removed outliers using Radial MR and reran analyses (S8 Fig and S13 Table). We identified no alteration in our primary results. Alongside that, we performed leave-one-out analyses, removing individual SNPs from the analyses and rerunning IVW MR. For ease, we show this for the first four outcomes only in S9 Fig, which show the removal of any individual SNP does not alter the main results. Finally, we restricted analyses to the seven SNPs within 10 kb of IL6R (S14 Table) and the rs2228145 SNP, which has a known functional effect on IL6R and downstream effects [36] (S10 Fig and S15 Table). Again, these effects were all concordant with our primary results, although with reduced power as expected by removing multiple SNPs from the meta-analysis.