key: cord-0725465-r58e9ojq
authors: Ahluwalia, Tarunveer S; Prins, Bram P; Abdollahi, Mohammadreza; Armstrong, Nicola J; Aslibekyan, Stella; Bain, Lisa; Jefferis, Barbara; Baumert, Jens; Beekman, Marian; Ben-Shlomo, Yoav; Bis, Joshua C; Mitchell, Braxton D; de Geus, Eco; Delgado, Graciela E; Marek, Diana; Eriksson, Joel; Kajantie, Eero; Kanoni, Stavroula; Kemp, John P; Lu, Chen; Marioni, Riccardo E; McLachlan, Stela; Milaneschi, Yuri; Nolte, Ilja M; Petrelis, Alexandros M; Porcu, Eleonora; Sabater-Lleal, Maria; Naderi, Elnaz; Seppälä, Ilkka; Shah, Tina; Singhal, Gaurav; Standl, Marie; Teumer, Alexander; Thalamuthu, Anbupalam; Thiering, Elisabeth; Trompet, Stella; Ballantyne, Christie M; Benjamin, Emelia J; Casas, Juan P; Toben, Catherine; Dedoussis, George; Deelen, Joris; Durda, Peter; Engmann, Jorgen; Feitosa, Mary F; Grallert, Harald; Hammarstedt, Ann; Harris, Sarah E; Homuth, Georg; Hottenga, Jouke-Jan; Jalkanen, Sirpa; Jamshidi, Yalda; Jawahar, Magdalene C; Jess, Tine; Kivimaki, Mika; Kleber, Marcus E; Lahti, Jari; Liu, Yongmei; Marques-Vidal, Pedro; Mellström, Dan; Mooijaart, Simon P; Müller-Nurasyid, Martina; Penninx, Brenda; Revez, Joana A; Rossing, Peter; Räikkönen, Katri; Sattar, Naveed; Scharnagl, Hubert; Sennblad, Bengt; Silveira, Angela; Pourcain, Beate St; Timpson, Nicholas J; Trollor, Julian; van Dongen, Jenny; Van Heemst, Diana; Visvikis-Siest, Sophie; Vollenweider, Peter; Völker, Uwe; Waldenberger, Melanie; Willemsen, Gonneke; Zabaneh, Delilah; Morris, Richard W; Arnett, Donna K; Baune, Bernhard T; Boomsma, Dorret I; Chang, Yen-Pei C; Deary, Ian J; Deloukas, Panos; Eriksson, Johan G; Evans, David M; Ferreira, Manuel A; Gaunt, Tom; Gudnason, Vilmundur; Hamsten, Anders; Heinrich, Joachim; Hingorani, Aroon; Humphries, Steve E; Jukema, J Wouter; Koenig, Wolfgang; Kumari, Meena; Kutalik, Zoltan; Lawlor, Deborah A; Lehtimäki, Terho; März, Winfried; Mather, Karen A; Naitza, Silvia; Nauck, Matthias; Ohlsson, Claes; Price, Jackie F; Raitakari, Olli; Rice, Ken; Sachdev, Perminder S; Slagboom, Eline; Sørensen, Thorkild I A; Spector, Tim; Stacey, David; Stathopoulou, Maria G; Tanaka, Toshiko; Wannamethee, S Goya; Whincup, Peter; Rotter, Jerome I; Dehghan, Abbas; Boerwinkle, Eric; Psaty, Bruce M; Snieder, Harold; Alizadeh, Behrooz Z
title: Genome-wide association study of circulating interleukin 6 levels identifies novel loci
date: 2021-01-30
journal: Hum Mol Genet
DOI: 10.1093/hmg/ddab023
sha: 5b619abd12e43d7309075cf86d63752d6bf30045
doc_id: 725465
cord_uid: r58e9ojq

Interleukin 6 (IL-6) is a multifunctional cytokine with both pro- and anti-inflammatory properties with a heritability estimate of up to 61%. The circulating levels of IL-6 in blood have been associated with an increased risk of complex disease pathogenesis. We conducted a two-staged, discovery and replication meta genome-wide association study (GWAS) of circulating serum IL-6 levels comprising up to 67 428 (n(discovery) = 52 654 and n(replication) = 14 774) individuals of European ancestry. The inverse variance fixed effects based discovery meta-analysis, followed by replication led to the identification of two independent loci, IL1F10/IL1RN rs6734238 on chromosome (Chr) 2q14, (P(combined) = 1.8 × 10(−11)), HLA-DRB1/DRB5 rs660895 on Chr6p21 (P(combined) = 1.5 × 10(−10)) in the combined meta-analyses of all samples. We also replicated the IL6R rs4537545 locus on Chr1q21 (P(combined) = 1.2 × 10(−122)). Our study identifies novel loci for circulating IL-6 levels uncovering new immunological and inflammatory pathways that may influence IL-6 pathobiology.

Dorret I. Boomsma 13, 14 , Yen-Pei C. Chang 12 , Ian J. Deary 25, 50 , Introduction Interleukin 6 (IL-6) is a multifunctional cytokine, which is involved in a wide range of immunomodulatory processes, from cellular migration and adhesion to proliferation and maturation (1, 2) . Interleukins are involved in immune cell differentiation and activation (3) . IL-6 is synthesized by a variety of different immune cells such as monocytes (4), B cells (5) and T cells (6) and also non-immune cells such as epithelial and smooth muscle cells (7) , adipocytes (8) , endothelial cells (9) and osteoblasts (10) .

Several factors have been implicated in circulating IL-6 levels. We have previously demonstrated that IL-6 levels decrease with age in children and increase with age in adults (11) . Also, increased levels of IL-6 have been observed in various diseases, not surprisingly in autoimmune diseases such as rheumatoid arthritis (12) and systemic juvenile idiopathic arthritis (13) , but also cardio-metabolic diseases like type 2 diabetes (14) , heart failure, coronary heart disease (15) and atherosclerosis (16) , as well in cancers (17) , atopic dermatitis (18) and psychological disorders like depression (19) . Due to its implications in the pathogenesis of different disorders, Il-6 has been used as an appropriate choice for drug targeting and used as a monitoring biomarker of disease progression and response to treatments (20) . The most illustrious IL-6 inhibitor is tocilizumab (21) , a monoclonal antibody binding the IL-6 receptor, which is already in use for treating patients with allergic asthma (22) , and immune system disorders like rheumatoid arthritis (23) and systemic juvenile idiopathic arthritis (24) , with high efficacy with some initial benefits towards respiratory illnesses like COVID-19 (25) .

IL-6 baseline levels are heritable with estimates from twin studies ranging between 15 and 61% (26) (27) (28) (29) . However, efforts to identify genetic variants associated with levels of IL-6 constituted relatively small-scale GWAS (30) (31) (32) (33) or sequencing-based candidate gene association studies (34) . To date, variants in the IL-6 receptor gene (IL-6R) and the gene encoding histo-blood group ABO system transferase (ABO) have been identified as statistically significant for an association to IL6-levels. Also, the genetic risk score constructed of IL-6 variants identified in the study by Shah and colleagues explained up to 2% of the variation in IL-6 levels (33), leaving a substantial part of its heritability unexplained. These seemingly sparse results and limited findings could be due to limitations in the study power caused by low sample size or a great inter-individual variability of IL-6 levels. One may speculate a substantial increase in the study size by increasing the number of participants, which would very likely lead to the identification of additional variants explaining IL-6 levels (35) (36) (37) .

The current study is the (till date) largest meta GWAS study including 67 428 individuals of European ancestry to identify genetic variants explaining the levels of circulating IL-6 and to understand underlying genetic mechanisms implicated in the pathophysiology of this cytokine.

A total of 52 654 individuals of European descent from 26 cohorts were included in the discovery GWAS meta-analysis with up to 2 454 025 autosomal SNPs passing quality control. Four cohorts (ALSPAC, MONICA/KORA, NTR and SardiNIA) identified genomewide significant associations in the ABO region, whereas none of the other 22 cohorts did, either individually or combined. These cohorts conditioned their results on their relevant top-SNP in ABO, the results of which were included in the discovery metaanalyses. The overall genomic control inflation factor (λ GC after correction) at the discovery stage meta-analysis was 1.0.

We identified 94 variants that were genome-wide significantly (P discovery < 5.0 × 10 −8 ; Supplementary Material, Table  S1 ) associated with IL-6 levels, representing two independent genetic loci on chromosomes 1q21 and 6p21. Two common SNPs (rs4537545 and rs660895), one per locus, Chr. 1q21 (IL6R), and Chr. 6p21 (HLA-DRB1/HLA-DRB5), showed the most significant association with IL-6 levels (index SNPs) and the third SNP (rs6734238) mapped on Chr. 2q14 (IL1F10/IL1RN) locus showed suggestive (5.0 × 10 −8 < P discovery < 1.0 × 10 −5 ) association in addition to 5 other loci (LHFPL3, LZTS1, GPC5/GPC6, USP32/APPBP2, STAU1; Supplementary Material, Table S2 ). Manhattan and QQ plots have been depicted in Figures 1A and 1B .

The minor alleles of IL6R rs4537545 * T (β = 0.091; P discovery = 8.39 × 10 −85 ), IL1F10/IL1RN rs6734238 * G (β = 0.025; P discovery = 1.45 × 10 −7 ) and HLA-DRB1/5 rs660895 * G (β = 0.036; P discovery = 1.80 × 10 −9 ) associated with increased circulating IL-6 levels ( Table 1) . Two additional genome-wide significant SNPs in the IL1R locus, rs11265618 (β = 0.047; P discovery = 1.21 × 10 −15 ) and rs10796927 (β = 0.034; P discovery = 1.24 × 10 −11 ), in low LD (r 2 < 0.25) with the lead SNP rs4537545 were carried forward for replication, and later conditional analysis as they seemed potential candidates as independent signals.

Overall, 12 SNPs spanning over 9 loci at a P discovery < 1 × 10 −5 in the discovery GWAS meta-analyses were selected for the replication stage (Supplementary Material, Table S2 ). This included the three index SNPs, two additional SNPs from the 1q21 locus (GWS but in low LD, r 2 < 0.25 with index SNP) plus an additional set of seven statistically suggestive SNPs with a P-value of 5 × 10 −8 < P < 1 × 10 −5 in the discovery meta-analyses (either in low LD, r 2 < 0.25 with the index SNP or independent loci). Additionally, 3 SNPs as negative controls and 3 SNPs in LD (r 2 > 0.25) with the Chr.1 index SNP, to control for possible genotyping errors of index SNP across replication cohorts, were also added to the replication list, yielding 18 SNPs for replication stage (Supplementary Material, Table S3 ).

Three loci including Chr.1q21 IL6R, Chr.6p21 HLA-DRB1/5 and Chr.2q14 ILF10/IL1RN replicated at P replication < 0.05, reaching GWS; 1q21 rs4537545, P combined = 1.20 × 10 −122 ; 6p21 rs660895, P combined = 1.55 × 10 −10 ; and 2q14 rs6734238, P combined = 1.84 × 10 −11 in the combined meta-analyses (Table 1 and Supplementary  Material, Table S3 ). Locus zoom plots available in Figures 1C,  1D , and 1E. The two additional signals at Chr.1q21 IL6R locus were replicated at P replication = 1.7 × 10 −4 for rs11265618 and P = 0.03 for rs10796927, reaching P combined = 2.5 × 10 −9 and P combined = 4.1 × 10 −13 , respectively (Supplementary Material, Table S3 ). The conditional analysis confirmed that rs11265618 and rs10796927 SNPs were not independent from (Supplementary Material, Table S4 ) but were driven by the index rs4537545 SNP.

In both, discovery and replication association analyses, the effect directionality was generally consistent across individual studies for GWS variants, while there was some evidence of borderline heterogeneity in one of the two novel loci (I 2 (P value) < 0.05) during the discovery and combined meta-analysis (Table 1 , Fig. 2) . The imputation quality scores (r 2 ) for the GWS (index) SNPs for each cohort (discovery and replication) are available in Supplementary Material, Table S5 . The other seven SNPs that showed suggestive association in the discovery stage, and expectedly the negative control SNPs did not reach GWS in the combined meta-analyses (Supplementary Material, Table S3 ). The three GWAS index SNPs when combined explained approximately 1.06% of the variance in circulating levels of IL-6 using data from the NESDA cohort. The phenotypic variance explained by all the common variants was estimated to be 4.45% using the SumVg method (38) .

IL6R was the only IL6 known locus that we replicate at GWS. IL1RN and HLA-DRB1, our primary findings have been reported as suggestive loci (1 × 10 −6 < P < 1 × 10 −4 ) by Shah et al. while some known/suggestive IL6 loci (ABO, BUD13, TRIB3 and SEZ6L) did not replicate (P discovery > 0.05) in the current study.

We looked up SNPs in LD with the index SNPs from the immunologically associated loci including IL-6R, rs4537545, 1q21; IL1F10, IL1RN, rs6734238, 2q14, intergenic; and HLA-DRB1/DRB5, rs660895, 6p21, intergenic. The search for functional/missense variants in high LD (r 2 > 0.8) with the lead SNPs led to the identification of only one nonsynonymous rs2228145 SNP in LD (r 2 = 0.95) with the rs4537545 index SNP from the IL6R locus. We used the Combined Annotation-Dependent Depletion (CADD) database to identify the functionality, i.e. deleterious, disease causal, pathogenicity, of rs2228145 in IL6R. CADD is an integrative annotation based on multiple genomic features scored into a single metric (39) . IL6R missense the rs2228145 variant has a CADD score of 15.98 (https://cadd.gs.washington.edu).

Genome-wide significant associations between IL6-associated top SNPs and other traits, and gene expression, data were mined using the Pheno Scanner v2 database (accessed, October 2020).

GWAS-based IL1F10/IL1RN rs6734238 * G allele has been associated with increased levels of serum C-reactive protein (CRP) and decreased fibrinogen levels, and blood cell traits in recent GWAS reports (40,41) (PMID:27863252; Supplementary Material, Table S6 ).

HLA-DRB1/DRB5 rs660895 * G allele is associated with increased risk of rheumatoid arthritis (RA) in Europeans and Asians (42) , IgA nephropathy in Asians (43) , while the decreased risk of ulcerative colitis and inflammatory bowel disease (IBD) (44) . IL-6R rs4537545 * T allele has been associated with increased circulating CRP levels (45), a decreased risk of RA (42) in mixed ancestries, while an increased risk of diabetes and asthma from the UK Biobank Neale's lab rapid GWAS (See Web Resources; Supplementary Material, Table S6 ). IL6R rs4537545T * allele is also associated with C-reactive protein, allergic disease, rheumatoid arthritis and coronary artery disease (Supplementary Material, Table S6 ).

Gene expression. IL1F10/IL1RN rs6734238 is associated with IL1F10/IL1RN expression levels in the skin, peripheral blood and whole blood (P < 5.0 × 10 −8 ; Supplementary Material, Table  S7 ). HLA-DRB1/DRB5 rs660895 has been associated with HLA-DRB1/DRB5/DRB6/DQB1/DQB2 expression levels in multiple tissues including peripheral blood, whole blood, monocytes, adipose tissue, thyroid, tibial artery, coronary artery, heart, lung, brain, colon, skeletal muscle, tibial nerve, skin and lymphoblastoid cell lines (P < 5.0 × 10 −8 ; Supplementary Material, Table S7 ). IL6R rs4537545 SNP is also associated with IL6R expression levels in peripheral and whole blood (P < 5.0 × 10 −8 ; Supplementary Material, Table S7 ).

Based on power calculator and assumptions mentioned under methods section, the estimated power for the 2 novel index SNPs was 98.3% rs6734238 (effect allele frequency, EAF: 0.42), and 76.9% rs660895 (EAF: 0.19), respectively.

We performed the largest (to date) GWAS meta-analysis for circulating IL-6 levels, which includes 66 341 individuals of European ancestry. We identified three loci associated with levels of circulating of IL-6 in the general population amongst which two are novel (Chr6p21, and Chr2q14), located in/nearby genes (HLA-DRB1 and IL1RN/IL-38) with inflammatory roles explaining up to 1.06% variance.

The strongest associated SNP, interleukin 6 receptor (IL-6R) rs4537545 at the 1q21 locus, is in high LD (r 2 = 0.95) with a missense IL-6R SNP rs2228145 (D358A) that results in an amino acid substitution at position 358 (Asp → Ala) on the extracellular domain of IL-6R and a high CADD score suggesting that the variant is pathogenic or functional or deleterious (among top 10% variants of the genome). The missense SNP is known to impair the responsiveness of cells targeted by IL-6 (46) by reducing IL-6R expression on cell surfaces (47) , and increasing levels of soluble IL-6R in individuals homozygous for this mutation (48, 49) . Recently it has been demonstrated that increased levels of sIL-6R induced by this variant can be explained by ectodomain shedding off IL-6R, a mechanism in which membrane-associated proteins are rapidly converted into soluble effectors whereby simultaneously cell surface expression of the same protein is reduced (50) . Increased levels of sIL-6R may act as a counterbalance to limit exaggerated IL-6 signaling and may explain the protective effect of the 358A allele for various cardiovascular diseases including coronary artery disease (CAD) (51-53), atrial fibrillation (54), lung function in asthmatics (55) and abdominal aortic aneurysm (56) as well as RA (57) . However, in contrast with this finding, the IL-6-sIL-6R complex itself is capable of transducing IL-6 signaling to non-IL-6R expressing cells, known as trans-signaling (58) , and it is this mechanism, as opposed to classic signaling, that is linked to chronic inflammatory disorders including IBD and RA (59) . Blocking IL-6 signaling cascades can be achieved by using an IL-6R specific inhibitor in the form of a monoclonal antibody, tocilizumab, which is a widely used therapy in the treatment of RA. Several variants in IL-6R, including rs2228145, may assist in the prediction of patient response to tocilizumab in RA (60) . The rs4537545 * T allele that is associated with IL6 levels is known to associate with increased circulating CRP levels (61) and a decreased risk of RA (42) in studies comprising mixed ancestries. Moreover, this SNP has been associated with IL6R expression in peripheral blood, skin, brain and adipose tissue (Supplementary Material, Table S7 ). The causal involvement of IL-6 levels in disease remains to be elucidated, but a recent study using a Mendelian randomisation (MR) approach did demonstrate that by using this SNP as instrumental variable, modelling the effects of tocilizumab, that IL-6R signalling has a causal effect on CAD (52) . On the other hand, pleiotropic nature of the IL-6R locus, influencing IL-6, CRP and fibrinogen levels, prohibits instrumental variable analysis and attribution of causality to one particular intermediate. Finally, several other genes encompass the 1q21 locus, including Src homology 2 domain containing E (SHE), and Tudor domain containing 10 (TDRD10), but have been ruled out to play a role at this locus (33) .

At the identified chromosome 2 locus the lead SNP, rs6734238, is intergenic and has also been associated with circulating CRP and fibrinogen levels (40, 41, 62) . The nearest genes to this locus are the interleukin 1 family member 10 (IL1F10, distance = 7.6 kb, currently known as IL-38) and interleukin 1 receptor antagonist (IL1RN, distance = 34.4 kb). IL1F10/IL-38 and IL1RN variants (rs6759676 and rs4251961) in partial LD with the lead SNP (r 2 LD :0.10 and 0.61) have been recently reported to be protective against the development of insulin resistance (63) . This further supports the molecular mechanisms behind IL-6mediated insulin secretion via glucagon-like peptide 1 (GLP-1) (64) contributing to type 2 diabetes (T2D) pathophysiology. For IL-6 specifically, it has been found that synthesis increases when dendritic cells are stimulated by bacterial lipopolysaccharides (LPS) in the presence of IL1F10 (65) . IL-1RN is another member of the interleukin 1 cytokine family, with suggestive evidence for involvement in determining IL-6 levels in the blood. One study found significant associations of IL-1RN rs4251961 with plasma CRP and IL-6 levels, albeit not independently replicated and not genome-wide significant (P = 1 × 10 −4 and P = 0.004) (66) . Our lead SNP was not in high LD (r 2 < 0.8) with variants in either neighboring genes and therefore in conjunction with its intergenic position, identifying a causal variant in this locus remains non-trivial.

The 6p21 rs660895, which was identified, resides within the HLA region, which forms one of the most complex genomic regions to study due to its large LD blocks and sequence diversity. This region has some population substructure in Europeans, which may have influenced the results; however, (1) each cohort population substructure adjustment was applied, followed by genomic correction for overall discovery stage metaanalyses. Thus, we reduced the chances that the population substructure may have had on this locus. The nearest genes to the index SNP, HLA-DRB1 (distance = 19.8 kb) and HLA-DQA1 (distance = 27.8 kb) are both histocompatibility complex genes encoding proteins that form cell surface complexes for certain immune system cells helping in antigen presentation to trigger an immune response. It is noteworthy that variations at this locus code for antigen-presenting complexes (APCs), which have been previously associated with diseases having a dysfunctional immune system; while we report for the first time that there exists also a strong association of this locus with circulating cytokine levels. Therefore, the association of this locus with the disease may corroborate through its effect through IL6 levels. One high-LD SNP (rs9272422, r 2 = 0.82 with our index SNP, rs660895) residing in the promoter region of HLA-DQA1 support this hypothesis and has been identified previously for systemic lupus erythematosus (67) and ulcerative colitis (68) . rs660895 * G allele is associated with increased risk of RA in Europeans and Asians (42) , IgA nephropathy in Asians (43) , while the decreased risk of ulcerative colitis and IBD. (44) Various studies aimed to identify genetic variation underlying levels of IL-6 (22-26) have found genome-wide significant associations in the IL-6R and ABO genes. The study performed by Shah and colleagues (33) found suggestive evidence (non-GWS; P = 3.8 × 10 −6 , respectively) for additional loci, including ABO, BUD13, TRIB3 and SEZ6L, none of which replicate in the current study (P discovery > 0.05) indicating that these might be false-positive findings due to low sample size (∼7800) or loci with sex-specific effects (associations based on women dominant population) or due to technical shortcomings with measurement assay (ABO locus).

It is surprising that even with increased statistical power (n discovery = 52 654; n replication = 14 774) in the current study (compared to the previous IL6 GWAS) (33) , we could identify three genetic loci (1q21, 2q14 and 6p21) accounting for ∼1% of the genetic variance for circulating IL-6 levels. According to the current estimates, the heritability levels for IL6 levels range between 15 and 61%, suggesting that an enormous increase in sample sizes would be required to identify additional variants explaining this remaining heritability. Multiple explanations for this socalled missing heritability phenomenon have been proposed in the past, which can be sought in rare or low frequency coding variants as observed for a similar metabolic quantitative trait by us (69) or can be explained by non-additive effects, which may cause inflated estimates of heritability. Plausible evidence for other sources of unexplained heritability that have been found are epigenetic changes, and haplotypes of common SNPs.

Collectively, our results provided additional insights into the biology of circulating IL-6. We identify new loci, limited by common variants in the Hap Map Reference panel. Albeit this is comparable to the 1000 genomes reference panel (70) but narrower compared to some newly available panels that show greater variant coverage in numbers and frequency range. Future studies are recommended to aim for identification of additional common but also rare variants, by firstly using richer imputation panels, such as UK10K project or the Haplotype Reference Consortium, a strategy that holds great promise, and secondly by making use of genetically isolated populations. Thirdly, we would like to stress the importance of phenotype harmonization. As we identified genome-wide variants in the ABO locus, in four studies participating in the discovery, but not in the remaining 22 cohorts, there is a strong indication that this locus may be assay specific. However, a proper demonstration of this hypothesis would require further testing, including repeating the GWAS in ABO-positive cohorts using a different IL-6 assay. Indeed it is emerging that the ABO locus has pleiotropic effects on many different traits and diseases (71) , which would suggest a more thorough analysis before disregarding this signal. Also, conventionally increasing sample sizes without correction for population substructures may raise heterogeneity within populations (72) , likely concealing the SNPs that affect particular subgroups. Future specific studies should counter the widely held assumption of unconditional risk alleles of complex traits and focus on the importance of studying more homogenous subgroups to, for example, investigate the age-dependent effect of genetic variants (73, 74) . Here, while further exploring the pleiotropic effect of IL-6-related variants, we identified phenotypes differentially regulated by diverse variants in the 1q21 locus. Biologic systems are dynamic complex networks and are evolving through lifespan and investigating the interrelationships existing between phenotypes as well as between genetic variations and phenotypic variations has the potential for uncovering the complex mechanisms. This is the case here for IL-6 and tailored methodologies should be devoted to the study of such traits, hopefully resulting in clinically significant breakthroughs. Future collaborative efforts therefore should strive to use well-calibrated assays, z-standardized protocols for sample handling, and processing (75) , though this will be difficult to achieve in practice. Lastly, we have attempted to perform formal association-based causal analysis to identify the likely causal loci, using the DEPICT approach; unfortunately, instead with only 2 novel GWS findings, our analyses were underpowered and thus not included. We also mine the gene expression and eQTL data for the identified SNPs using established databases; however, we were unable control for random co-localization signals or other confounders as we had limited access to summary level data. In conclusion, we identify two novel common genetic variants associated with circulating IL-6 levels that may influence the pathophysiology of complex cardio-metabolic, psychiatric and immunological traits, among individuals of European ancestry. This is a step further towards unravelling new biological pathways and potential therapeutic targets that can be developed for the IL-6-related disorders, while suggesting looking deeper into the genome for coding variants (rare and common) having larger individual effects (Figs 1 and 2 ).

Study populations. The overall study design ( Supplementary  Material, Fig. S1 ) involved the discovery cohorts with 53 893 individuals. After overlapping individuals with available genotype and phenotype data, the discovery stage included 52 654 individuals from 26 cohorts of European ancestry listed under Supplementary Material, Table S9 and in Supplementary Material, Text S1 and S2.

Each participating cohort performed genome-wide genotyping using a variety of genotyping platforms and applied a predefined quality control (QC) of genotype data (Supplementary Material, Table S11) followed by performing imputation of non-genotyped genetic variants, on the backbone of haplotypes inferred from the Hapmap Phase II reference panel (NCBI Build 36), and using statistical software such as IMPUTE (75) Table S11 ). Each cohort was recommended a set of general SNP quality filters including MAF < 0.01; Hardy Wienberg Equilibrium (HWE) P ≤ 10 −6 ; imputation quality r 2 ≤ 0.3; and genotyping call rate <0.95 (Supplementary Material, Fig. S1 ). Once we received summary results from each participating study, we ran a series of QC checks. Firstly, these included a set standard checks, including the imputation quality filters (basis the imputation program used and/or r 2 imputation < 0.3 were excluded), and then checks for genomic inflation (quantile-quantile or QQ plots). We adapted filter thresholds per cohort to reduce any observed deviation from the null while missing SNP loss due to the QC process. Finally, ∼2.45 million (2 454 025) common SNPs were part of the discovery meta-analysis. . The blue dashed line represents the null, and λ GC value represents the genomic inflation factor lambda. Each data point represents the observed versus the expected P-value of a variant included in the association analyses; (C-E) Regional association plots for each of the three genome-wide significant loci, 1q21, 2q14 and 6p21, respectively. Pairwise LD (r 2 ) with the lead SNP is indicated following a color-coded scale. Horizontal axis: relative genomic position of variants within the locus, vertical axis: −log10 P-value of each SNP. 

Each study conducted an independent GWAS analysis between SNPs and natural log-transformed values of serum IL-6 levels following a predefined analysis plan (Supplementary Material, Methods S4). Association analyses were conducted using linear regression model, or linear mixed effect models to account for familial correlation when warranted, with additive genetic effects, accounting for imputation uncertainty while adjusting for age, sex, population substructure (through study-specific principal components) and/or study-specific site, when necessary. GWAS summary result obtained from each cohort underwent a series of QC checks using the QCGWAS package in R (78) (Supplementary Material, Text S3 and Supplementary Material, Fig. S1 ). Being aware of the potential false-positive association in the ABO region on chromosome 9 (28, 30) , while using an R&D systems high-sensitivity assay kit to measure IL6 levels (R&D systems, Minneapolis, MN, USA), four (out of 22) discovery cohorts that observed genome-wide significant results in the ABO locus were asked to rerun the GWAS analysis conditional on the top ABO SNP (i.e. rs8176704) before including them in the final discovery meta-analysis (Supplementary Material, Text S3.3).

Individual GWAS results from 26 European studies were meta-analyzed using the inverse variance weighted, fixed-effects method as implemented in GWAMA while applying the double genomic control (GC) correction for population stratification, i.e. first to each study individually and subsequently also to the pooled results after meta-analysis (79).

Regional association plots for the discovered loci were generated through the LocusZoom (78) tool. We used the SNAP tool (80) to perform the pairwise LD checks (HapMap release 22 data) and to verify low LD with secondary signals. All SNPs selected for the replication stage had to fulfill the following criteria: (1) having an association P discovery ≤ 1 × 10 −5 and being in very low LD with the index SNP (r 2 < 0.2) and (2) available in at least 50% of study cohorts. Table S11 ) and statistical checks were made as in the discovery stage.

Venous blood samples (serum or plasma) were collected and stored at −80 • C. Serum/plasma IL-6 levels (pg/ml) were measured using various immunoassay methods described in Supplementary Material, Table S10 . Each cohort tested the selected SNPs using the same statistical model as for the discovery association analyses (Supplementary Material, Fig.  S1 ). Effect size estimates of all replication SNPs from each replication study were compared with the effect size estimates from the discovery meta-analyses. When effect sizes from individual cohorts did not align, we excluded these cohorts from the replication meta-analyses (n cohorts = 3). To account for the inter-study assay differences insensitivity, we combined results across the replication studies using a fixed effect sample size weighted Z-score meta-analysis as implemented in the METAL package (https://genome.sph.umich.edu/wiki/METAL) (81) .

The summary GWAS meta-analyses result from the discovery and replication stages were then used to perform a combined (discovery + replication) GWAS meta-analysis using a sample size weighted Z-score method. Test for heterogeneity was also performed as part of the meta-analysis package using METAL where I 2 statistic denoting the percentage of variation across studies was estimated (I 2 = 100% × (Q − df)/Q) where Q is the Chi-Square statistic. Significance for heterogeneity was denoted by the heterogeneity (or Het P ) P values. Variants that were significant in the replication meta-analysis at P < 0.05 and had an overall P combined < 5 × 10 −8 in the combined meta-analysis were considered statistically GWS. SNPs within the range of 1 Mb (or 10 6 bases) on either side of the most significant (i.e. index) SNP (with LD, r 2 > 0.25) were considered part of the same locus, whereas those in low LD (r 2 < 0.25) were tested if they were independent using conditional analysis.

We performed an approximate joint conditional analysis to identify distinct signals in a specific chromosomal region as implemented in GCTA (82) using high-quality genome-wide genotyped/imputed reference data from two studies (NEtherlands Study of Depression and Anxiety (NESDA) from the Netherlands and/or Genetics of Obesity in Young Adults (GOYA) from Denmark) to estimate linkage disequilibrium (LD) (83) between SNPs.

Conditional analysis for identification of independent signals was performed on GWS SNPs ( ±1 Mb to the index SNP and having low LD, r 2 < 0.2 with the index SNP) using summary statistics from the discovery GWAS meta-analysis data (-COJO option in GCTA) after confirming the GWS loci from the combined meta-analysis.

Heritability estimates. We approximated the variance explained by all distinct lead SNPs from the meta-analysis using the following formula:

where EAF is the effect allele frequency, β i the effect size of the individual variants, and n is the total number of lead variants. The current formula may overestimate the variance to a small extent as some level of SNP correlation was existent (LD r 2 < 0.25). The variance of the residuals of ln (IL-6) was calculated using data from the NESDA cohort (n = 2517). The total common SNP heritability of serum IL-6 levels explained by all GWAS variants was estimated using the observed Z-statistics from the discovery analyses for a subset of pruned SNPs. Following the original method (SumVg) (38) , we pruned the imputed (based on the 1000G Phase1 Integrated Release, Version 3, 2012 .04.30 reference panel) genotypes of the NESDA cohort using PLINK v1.07 (84) , by removing correlated SNPs at r 2 > 0.25 within a 100-SNP sliding window, and with a step size of 25 SNPs per forwarding slide. This resulted in a pruned set of 163 459 SNPs.

SNP mapping and functionality. We searched for variants in high LD (r 2 > 0.8) within a 1 Mb region on either side of the lead SNPs using 1000 Genomes sequence data (Phase1 Integrated Release, Version 3, 2012.04.30), utilizing tools available in Liftover (85) , VCFtools (86) and clumping in PLINK (84) . We subsequently annotated these variants using ANNOVAR (87) with the RefSeq (88) database for variant function and genic residence or distance. We used the Combined Annotation-Dependent Depletion (CADD) database to identify the functionality, i.e. deleterious, disease causal, pathogenicity, for the index SNPs.

Scanner v2 (89) data mining tool was used (Access date October 2020) to identify existing GWS (at P < 5 × 10 −8 ) associations between IL6 identified SNPs and other traits, and gene expression/eQTLs data.

We used GWAs power estimator (see Web Resources) by assuming a relative risk of 1.10 (or an effect estimate of 0.10), given N = 66 000, alpha (P-value) = 5 × 10 −8 (also GCTA power calculator, Supplementary Text S3.5).

Supplementary Material is available at HMG online. GWAS summary statistics data is available upon request to corresponding authors T.S.A. and/or B.Z.A. The UK Biobank Neale's lab rapid GWAS: (http://www.neale lab.is/blog/2017/7/19/rapid-gwas-of-thousands-of-phenotypesfor-337000-samples-in-the-uk-biobank).

The guarantors of the study are T.S.A., B.P. and B.Z.A., from conception and design to conduct of the study and acquisition of data, data analysis and interpretation of data. T.S.A. and B.P. have written the first draft of the manuscript. All co-authors have provided important intellectual input and contributed considerably to the analyses and interpretation of the data. All authors guarantee that the accuracy and integrity of any part of the work have been appropriately investigated and resolved and all have approved the final version of the manuscript. BP had full access to the data and the corresponding authors T.S.A. and B.Z.A. had the final responsibility for the decision to submit for publication.

IL-6: from its discovery to clinical applications

Interleukin 6: from bench to bedside

Evolutionary divergence and functions of the human interleukin (IL) gene family

Highly purified human peripheral blood monocytes produce IL-6 but not TNFalpha in response to angiotensin II

Autostimulatory effects of IL-6 on excessive B cell differentiation in patients with systemic lupus erythematosus: analysis of IL-6 production and IL-6R expression

Induction of IL-6 release from human T cells by PAR-1 and PAR-2 agonists

Human intestinal epithelial and smooth muscle cells are potent producers of IL-6

Release of interleukins and other inflammatory cytokines by human adipose tissue is enhanced in obesity and primarily due to the nonfat cells

Human endothelial cells produce IL-6. Lack of responses to exogenous IL-6

Mechanical loading highly increases IL-6 production and decreases OPG expression by osteoblasts

Biological variations, genetic polymorphisms and familial resemblance of TNFalpha and IL-6 concentrations: STANISLAS cohort

Serum interleukin 6 levels in rheumatoid arthritis: correlations with clinical and laboratory indices of disease activity

Correlation of serum interleukin-6 levels with joint involvement and thrombocytosis in systemic juvenile rheumatoid arthritis

Editorial: novel biomarkers for type 2 diabetes

Inflammatory markers and onset of cardiovascular events: results from the health ABC study

IL-6 in diabetes and cardiovascular complications

Versatile functions for IL-6 in metabolism and cancer

Protein-coding variants contribute to the risk of atopic dermatitis and skin-specific gene expression

Identification of IL6 as a susceptibility gene for major depressive disorder

IL-6 biology: implications for clinical targeting in rheumatic disease

Effects of interleukin-6 receptor blockade on allergen-induced airway responses in mild asthmatics

Efficacy of tocilizumab in patients with moderate to severe active rheumatoid arthritis and a previous inadequate response to disease-modifying antirheumatic drugs: the ROSE study

Efficacy and safety of tocilizumab in patients with systemic-onset juvenile idiopathic arthritis: a randomised, double-blind, placebo-controlled, withdrawal phase III trial

Association between early treatment with tocilizumab and mortality among critically ill patients with COVID-19

Genetic and environmental contributions to plasma Creactive protein and interleukin-6 levels-a study in twins

The age-dependency of genetic and environmental influences on serum cytokine levels: a twin study

Genetic architecture of the proinflammatory state in an extended twin-family design

Genetic and environmental determinants of population variation in interleukin-6, its soluble receptor and C-reactive protein: insights from identical and fraternal twins

A genome-wide association study identifies protein quantitative trait loci (pQTLs)

A genomewide association scan on the levels of markers of inflammation in Sardinians reveals associations that underpin its complex regulation

Novel genetic loci identified for the pathophysiology of childhood obesity in the Hispanic population

Gene-centric analysis identifies variants associated with interleukin-6 levels and shared pathways with other inflammation markers

Genome sequencing elucidates Sardinian genetic architecture and augments association analyses for lipid and blood inflammatory markers

PheWAS: demonstrating the feasibility of a phenome-wide scan to discover gene-disease associations

The use of phenomewide association studies (Phe WAS) for exploration of novel genotype-phenotype relationships and pleiotropy discovery

Phenome-wide association study (Phe WAS) for detection of pleiotropy within the population architecture using genomics and epidemiology (PAGE) network

Uncovering the total heritability explained by all true susceptibility variants in a genome-wide association study

CADD: predicting the deleteriousness of variants throughout the human genome

Multiethnic meta-analysis of genomewide association studies in >100 000 subjects identifies 23 fibrinogen-associated loci but no strong evidence of a causal association between circulating fibrinogen and cardiovascular disease

Meta-analysis of genome-wide association studies in >80 000 subjects identifies multiple loci for Creactive protein levels

Genetics of rheumatoid arthritis contributes to biology and drug discovery

A genomewide association study in Han Chinese identifies multiple susceptibility loci for IgA nephropathy

Association analyses identify 38 susceptibility loci for inflammatory bowel disease and highlight shared genetic risk across populations

Genetic loci associated with C-reactive protein levels and risk of coronary heart disease

Functional IL6R 358Ala allele impairs classical IL-6 receptor signaling and influences risk of diverse inflammatory diseases

Interleukin-6 receptor polymorphism is prevalent in HIV-negative Castleman disease and is associated with increased soluble interleukin-6 receptor levels

A common variant of the interleukin 6 receptor (IL-6r) gene increases IL-6r and IL-6 levels, without other inflammatory effects

Polymorphisms in the IL-6 receptor (IL-6R) gene: strong evidence that serum levels of soluble IL-6R are genetically influenced

Molecular and cellular mechanisms of ectodomain shedding

Interleukin-6 receptor pathways in coronary heart disease: a collaborative metaanalysis of 82 studies

The interleukin-6 receptor as a target for prevention of coronary heart disease: a mendelian randomisation analysis

Large-scale association analysis identifies new risk loci for coronary artery disease

Large-scale candidate gene analysis in whites and African Americans identifies IL6R polymorphism in relation to atrial fibrillation: the National Heart, Lung, and Blood Institute's candidate Gene Association resource (CARe) project

The IL6R variation asp (358) ala is a potential modifier of lung function in subjects with asthma

Interleukin-6 receptor pathways in abdominal aortic aneurysm

High-density genetic mapping identifies new susceptibility loci for rheumatoid arthritis

Interleukin-6 trans-signalling in chronic inflammation and cancer

Interleukin-6 trans-signaling and colonic cancer associated with inflammatory bowel disease

Interleukin-6-receptor polymorphisms rs12083537, rs2228145, and rs4329505 as predictors of response to tocilizumab in rheumatoid arthritis

Genome analyses of >200,000 individuals identify 58 loci for chronic inflammation and highlight pathways that link inflammation and complex disorders

Genomewide association and population genetic analysis of C-reactive protein in African American and Hispanic American women

Genetic determinants of circulating interleukin-1 receptor antagonist levels and their association with glycemic traits

Interleukin-6 enhances insulin secretion by increasing glucagon-like peptide-1 secretion from L cells and alpha cells

IL-38 binds to the IL-36 receptor and has biological effects on immune cells similar to IL-36 receptor antagonist

Polymorphisms of the IL1-receptor antagonist gene (IL1RN) are associated with multiple markers of systemic inflammation

Association of systemic lupus erythematosus with C8orf13-BLK and ITGAM-ITGAX

Host-microbe interactions have shaped the genetic architecture of inflammatory bowel disease

A novel rare CUBN variant and three additional genes identified in Europeans with and without diabetes: results from an exome-wide association study of albuminuria

Diagnostics and treatment of the diabetic foot

Detection and interpretation of shared genetic influences on 42 human traits

Next-generation genome-wide association studies: time to focus on phenotype?

Post-genome-wide association study challenges for lipid traits: describing age as a modifier of gene-lipid associations in the population architecture using genomics and epidemiology (PAGE) study

Birth cohort, age, and sex strongly modulate effects of lipid risk alleles identified in genome-wide association studies

A new multipoint method for genome-wide association studies by imputation of genotypes

Fast and accurate genotype imputation in genome-wide association studies through pre-phasing

Imputation-based analysis of association studies: candidate regions and quantitative traits

QCGWAS: a flexible R package for automated quality control of genome-wide association results

GWAMA: software for genome-wide association meta-analysis

Analyses and comparison of imputation-based association methods

METAL: fast and efficient meta-analysis of genomewide association scans

Common SNPs explain a large proportion of the heritability for human height

The Netherlands study of depression and anxiety (NESDA): rationale, objectives and methods

PLINK: a tool set for whole-genome association and population-based linkage analyses

The UCSC genome browser database: update

The variant call format and VCFtools

ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data

NCBI reference sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins

Pheno scanner: a database of human genotypephenotype associations

We acknowledge the participants of this study and the funding resources that have contributed to achieving this effort. Detailed acknowledgments can be found in Supplementary Material, Table S12 .

There is no conflict of interest from any of the participating co-authors. T.S.A. is a shareholder with Novo Nordisk A/S and Zealand Pharma A/S.

Novo Nordisk Foundation (NNF18OC0052457 to T.S.A.). Detailed Funding statements from participating study cohorts can be found in Supplementary Material, Table S12 .