Genetic Evidence on the Association of Interleukin (IL)-1-Mediated Chronic Inflammation with Airflow Obstruction: A Mendelian Randomization Study

Abstract

Preclinical studies suggest interleukin (IL)-1α/β is involved in the pathogenesis of chronic obstructive pulmonary disease (COPD). However, recent trials of anti-IL-1 therapies showed limited benefit for COPD. To clarify, we primarily examined total and direct effects of IL-1 and its receptors/coreceptors/receptor antagonists (IL-1/IL-1Rs) on airflow obstruction (AO) using Mendelian randomization (MR), and secondarily explored reverse causation using bidirectional MR. We selected independent cis protein quantitative trait loci (cis-pQTLs) as genetic instruments for IL-1/IL-1Rs from two proteomic genome-wide association studies (n = 11,594) of European ancestry (mean age ∼47 years). We applied those cis-pQTLs to the International COPD Genetics Consortium (n = 15,256 cases, 47,936 controls) of ∼81.9% European descent (∼57 years). No IL-1/IL-1Rs were significantly associated with AO after correction for multiple testing. However, a higher genetically predicted IL-1 receptor antagonist (IL-1Ra) was nominally associated with a 20% reduction in AO risk using univariable MR, with a larger direct effect (∼31%, i.e. not via IL-1α/β) using multivariable MR. Furthermore, higher total IL-18 binding protein (IL-18BP) was nominally associated with lower AO. Nominal total effects were also noted for higher IL-1α with lower AO and higher IL-1R1 with higher AO. Higher IL-1Ra and IL-18BP might have a role in preventing AO, but need to be contextualized.

Supplemental data for this article is available online at https://doi.org/10.1080/15412555.2021.1955848 .

Keywords:

• Airflow obstruction
• IL-1
• IL-1 receptors
• IL-1 receptor antagonists
• Mendelian randomization

Introduction

Chronic obstructive pulmonary disease (COPD) is a leading cause of morbidity and mortality worldwide [1]. In addition to airway remodeling, COPD involves chronic airway inflammation, which is characterized by persistent respiratory symptoms and airflow obstruction (AO) due to the abnormal release of inflammatory cytokines and dysregulated immune cell activity [2]. Damage from inhaled noxious particles, such as from smoking, air pollution, and fuel combustion, induces interleukin (IL)-1 release, comprising IL-1α/β, IL-18, IL-33, IL-36α/β/γ, IL-37 [3], along with other cytokines. IL-1 is thought to be implicated in COPD by interacting with their receptors/coreceptors and receptor antagonists [2, 4]. COPD is incurable and existing treatments only relieve symptoms. Given anti-IL-1/IL-1Rs therapies are widely used as anti-inflammatories condition, such as rheumatoid arthritis, a better understanding of their etiological roles in COPD could inform their potential for preventing COPD [5].

Clinical trials of anti-inflammatory therapeutics, such as MEDI8969 [IL-1R1 monoclonal antibody] [6], canakinumab [IL-1β monoclonal antibody] [7], mepolizumab [IL-5 monoclonal antibody] [8], benralizumab [IL-5 receptor monoclonal antibody] [9], CNTO6785 [IL-17 monoclonal antibody] [10], MK7123 [cytokine receptor CXCR2 antagonist] [11], and tumor necrosis factor α antagonist [12], have shown limited benefit in preventing the development or exacerbation of COPD. The relatively small number of COPD cases in these trials, competing risks, and the inclusion of late-stage COPD patients may partially explain the null results [2]. Patients with stable COPD may have eosinophilic or neutrophilic inflammation [13,14], and elevated IL-1β and IL-17 [15], making it unclear whether IL-1-mediated chronic inflammation increases the risk of developing COPD or alternatively whether COPD promotes the dysregulation of IL-1/IL-1Rs.

Observationally, IL-1/IL-1Rs are related to the risk of developing COPD [2,4,16,17]. For example, a recent study summarizing preclinical evidence has suggested that IL-1α and IL-1β, dysregulating IL-1 signaling, play important roles in the pathogenesis of COPD [2]. Moreover, IL-1α, regulating the airway epithelium and lung fibroblasts, appears to be associated with COPD risk [17]. However, these studies are open to confounding by factors such as health status or smoking. Mendelian randomization (MR) provides an alternative approach by using genetic variants to proxy an exposure, so as to obtain less confounded estimates. We used MR to obtain total effects and direct effects of 12 interleukins (IL-1α/β, IL-1Ra, IL-1R1, IL-1Racp, IL-18, IL-18BP, IL-18Rα, IL-36α/β/γ, and IL-37) on AO; and secondarily to assess whether liability to AO affects these interleukins. Here, we considered AO as the outcome, which helps identify people at risk of COPD exacerbations and mortality [18], because only pre-bronchodilator spirometry is used by the International COPD Genetics Consortium (ICGC) to define COPD [19]. Pre-bronchodilator measurements are often used in clinical practice and epidemiological studies [1,8,9], although post-bronchodilator spirometry measurements are recommended by the Global Initiative for Chronic Obstructive Lung Disease (GOLD) [20]. We also used tuberculosis as a negative control outcome to assess the credibility of the MR estimates.

Material and methods

Study design

Here, we used two-sample univariable and multivariable MR to assess three possible relations between IL-1/IL-1Rs and AO. First, IL-1/IL-1Rs might affect the development of COPD (Figure 1(A)). Second, patients with AO might have higher IL-1/IL-1Rs than other people (Figure 1(B)). Third, the effect of IL-1 on AO could be partially mediated by its receptors and receptor antagonists (i.e. IL-1Rs, Figure 1(C)). As such, IL-1 would have an indirect effect on AO, perhaps via IL-1Rs. MR rests on three key assumptions. First, relevance, i.e. the genetic instruments predict the exposure. Second, independence, i.e. the genetic instruments are not confounded. Third, exclusion restriction, i.e. the genetic instruments are only linked to the outcome via affecting the exposure.

Figure 1. Directed acyclic graphs depicting three possible scenarios that could explain cause effects between interleukin (IL)-1/IL-1 receptors/receptor antagonists (IL-1Rs) and airflow obstruction (AO). (A) IL-1/IL-1Rs have effects on AO. (B) AO has reverse effects on IL-1/IL-1Rs. (C) IL-1 has an indirect effect on AO, perhaps via its receptors (IL-1Rs). Under the multivariable MR framework, the direct effect of IL-1 on AO is estimated as β, whilst its indirect effect is estimated as αγ. The total effect can be estimated as β+αγ $\mathrm{\beta }+\mathrm{\alpha }\mathrm{\gamma }$.

Data source

Proteomics GWAS studies

Two proteomic genome-wide association studies (GWAS) studies recruited in up to 11,594 European participants (including 3301 from the INTERVAL study [21] and up to 8293 from YFS and FINRISK survey [22]) with an average age of ∼47.0 years. The proteomic GWAS was adjusted for age, sex, body mass index, time between blood draw and processing, and the first three or ten genetic principal components. Genetic predictors of cis protein quantitative trait loci (cis-pQTLs, which affect protein abundance with little or no attenuated effect on messenger RNA or ribosome levels concerning trans-pQTLs [23]) were obtained from summary statistics for IL-1 (i.e. IL-1α [log-transformed relative fluorescent unit, log (RFU) [24]; IL-1β [standard-deviation (SD, ∼2.2 pg/ml in YFS)], IL-1Ra [SD, ∼3.5 ng/ml in FINRISK2002], IL-18 [SD, ∼38.0 pg/ml in YFS], IL-36α [log(RFU)], IL-36β [log(RFU)], IL-36γ [log(RFU)], and IL-37 [log(RFU)]) and IL-1Rs (i.e. IL-1R1 [log(RFU)], IL-1Racp [log(RFU)], IL-18Rα [log(RFU)], IL-18 binding protein [IL-18BP, log(RFU)]).

ICGC

ICGC recruited 82,438 participants (24,754 COPD cases and 57,684 controls) into a two-stage genome-wide association meta-analysis [19]. Considering that the definitions of COPD cases and non-COPD controls in ICGC were based on pre-bronchodilator spirometry data: forced expiratory volume in 1 s (FEV₁) < 80% and FEV₁ to forced vital capacity (FVC) ratio of <0.7 for COPD cases and FEV₁ > 80% and FEV₁/FVC > 0.7 for non-COPD controls, misdiagnosis of COPD may happen according to recommendations in GOLD [20]. Hence, we considered cases with AO and controls based on the pre-bronchodilator spirometry data in ICGC [18]. Specifically, in stage 1 (discovery), a total of 22 studies with genome-wide associations and 4 COPD case-control or cohort studies with spirometry, including 15,256 AO cases and 47,936 normal controls, were included in a fixed-effects meta-analysis after adjustment for age, sex, pack-years of smoking, ever and current smoking status, and ancestry-based principal components. In the stage 1 analysis, the mean age of participants was ∼56.4 years, ∼50.7% were women, and ∼81.9% were of European descent. Also, 71.0% were ever smokers with an average 22.2 pack-years. In the stage 2 analysis which included the UK BiLEVE study as a validation of the stage 1 findings, the mean age of participants was ∼57.9 years, ∼48.1% were women, 55.3% ever smokers with an average 22.5 pack-years. More detail on baseline characteristics of participants included in ICGC is described elsewhere [19]. We selected genetic variants for AO using ICGC stage 1 analysis to avoid potential bias due to possible population stratification [25], given the INTERVAL study and ICGC might have different allele frequencies for those cis-pQTLs and thereby confound the instrument-outcome association. We conducted MR analyses using the publicly available summary statistics from dbGaP under study accession: phs000179.v6.p2.

Statistical analysis

Instruments selection

Figure 2 shows the analytical framework of two-sample MR used in this study. We selected pQTLs potentially predicting IL-1/IL-1Rs on the basis of a threshold of p <5 × 10–6 $P<5\times {10}^{-6}$ [26]. We checked these pQTLs for accurate identification and description in RegulomeDB (https://regulomedb.org/regulome-search/) to reduce random variability [27]. Cis-pQTLs were identified by excluding pQTLs with expression quantitative trait loci (eQTLs) [23] in either RegulomeDB or PhenoScanner (http://www.phenoscanner.medschl.cam.ac.uk/) to avoid unknown pleiotropy. We identified independent cis-pQTLs (r2 < 0.01) based on the 1000 genomes European reference panel obtained from LDlink (https://ldlink.nci.nih.gov/). To reduce the possibility of selection bias and unknown pleiotropy, we further excluded SNPs associated with other causes of death and possible confounders of IL-1/IL-1Rs on COPD (e.g. smoking and eosinophils) in PhenoScanner at the genome-wide significance $(5\times {10}^{-8}$). We calculated the F-statistic as the square of cis-pQLTs on exposure divided by the square of its standard error and excluded those with an F-statistic less than the rule of thumb of 10 to minimize weak instrument bias [28].

Figure 2. The analytical framework of two-sample Mendelian randomization (MR) used in this study. Cis-pQTLs: cis protein quantitative trait loci; COPD: chronic obstructive pulmonary disease; ICGC: International COPD Genetics Consortium.

Two-sample MR analyses

We applied the identified cis-pQLT genetic variants to summary statistics from ICGC. We replaced poorly genotyped cis-pQLTs with available proxy variants in high linkage disequilibrium (r2 < 0.80) obtained using LDlink. We aligned genetic variants for IL-1/IL-1Rs and COPD on the same effect allele by flipping the non-palindromic strand variants, and removed palindromic (MAF threshold = 0.3) and incompatible variants. MR estimates for IL-1/IL-1Rs on AO were obtained by pooling Wald estimates (i.e. the ratio of the genetic variant-outcome effect and the genetic variant-exposure effect) using inverse-variance weighting with random-effects meta-analysis, to obtain total effects. We detected heterogeneity using Cochran’s Q-test.

In sensitivity analyses, fixed-effects with inverse-variance weighting, a weighted median [29], MR-Egger [30], and MR robust adjusted profile score (MR-RAPS) [31] were used. The weighted median method provides consistent estimates when at least 50% of the weight is from valid instrumental variants [29]. MR-Egger detects the presence of genetic pleiotropy using the nonzero MR Egger intercept [30]. MR-RAPS addresses systematic and idiosyncratic pleiotropy (i.e. balanced pleiotropic effects with an allowance of some large values), where no genetic variant satisfies the exclusion restriction assumption, and yields robust estimates in the presence of weak instruments [31].

We used bidirectional MR [32] to check the direction of causality (Figure 1(B)); and multivariable MR based on robust weighted linear regression [33] to investigate whether IL-1/IL-1Rs directly effects the risk of developing AO after adjustment for genetically predicted effects of other IL-1/IL-1Rs (Figure 1(C)). In the multivariable MR, only genetic variants independent across the relevant exposures were used. We calculated the likely lower bound of conditional F-statistic to assess multivariable instrument strength, assuming no correlation between IL-1 and IL-1Rs [34]. We also used multivariable MR Egger oriented on the exposure of interest and the modified Q-statistic to examine the presence of unknown pleiotropy with an assumption about no correlation between exposure, which yields conservative causal estimates [34]. Furthermore, we validated our results using summary statistics for doctor-diagnosed COPD from the UK Biobank [35].

Finally, we used tuberculosis from the FinnGen study [36] as a negative control outcome to cross-validate the MR estimates of IL-1/IL-1Rs on AO, which are particularly vulnerable to biases [37–40]. Tuberculosis is a negative control outcome because it not only shares the same potential biases but also has a known cause (mycobacterium tuberculosis infection) [41,42]. That is, IL-1/IL-1Rs would not be expected to affect the risk of developing tuberculosis in the absence of mycobacterium tuberculosis infection. As such, any associations found for IL-1/IL-1Rs with tuberculosis would indicate biased MR estimates.

Power analysis

Considering that genetic instruments explained a median variation in IL-1/IL-1Rs of 5.8% in the proteomics GWAS [21,22], the study has ∼80% power at 5% alpha to detect an odds ratio of 0.90 for AO and 0.61 for tuberculosis using the online calculator (https://shiny.cnsgenomics.com/mRnd/) [43,44].

Statistical software

All statistical analyses were conducted using R version 3.6.2 software platform (https://cran.r-project.org/) with the HalopR (https://github.com/cluoma/haloR), ieugwasr (https://mrcieu.github.io/ieugwasr/index.html), TwoSampleMR (https://mrcieu.github.io/TwoSampleMR/index.html), and robustMVMR (https://cran.r-project.org/web/packages/robustMVMR/index.html) packages. We used p values at the Bonferroni-corrected threshold to define statistical significance, which was computed as 0.05/3 (for three sub-IL-1 families)/2 (for two outcomes)/2 (for the bidirectional MR) = 0.004. p values between 0.004 and 0.05 were reported as nominal. We adhered to STROBE-MR: Guidelines for strengthening the reporting of MR studies for reporting our results [38].

Results

A total of 123 genetic variants were identified as instruments for IL-1/IL-1Rs from the two proteomic GWAS (Supporting Information Table S1), including one variant (rs61335305) that predicted both IL-1β and IL-1Ra. None of these genetic variants were associated with potential confounders or other causes of death at 5 × 10−8$5\times {10}^{-8}$. Of these 123 genetic variants, 57 remained after excluding palindromic variants and variants incompatible with COPDGene, as shown in Supporting Information Tables S2-S3. The average F-statistic ranged from 14.6 to 215.3, indicating lack of weak instruments. A total of 13 genetic variants were identified as instruments for liability to AO from ICGC stage 1. From which, we excluded rs17486278 because it is associated with cigarettes per day [45] a potential confounder of IL-1/IL-1Rs on COPD and rs754388 because it is associated with eosinophils count [46] and may be pleiotropic, leaving 11 genetic instruments for liability to AO (Supporting Information Tables S4-S5). The average F-statistic varied from 29.0 to 104.8, indicating instrument strength.

Total and direct effects of IL-1/IL-1Rs on AO

Figures 3–5 show estimates for 12 genetically predicted IL-1/IL-1Rs-mediated inflammation cytokines on AO and the negative control outcome of tuberculosis using both univariable and robust multivariable MR. No significant associations of IL-1/IL-1Rs with AO were noted after correction for multiplicity. However, genetically predicted higher IL-1α was nominally associated with a 13% reduction in risk of developing AO, with no heterogeneity identified (Q-statistic = 3.13, p = 0.792). The MR Egger intercepts showed no evidence of pleiotropy (Supporting Information Table S6). However, such an association diminished after adjustment for genetically predicted effects of other IL-1/IL-1Rs. Multivariable MR Egger intercepts indicated possible pleiotropy (p = 0.004), but not the modified Q-statistic (p = 0.675). The conditional F-statistic suggested weak instrumental strength, as shown in Supporting Information Table S7. Total effects based on the conventional multivariable MR [47,48] yielded similar estimates (Supporting Information Table S8).

Figure 3. Causal estimates of genetically predicted interleukin (IL)-1$\mathrm{\alpha }$α/β, IL-1 receptor antagonist (IL-1Ra), IL-1 recetpor-type 1 (IL-R1), and IL-1 receptor accessory protein (IL-1Racp) with airflow obstruction (AO) and the negative control outcome of tuberculosis using univariable and multivariable Mendelian randomization analyses. Multivariable MR* adjusted either IL-1$\mathrm{\alpha }$α/β or IL-1Ra; and multivariable MR** adjusted either IL-1$\mathrm{\alpha }$α/β, IL-1Ra, IL-1R1, IL-1Racp, IL-18, IL-18BP, IL-18Rα$\mathrm{\alpha }$, IL-36α/β/γ$\mathrm{\alpha }/\mathrm{\beta }/\mathrm{\gamma }$, or IL-37. The modified Q-statistic for multivariable MR* estimates was 6.63 with p = 0.675, and the corresponding value for multivariable MR** estimates was 33.48 with p = 0.945. IL-18: interleukin-18; IL-18BP: interleukin-18 binding protein; IL-18Rα: interleukin-18 receptor 1; IL-36α/β/γ: interleukin-36 α/β/γ; IL-37: interleukin-37.

Figure 4. Causal estimates of genetically predicted interleukin (IL)-18, IL-18 binding protein (IL-18BP), IL-18 receptor 1 $\mathrm{\alpha }$ (IL-18Rα$\mathrm{\alpha }$), and IL-37 with airflow obstruction (AO) and the negative control outcome of tuberculosis using univariable and multivariable Mendelian randomization analyses. Multivariable MR* adjusted either IL-18, IL-18BP, IL-18Rα$\mathrm{\alpha }$, or IL-37; and multivariable MR** adjusted either IL-1$\mathrm{\alpha }$α/β, IL-1Ra, IL-1R1, IL-1Racp, IL-18, IL-18BP, IL-18Rα$\mathrm{\alpha }$, IL-36α/β/γ$\mathrm{\alpha }/\mathrm{\beta }/\mathrm{\gamma }$, or IL-37. The modified Q-statistic for multivariable MR* estimates was 4.36 with p = 0.824, and the corresponding value for multivariable MR** estimates was 33.48 with p = 0.945. IL-1α: interleukin-1 α; IL-1β: interleukin-1 β; IL-1Ra: interleukin-1 receptor antagonist; IL-1R1: interleukin-1 receptor-type 1; IL-37: interleukin-37.

Figure 5. Causal estimates of genetically predicted interleukin (IL)-36α$\mathrm{\alpha }$, IL-36β$\mathrm{\beta }$, and IL-36γ$\mathrm{\gamma }$ with airflow obstruction (AO) and the negative control outcome of tuberculosis using univariable and multivariable Mendelian randomization analyses. Multivariable MR* adjusted either IL-36α$\mathrm{\alpha }$, IL-36β$\mathrm{\beta }$, or IL-36γ$\mathrm{\gamma }$; and multivariable MR** adjusted either IL-1$\mathrm{\alpha }$α/β, IL-1Ra, IL-1R1, IL-1Racp, IL-18, IL-18BP, IL-18Rα$\mathrm{\alpha }$, IL-36$\mathrm{\alpha }/\mathrm{\beta }/\mathrm{\gamma }$α/β/γ, or IL-37. The modified Q-statistic for multivariable MR* estimates was 1.52 with p = 0.997, and the corresponding value for multivariable MR** estimates was 33.48 with p = 0.945. IL-1α: interleukin-1 α; IL-1β: interleukin-1 β; IL-1Ra: interleukin-1 receptor antagonist; IL-18: interleukin-18; IL-18BP: interleukin-18 binding protein; IL-1R1: interleukin-1 receptor-type 1; IL-18Rα: interleukin-18 receptor 1; IL-37: interleukin-37.

Genetically predicted IL-1Ra was nominally associated with a 20% decreased risk of AO (Figure 3), with no heterogeneity (Q-statistic = 3.15, p = 0.207) or pleiotropy identified (p = 0.708). Such an association remained after adjustment for genetically predicted effects of IL-1α and IL-1β but diminished after adjustment for other IL-1/IL-1Rs. Genetically predicted IL-1β nominally increased risk of AO after adjustment for genetically predicted effects of other IL-1/IL-1Rs, but with no total effect. Genetically predicted higher IL-1R1 was nominally associated with a 25% higher risk of AO with no heterogeneity (Q-statistic = 0.84, p = 0.659) or pleiotropy (p = 0.536). However, such an association diminished after adjustment for genetically predicted other IL-1/IL-1Rs. Multivariable MR Egger intercepts (p = 0.913) and the modified Q-statistic (p = 0.945) did not indicate pleiotropy. Our results are consistent with the total effect of IL-1Ra on the risk of doctor-diagnosed COPD obtained from summary-level statistics in UK Biobank (Supporting Information Table S9).

Genetically predicted IL-18BP was nominally associated with a 27% decreased risk of AO (Figure 4), with no heterogeneity (Q-statistic = 0.07, p = 0.966) or pleiotropy (p = 0.867). Such an association remained after adjustment for genetically predicted effects of IL-18 and IL-37 but diminished after adjusting for other IL-1/IL-1Rs. Multivariable MR Egger intercepts (p = 0.041, Supporting Information Table S7) suggested possible pleiotropy, but not the modified Q-statistic (p = 0.945). Sensitivity analyses yielded consistent results (Supporting Information Table S6).

Total effects of liability to AO on IL-1/IL-1Rs

No associations of liability to AO on IL-1/IL-1Rs were observed (Figure 6). The MR Egger intercepts did not indicate pleiotropy. In addition, sensitivity analyses yielded consistent results, as shown in Supporting Information Table S10.

Figure 6. Causal estimates of genetically predicted airflow obstruction (AO) on interleukin (IL)-1 and their receptors and receptor antagonists using univariable Mendelian randomization analyses. IL-1α: interleukin-1 α; IL-1β: interleukin-1 β; IL-1Ra: interleukin-1 receptor antagonist; IL-18: interleukin-18; IL-18BP: interleukin-18 binding protein; IL-1R1: interleukin-1 receptor-type 1; IL-18Rα: interleukin-18 receptor 1; IL-36α/β/γ: interleukin-36 α/β/γ; IL-37: interleukin-37.

Negative control analyses

No associations of IL-1/IL-1Rs with tuberculosis were observed (Figures 3–5 and Supporting Information Table S6).

Discussion

Principal findings

This first MR study comprehensively assessing the roles of IL-1/IL-1Rs in the risk of developing AO provided genetic evidence for some suggestive causal associations of higher IL-1Ra and higher IL-18BP with decreased risk of AO. There was also some evidence for possible associations of higher IL-1α with decreased risk of AO and higher IL-1R1 with increased risk of AO. Furthermore, higher IL-1β was weakly associated with increased risk of AO independent of other IL-1/IL-1Rs.

Comparison with other studies

Our results are consistent with a previous review showing neutralizing antibodies to IL-18 (i.e. mimicking IL-18BP to block the effect of IL-18) have efficacy in preclinical models of inflammation and COPD or AO [49]. Evidence of higher IL-18BP decreased risk of AO was particularly strong after adjusting for genetically predicted effects of IL-18 and IL-37 (Figure 3). Moreover, IL-18 is associated with some auto-immune diseases, such as inflammatory bowel disease [50] and immune-mediated diseases [51], further supporting our results. Nevertheless, more testing of the IL-18 related pathways might be helpful to understand its role in AO better.

Strengths and limitations

This MR study considered several possible scenarios so as to obtain causal effects of IL-1/IL-1Rs on AO. We also used tuberculosis as the negative control outcome to assess the credibility of our instruments. Given the highly correlated IL-1 family, we used robust multivariable MR, which is less subject to conditionally weak instrument bias [52]. Furthermore, our findings are less likely to open to selection bias because of the relatively young average age (∼57 years) of participants in ICGC (i.e. few deaths would have occurred before the recruitment). The opposite direction concerning direct effects of IL-1β and IL-1Ra with AO further substantiates our findings, given the antagonistic effects of IL-1β and IL-1Ra on AO [2].

However, several limitations are noted. First, MR requires stringent assumptions of relevance, independence, and exclusion-restriction. To satisfy the relevance assumptions, we restricted the instruments to well-genotyped and curated cis-pQTLs strongly predicting IL-1/IL-1Rs. Any “winner’s curse” due to selecting instruments at a significance level of $5\times {10}^{-6}$ 5 × 10−6 might bias toward null [53]. To satisfy independence, we checked for potential confounders. Use of non-overlapping GWAS make our findings less likely biased by population stratification. To satisfy exclusion-restriction, we checked for known pleiotropy, but unknown pleiotropy might exist because of the limited understanding of these selected cis-pQTLs for predicting IL-1/IL-1Rs. Selection bias due to inevitably selecting survivors of their genetic make-up, AO and any competing risk of AO may remain. The multivariable MR Egger intercepts also suggest no genetic pleiotropy (Supporting Information Table S7). Second, genetic instruments typically only explain a small proportion of the variance in the exposure, as here. This might result in underpowered results, although multivariable MR often has more power than univariable MR by including more instruments at the cost of relatively more stringent assumption on horizontal pleiotropy. Third, the genetic instruments were adjusted for body mass index meaning that any effects of interleukins on AO via body mass index may not be evident [54]. However, such an effect might be limited as participants being unwell or with illness or infection were ineligible for donation in the INTERVAL study [55]. Fourth, canalization buffering genetic factors may exist and bias the causal estimates with unknown magnitude and direction. Fifth, effects of IL-1/IL-1Rs may vary by COPD subtype, including small airway disease and emphysema. We could not conduct subtype-specific analyses due to the lack of sub-type specific COPD GWAS. Sixth, given these GWAS mainly included healthy people of European descent, our findings may not extend to other populations, although causes are always consistent but may not be relevant in different populations [56]. Seventh, MR estimates proxy lifetime exposure to natural genetic variation rather than effects of short-term interventions.

In addition, we considered cases with AO and controls based on pre-bronchodilator spirometry in ICGC with fixed values of 80% for FEV₁ and 0.7 for FEV₁/FVC, instead of the 5th percentile of post-bronchodilator spirometry [20,57,58]. The mean age of participants in stage 1 in ICGC was ∼56.4 years (ranging from 40.5 to 79.5 years in sub-studies) [19]. We cannot rule out the possibility that using pre-bronchodilator spirometry may misclassify some asthmatic participants as having AO, whereas overdiagnosis in the elderly and underdiagnosis of AO in younger participants could occur due to the use of age-biased fixed ratio FEV1/FVC instead of the 5th percentile as an unbiased lower limit of normal [57,58]. Although, the outcome was classified in ICGC independent of the exposure (interleukin), inevitably, some non-differential misclassification of AO exists, which may bias toward the null [59]. To further validate our results, we applied the selected genetic instruments for IL-1/IL-1Rs to doctor-diagnosed COPD (1658 cases and 110,925 controls) in the UK Biobank and to early-onset COPD (diagnosed age <65 years; 1231 cases and 67,360 controls) and late-onset COPD (diagnosed age 60+ years; 1030 cases and 67,360 controls) in FinnGen using univariable MR. Consistent results were found, as shown in Supporting Information Table S9.

Mechanistically, IL-1Ra, competing for binding to IL-1R1 with IL-1α/β, is thought to regulate the function of IL-1α/β (e.g. blocking the effects of IL-1α/β), which contributes to repair, remodeling, and fibrosis of the lung [60]. IL-18 is partially regulated by IL-18BP in Th1 and Th2 cell development, inducing pulmonary inflammation and lung injury, and involving the pathogenesis of AO or COPD [61]. The alternative biological pathway (e.g. T helper lymphocytes-mediated immune response [62,63] and pattern recognition receptors-mediated innate immune response [64]) from IL-1/IL-1Rs to AO or COPD might exist as indicated by the consistent total effects of increased IL-1α and increased IL-1Ra on the decreased risk of AO. For instance, IL-1/IL-1Rs has been shown to be strongly associated with increased neutrophils and eosinophils as well as COPD exacerbations in COPD patients [14,65]. As such, despite the opposite effects of IL-1α and IL-1Ra on AO as would be expected, the interplay between IL-1 and its receptors/receptor antagonists are likely complex.

Public health implications

Primary prevention is important to reduce the COPD-related disease burden, but management and treatment of COPD remains challenging, largely due to being a heterogeneous, late-onset disease with complex pathophysiology [66–68]. Clarifying the role of interleukins provides additional insight. Replication using non-Western populations would be valuable to cross-validate our findings when such large-scale GWAS become available. IL-1Ra is increased by the drug anakinra, which is used as a treatment for rheumatoid arthritis, and possibly acts on immune function via raising testosterone in men [69], especially among those with metabolic syndrome and low testosterone [70]. Anakinra was also considered as a possible treatment for cardiovascular disease, but genetically mimicking anakinra appears to increase risk [71,72], which might need to be taken into account.

Conclusion

This MR study provides evidence suggesting increasing IL-1Ra and IL-18BP might have roles in preventing AO in people of European descent, perhaps via IL-1α/β, IL-18, and IL-37. A better understanding of the etiological roles of IL-1/IL-1Rs in chronic diseases could be relevant to reducing the worldwide burden of COPD.

Ethics approval and consent to participate

Each study has been specifically approved by the Ethical Committees of the original studies and all the participants provided a written informed consent. This analysis of publicly available data does not require ethical approval.

Authors contribution

ZY and MKK analyzed and interpreted the data regarding COPD, and were major contributors in writing the manuscript. All authors read and approved the final manuscript.

Consent for publication

Not applicable.

AO	=	airflow obstruction
cis-pQTLs	=	cis protein quantitative trait loci
CI	=	confidence interval
COPD	=	chronic obstructive pulmonary disease
FVC	=	forced vital capacity
eQTLs	=	expression quantitative trait loci
GWAS	=	genome-wide association study
ICGC	=	International COPD Genetics Consortium
IL-1	=	interleukin-1
IL-1α	=	interleukin-1α
IL-1β	=	interleukin-1β
IL-1R1	=	interleukin-1 receptor-type 1
IL-1Ra	=	interleukin-1 receptor antagonist
IL-1Racp	=	interleukin-1 receptor accessory protein
IL-18	=	interleukin-18
IL-18BP	=	interleukin-18 binding protein
IL-18Rα	=	interleukin-18 receptor 1
IL-36α	=	interleukin-36α
IL-36β	=	interleukin-36β
IL-36γ	=	interleukin-36γ
IL-37	=	interleukin-37
MR	=	Mendelian randomization
MVMR	=	multivariable Mendelian randomization
OR	=	odds ratio
RFU	=	relative fluorescent unit
SD	=	standard deviation

Acknowledgment

This research used data generated by the COPDGene study, which was supported by NIH grants U01 HL089856 and U01 HL089897. The COPDGene project is also supported by the COPD Foundation through contributions made by an Industry Advisory Board comprised of Pfizer, AstraZeneca, Boehringer Ingelheim, Novartis, and Sunovion. This research has been conducted using the data (dbGaP study accession: phs000179.v6.p2) granted through the Database of Genotypes and Phenotypes (dbGaP) under application number 23474. The authors also thank the participants and investigators of the INTERVAL and FinnGen studies.

Declaration of interest

The authors have declared no conflicts of interest.

Funding

This research received no specific grant from any funding agency in public, commercial, or not-for-profit sectors.

Data availability

The summary statistics for genetic associations of IL-1α, IL-1R1, IL-1Racp, IL18BP, IL18Rα, IL-36α, IL-36β, IL-36γ, IL-37 are available in the INTERVAL study (http://www.phpc.cam.ac.uk/ceu/proteins/). The summary statistics for genetic associations of IL-1β, IL-1Ra, IL-18 are available in YFS and FINRISK survey (https://grasp.nhlbi.nih.gov/FullResults.aspx). The summary genetic statistics for genetic associations of chronic obstructive pulmonary disease are available in the Genetic Epidemiology of COPD (dbGaP study accession: phs000179.v6.p2). The annotated QTLs are available in RegulomeDB (https://regulomedb.org/regulome-search/) and PhenoScanner (http://www.phenoscanner.medschl.cam.ac.uk/). The 1000 genomes European reference panel is available via LDlink (https://ldlink.nci.nih.gov/).

Code availability

Software codes and data are available on request.