Study design

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.

Included populations

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.

Definition of outcomes

Our primary outcomes were the incidence of sepsis, sepsis requiring critical care utilisation, and 28-day mortality after an episode of sepsis or sepsis requiring critical care admission, measured in UK Biobank.

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.

For the definition of critical illness related to infection, we utilised critical care admission data provided in HES. We considered any critical care episode during the index infection admission as a critical care admission related to that infection. Controls were defined as all other participants of UK Biobank. As HES data are only available for English participants of UK Biobank, we excluded all participants who were not recruited at a recruitment centre in England for the critical care portion of the study.

We generated 28-day survival outcomes for sepsis (with or without critical care admission) and pneumonia (only in those admitted to critical care). Dates of admission to hospital were extracted from the HES data and matched to national registry death data supplied by UK Biobank.

In order to attempt to avoid issues of bias relating to the structure of the data included in primary analysis (in particular, collider bias where selection on case status can induce associations between traits related to case status), we included all participants in our analysis of sepsis risk. For our analysis, we chose to use the whole UK Biobank population as a control for each outcome, even for the critical care related outcomes (e.g., we compared critically ill patients with sepsis versus the whole population, rather than critically ill patients with sepsis versus all patients with sepsis).

GWAS of infection outcomes

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.

Mendelian randomisation—Definition of instruments

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.

Statistical analysis

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].

For the UK Biobank and FinnGen outcomes, we performed fixed effects meta-analysis across each specific infection to generate summary estimates for each condition using the R package meta.

cisCRP variants

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].

After weighting, we performed two-sample MR analyses using this instrument on each outcome in turn, comparing results to that generated using the cisIL6R instrument.

cisgp130 instruments

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.

Sensitivity analyses

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].

Reporting guidelines

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

Ethics

Most of the data in this study were publicly available and nonidentifiable. Therefore, no ethical approval was required to access it. Access to UK Biobank was arranged by the UK Biobank IDAC. This study was performed under application number 56243.

Identification of sepsis cases in UK Biobank

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. IVW meta-analysis of MR estimates of IL6R blockade for UK Biobank and COVID-19 HGI outcomes.

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

Mendelian randomisation

We performed IVW MR on each outcome in turn. As our instruments are weighted by hsCRP, odds ratios (OR) are on the scale of natural log hsCRP decrease. For reference, we call this IL6R blockade, with ORs of more than one representing increasing risk with greater interference and ORs of less than one representing decreasing risk.

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

These severity-dependent effects were mirrored when we performed IVW MR on COVID-19-related outcomes from the COVID-HGI, with evidence suggesting that IL6R blockade is more protective from critical respiratory illness, defined as those who required respiratory support or who died during hospitalisation (OR 0.69, 95% CI 0.57 to 0.84) than from hospitalisation alone (OR 0.83, 95% CI 0.74 to 0.93). Notwithstanding uncertainty around the point estimates, MR estimates for IL6R blockade were larger in sepsis than in COVID-19, for the comparable disease severity.

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Follow-up analysis

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 specific sepsis definitions available in FinnGen—sepsis due to pneumonia and streptococcal sepsis—mortality data were not available, but these were rarer than the combined sepsis outcome (0.5% and 0.2%, respectively).

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.

Additional survival outcomes

Additional survival outcomes were available from two previous GWAS performed in the GAiNS and GenOSept consortium, both of which recruited patients with sepsis in critical care but also included a subset of patients with confirmed CAP. In both studies, case–control GWAS was performed with the outcome of 28-day survival.

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.

Alternative genetic instruments

Although our genetic instrument includes variants around IL6R (acting as a proxy for the target drug effect), we weighted it by the effect on hsCRP as this represents an appropriate readout of IL6R function and the effect on downstream IL-6 signalling. However, this does not mean that hsCRP is itself necessarily part of the causal pathway, and it may simply represent a measurable marker of IL6R blockade. It is also plausible that the beneficial effect of IL-6 blockade in sepsis is actually mediated by reductions in hsCRP. In an attempt to investigate this potential mediation and elucidate a potential mechanism of effect of IL6R blockade, we reran analyses using established cis variants around CRP known to be associated with circulating hsCRP levels.

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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

Outside of sepsis, there was little evidence of correlation between IVW MR results from cisCRP- and cisIL6R-related instruments (Pearson’s R 0.06), and, in addition, there was no strong evidence for an association with COVID-19 outcomes.

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).

The OR for sepsis with each SD increase in gp130 levels was 0.97 (95% CI 0.91 to 1.04), with similar estimates for the sepsis outcome in FinnGen (OR 1.04, 95% CI 0.96 to 1.14). Results for other infection outcomes were similar, with little clear evidence of an effect. The strongest evidence of any effect was in death from critical care pneumonia (OR 1.58; 95% CI 1.14 to 2.19) and death from critical care sepsis (OR 1.57; 95% CI 1.06 to 2.34), although effect estimates on death from hospitalised sepsis were close to the null (OR 0.95; 95% CI 0.81 to 1.13).

Noise-augmented clustering

As there was visual evidence of potential clustering of SNP effect estimates for our primary instrument (26 cis IL6R SNPs), we performed a form of mixture modelling (noise-augmented directional clustering) to identify clusters of instruments.

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.

Sensitivity analyses

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.