The Role of Inflammatory Cytokines as Intermediates in the Pathway from Increased Adiposity to Disease

Objective This study aimed to investigate the role of cytokines as intermediates in the pathway from increased adiposity to disease. Methods BMI and circulating levels of up to 41 cytokines were measured in individuals from three Finnish cohort studies (n = 8,293). Mendelian randomization (MR) was used to assess the impact of BMI on circulating cytokines and the impact of BMI‐driven cytokines on risk of obesity‐related diseases. Results Observationally, BMI was associated with 19 cytokines. For every SD increase in BMI, causal effect estimates were strongest for hepatocyte growth factor, monocyte chemotactic protein‐1 (MCP‐1), and tumor necrosis factor–related apoptosis‐inducing ligand (TRAIL) and were as ratios of geometric means 1.13 (95% CI: 1.08‐1.19), 1.08 (95% CI: 1.04‐1.14), and 1.13 (95% CI: 1.04‐1.21), respectively. TRAIL was associated with a small increase in the odds of coronary artery disease (odds ratio: 1.03; 95% CI: 1.00‐1.06). There was inconsistent evidence for a protective role of MCP‐1 against inflammatory bowel diseases. Conclusions Observational and MR estimates of the effect of BMI on cytokine levels were generally concordant. There was little evidence for an effect of raised levels of BMI‐driven cytokines on disease. These findings illustrate the challenges of MR when applied in the context of molecular mediation.


The Cardiovascular Risk in Young Finns Study
The Cardiovascular Risk in Young Finns Study (YFS) is a multicentre follow-up study with randomly chosen subjects from the Finnish cities of Helsinki, Kuopio, Oulu, Tampere, and Turku and their rural surroundings.
The study began in 1980 when 3,596 children and young adults participated in the first cross-sectional survey. The follow-up visits have been conducted in 1983,1986,1989,2001,2007, and 2011. The present cross-sectional study includes 1980 unrelated individuals who participated in the 2007 follow-up and who had both cytokine measurements and genotype data available. All participants gave written informed consent and the study was approved by local ethics committees 1 .

Laboratory analyses
The cytokine quantification was performed from semi-fasting (minimum 4 hours fasting) EDTA plasma in FINRISK 1997, from semi-fasting heparin plasma in FINRISK 2002 and from fasting serum in YFS. In YFS and FINRISK 2002, a total of 48 cytokines were measured by using Bio-Rad's premixed Bio-Plex Pro Human Cytokine 27-plex and 21-plex assays, as previously described [2][3][4] . In FINRISK 1997, a custom selected 20-cytokine multiplex assay was run. Cytokines were quantified for 2200, 2775 and 7906 individuals in YFS, FINRISK 2002 and FINRISK 1997. After merging cytokine, genotype and BMI data, final sample sizes used in the three cohorts were 1980, 1705 and 4608 individuals, respectively. Seven cytokines (out of the 48) were excluded from all subsequent analyses due to excessive missingness (>90%). This left a total of 41 measured cytokines in YFS and FINRISK 2002, 17 of which were also available in FINRISK 1997. For YFS, FINRISK 1997 and 2002, CRP was measured using high-sensitive turbidimetric immunoassay kit with an automated analyser.

Cytokine GWAS
As part of this work, a revised version of the previously published GWAS of 41 inflammatory cytokines 2 , without fitting BMI as a covariate in the model, was generated. The re-run was conducted on the same data as used previously (the Finnish cohorts described above) and all other aspects of the analysis were as described previously 2 . Leaving BMI out of the model increased the sample size by 44 individuals (39 in YFS and 5 in FINRISK 2002) who were excluded from the original GWAS and our Step 1 analyses because of having missing data on BMI.

Potential confounders
Description: Socioeconomic status was assessed using participant's educational background. In FINRISK 1997 and 2002 educational status was defined based on total years of full-time education and divided into three categories (low, average, high) and in YFS based on highest obtained degree (comprehensive school, secondary non-academic education and academic education). Smoking status was assessed using two categories, current smokers vs. ex and never smokers in FINRISK 1997 and 2002 and daily smokers vs. occasional and non-smokers in YFS. Alcohol consumption was used as a continuous covariate (average grams per week). To obtain approximate normality alcohol consumption was natural logtransformed prior to all analyses.
Part 1A -supplementary analyses: Associations of BMI, cytokines and CRP with potential confounders were estimated using linear regression. For BMI the confounder associations were adjusted for age and sex whereas for cytokines no additional adjustments were made due to the transformation process applied before analyses. Confounders tested were smoking status, alcohol consumption and socioeconomic status.
Associations were estimated separately for each potential confounder (used as the independent variable) and BMI, cytokines and CRP (as dependent variables). In addition, sensitivity analyses were conducted for the observational associations between cytokines and BMI in which potential confounders were fitted as additional covariables in the model.

Selection of genetic variants for cytokines
The genetic instruments for the BMI-associated cytokines taken forward to analysis step 2 were selected based on the results of our re-run of the previously published GWAS of 41 inflammatory cytokines (n=8337 from 3 studies) 2 . To maximise the number of available instruments for each cytokine, a list of cytokineassociated SNPs (p<5x10 -8 ) was first cross-matched to the corresponding GWAS summary statistics for each disease outcome in turn. During this harmonization process, SNPs for which strand alignment was ambiguous (A/T or G/C variants with intermediate allele frequencies) were removed. In cases where the associated SNP was not present, substitute SNPs (i.e. those in linkage disequilibrium with the associated SNP) were used where r 2 >0.8 (based on the European subset of 1000 Genomes). Then, taking the subset of SNPs that were present (or substituted) in both datasets, we performed linkage disequilibrium (LD) clumping using an r 2 <0.01 threshold in order to generate an independent set of instruments. This procedure was implemented using the TwoSampleMR R package 5 (version 0.4.26).

Derivation of body mass index (BMI) genetic risk score (GRS)
A weighted genetic risk score (GRS) for BMI was generated for use in the one-sample MR analysis (analysis step 1B) based on the 97 SNPs selected as instruments for BMI. A list of BMI SNPs and weights (i.e. betas extracted from the published meta-GWAS 6 ) that contributed to the BMI GRS instrument is given in Table S8. The sum of genotypes (coded 0,1,2) was multiplied by the strength of the effect of each SNP on BMI (in SD units) derived from the same BMI GWAS used to select the instruments 6 . The validity of the weighted GRS as a genetic instrument was assessed by examination of: i) first stage F-statistic which is a measure of the strength of the association between the instrument and exposure (calculated from the regression of BMI on the weighted GRS instrument in the first stage of the two stage least squares), and ii) associations of the GRS with potential confounding factors (smoking status, alcohol consumption and socioeconomic status).

Causal associations between BMI and levels of inflammatory variables: two-sample MR
For completeness, the causal associations between BMI and levels of inflammatory variables assessed in a one-sample framework in the primary analysis (as presented in the main text) were also evaluated using a two-sample framework. The genetic instruments used were the same as those selected for the one-sample MR described in the main text (Table S8). SNP-exposure associations were extracted from the same BMI GWAS as used for instrument selection 6 and SNP-outcome associations from our re-run of a GWAS of 41 inflammatory cytokines (n=8337) and results from a recent meta-GWAS of CRP (n>200 000); rs2245368 was not present in the CRP dataset leaving 96 SNPs as instruments 7 . For TNF-b four SNPs (rs16851483, rs17024393, rs13107325, rs11847697) were missing from the 97 BMI-related instruments. They were excluded from the cytokine meta-analysis because of low minor allele count in one of the two cohorts for which TNF-b was available and only SNPs for which results from more than one cohort was available, were included in the final summary statistics and used in the analyses.
For each instrument (SNP), causal estimates were derived using the Wald ratio method and estimates combined using the inverse-variance weighted (IVW) method to provide a single causal estimate of the exposure-outcome association 8 . To assess the potential for an effect of pleiotropic SNPs on the results, we performed the following sensitivity analyses: 1) MR-Egger regression method 9 to test overall directional pleiotropy and provide a valid causal estimate, taking into account the presence of pleiotropy; and 2) weighted median method 9 which provides a consistent estimate of causal effect if at least 50% of the information in the analysis comes from variants that are valid instrumental variables.

Power calculations
Retrospective power analyses were conducted to aid interpretation of null results. For the first step MR, power calculations were conducted using the online calculator available at (https://shiny.cnsgenomics.com/mRnd/) 10 . The parameters were as follows: alpha = 0.05, betaOLS = 0.08 (this was the median effect size seen in the observational analysis where p<0.05), R 2 xz = 0.027 (the variance in BMI explained by the GRS), variance of the exposure (X) = 1 and variance of the outcome (Y) = 1. The betayx parameter that gave power = 80% was then determined for the two relevant sample sizes (N=8,000 and N=3,500). For the second step MR, power calculations were calculated using the following formula 11 : Power = (pnorm(sqrt(N*R 2 xz*(ratio/(1+ratio))*(1/(1+ratio)))*betayx-qnorm(1-alpha/2)))*100 The parameters were as follows: sample size (N) = 50,000, ratio= 1:2 (ratio of cases to controls), alpha = 0.05 and R 2 xz (variance in trait explained by genetic variants) = 0.02, 0.05 and 0.025 for HGF, MCP1 and TRAIL, respectively. The betayx parameter that gave 80% power was then determined given N=50,000 (our threshold for inclusion) and assuming a 1:2 ratio of cases to controls.

Causal associations between BMI and levels of inflammatory variables: two-sample MR
Causal effect estimates of BMI on the 41 cytokine measures estimated from two-sample MR are shown in Figure S5. In line with the one-sample results presented in the main text, there was evidence of a causal effect of BMI on CRP and three cytokines: HGF, MCP-1 and TRAIL (p<0.001) with suggestive evidence for association with a further seven cytokines (p<0.05). There was no evidence of heterogeneity across the 97 BMI-associated variants (Q-statistic p>0.05) for HGF, MCP-1 and TRAIL ( Figure S6).
Using the 96 available SNPs as instruments, a one SD increase in BMI was associated with a 0.35 (95%CI: 0.24, 0.46; p<1.61x10 -10 ) increase in natural log-transformed CRP levels. Including the 93 available SNPs for HGF, a one SD increase in BMI was associated with a 0.32 normalized SD (95% CI: 0.17, 0.46; p<1.25x10 -5 ) increase in HGF. Using all 97 BMI-associated SNPs, a one SD increase in BMI was associated with a 0.26 (95% CI: 0.13, 0.40; p<0.0001) and 0.24 (95% CI: 0.11, 0.37; p<0.0004) normalized SD increase in MCP-1 and TRAIL, respectively. Results from the IVW analyses were consistent with the weighted median and MR Egger analyses for all four associated inflammatory variables (Figures S5 and Figure S1. Comparison of 27 SNP-cytokine association estimates with and without adjusting for BMI.

Supplemental Figures
BMI adjusted SNP-cytokine associations were obtained from previously published GWAS by Ahola-Olli et al. (2). Corresponding SNP-cytokine estimates were extracted from the re-run of the same GWAS without adjusting for BMI.      Figure S6. Mendelian randomization associations of body mass index (BMI) and four inflammation related variables (CRP, HGF, MCP-1, TRAIL) using a two-sample approach.
Effect estimates are given in normalized SD units per 1-SD higher BMI. Estimates derived from three alternative methods are presented. Q refers to the Cochran's Q heterogeneity statistic and p-value to the corresponding p-value.  Figure S7. Forest plots of SNP-specific Mendelian randomization associations between circulating TRAIL and odds of coronary artery disease.

SNP Combined
Heterogeneity: I 2 = 0%, 11   Interleukin-12p40 IL-12p40 21-plex assay a Custom-made multiplex assay data not included in our analyses for these cytokines b Cytokines with >90% of values missing were excluded Table S2. List of disease outcomes used in two-sample Mendelian randomization between BMI-driven cytokines and disease outcomes.   BMI was used as a dependent variable and association with each confounder (as independent variable) was analysed separately using linear model. All models are adjusted for age and sex. Low socioeconomic status was used as a reference level in the association between BMI and socioeconomic status. Non-smoking status (ex, never and occasional smokers) was used as the reference for the association between BMI and smoking; Alcohol consumption was measured in grams per week and natural log-transformed for the analyses (one was added to zero values). Betas reflect change in BMI in SD units. CI, confidence interval.

5.00E-07
Rank-normal transformed cytokines were used as a dependent variable and smoking status was treated as independent variable. Analyses were run using linear regression without further adjustments. Non-smoking status (ex, never and occasional smokers) was used as reference level. Betas reflect change in cytokine levels in normalized SD units. CI, confidence interval. Rank-normal transformed cytokines were used as a dependent variable and alcohol consumption was treated as independent variable. Analyses were run using linear regression without further adjustments. Alcohol consumption was measured in average grams per week and natural log-transformed for the models (one was added to zero values). Betas reflect change in cytokine levels in normalized SD units. CI, confidence interval Rank-normal transformed cytokines were used as a dependent variable and socioeconomic status was treated as independent variable. Analyses were run using linear regression without further adjustments. Low socioeconomic status was used as a reference level. Betas reflect change in cytokine levels in normalized SD units. CI, confidence interval  Measured BMI was used as a dependent variable and BMI related GRS as an independent variable when looking at association between BMI and BMI-GRS in linear model. Beta reflects change in BMI per unit change in GRS. BMI related genetic risk score (GRS) was used as a dependent variable and association with each confounder (smoking, alcohol consumption and socioeconomic status) as independent variable was analysed separately using linear model. All models are adjusted for age and sex. Low socioeconomic status was used as a reference level in the association between BMI related GRS and socioeconomic status. Non-smoking status (ex, never and occasional smokers) was used as a reference for the association between BMI related GRS and smoking; Alcohol consumption was measured in grams per week and natural log-transformed for the analyses (one was added to zero values). Betas reflect change in GRS. CI, confidence interval; GRS, Genetic risk score  Table S11. Results from one-sample Mendelian randomization analysis using natural log-transformed cytokine concentrations for a subset of cytokines that showed evidence of association with body mass index (BMI). Estimates correspond to change in natural log-transformed cytokine concentration (pg/ml) and CRP concentration (mg/l) per 1-SD change in BMI. Analyses were done using linear regression and are adjusted for age and sex. One was added to each cytokine value prior to transformation to account for zero values. Formula used to calculate F-statistic: (beta exposure) 2 / (SE beta exposure) 2 Formula used to calculate F-statistic (n,r2,k): ((n -k -1) / k ) x (r 2 / (1 -r 2 )), where n is sample size, k is number of instruments, and r 2 refers to variance explained by the instrument. Formula used to calculate variance explained: (beta exposure) 2 x (2 x frequency (1 -frequency)), where frequency refers to allele frequency. Cis variants for each inflammation related variable are bolded in the table.    Meta-analysis results of gender-specific SNP-inflammatory variable associations across the three cohorts obtained using GWAMA software (https://genomics.ut.ee/en/tools/gwama) with -sex option enabled. Gender differentiated pvalue refers to combined p-value of females and males assuming different effect sizes between genders and gender heterogeneity p-value to heterogeneity between genders. SNPs in bold are cis variants for corresponding inflammatory related variable. CHR, chromosome; SNP, single nucleotide polymorphism; se, standard error; Q, Cochran's heterogeneity statistic; I 2 , Higgins' heterogeneity index