Mendelian randomization prioritizes abdominal adiposity as an independent causal factor for liver fat accumulation and cardiometabolic diseases.

Background: Observational studies have linked adiposity and especially abdominal adiposity to liver fat accumulation and non-alcoholic fatty liver disease. These traits are also associated with type 2 diabetes and coronary artery disease but the causal factor(s) underlying these associations remain unexplored.
Methods: We used a multivariable Mendelian randomization study design to determine whether body mass index and waist circumference were causally associated with non-alcoholic fatty liver disease using publicly available genome-wide association study summary statistics of the UK Biobank (n = 461,460) and of non-alcoholic fatty liver disease (8434 cases and 770,180 control). A multivariable Mendelian randomization study design was also used to determine the respective causal contributions of waist circumference and liver fat (n = 32,858) to type 2 diabetes and coronary artery disease.
Results: Using multivariable Mendelian randomization we show that waist circumference increase non-alcoholic fatty liver disease risk even when accounting for body mass index (odd ratio per 1-standard deviation increase = 2.35 95% CI = 1.31-4.22, p = 4.2e-03), but body mass index does not increase non-alcoholic fatty liver disease risk when accounting for waist circumference (0.86 95% CI = 0.54-1.38, p = 5.4e-01). In multivariable Mendelian randomization analyses accounting for liver fat, waist circumference remains strongly associated with both type 2 diabetes (3.27 95% CI = 2.89-3.69, p = 3.8e-80) and coronary artery disease (1.66 95% CI = 1.54-1.8, p = 3.4e-37). Conclusions: These results identify waist circumference as a strong, independent, and causal contributor to non-alcoholic fatty liver disease, type 2 diabetes and coronary artery disease, thereby highlighting the importance of assessing body fat distribution for the prediction and prevention of cardiometabolic diseases.

Introduction

Non-alcoholic fatty liver disease (NAFLD) is characterized by hepatic lipid accumulation ranging from steatosis (>5% of liver weight is lipids) to non-alcoholic steatohepatitis (NASH, presence of inflammation)[1]. Although liver steatosis may be relatively benign in most cases, more severe forms of NAFLD such as NASH and hepatic fibrosis can lead to liver cirrhosis and hepatocellular carcinoma. Approximately 25% of the adult population globally is affected by NAFLD with the prevalence rapidly increasing and potentially becoming the leading cause of liver failure in the United States by 2025[2,3]. Adiposity and body fat distribution are closely linked with NAFLD[4]. In observational studies such as the INSPIRE ME study, a large international imaging study using computed tomography, waist circumference was closely associated with liver fat accumulation independently of body mass index (BMI)[5].

Studies have also shown that both liver fat accumulation/NAFLD and waist circumference are associated with CAD and T2D[6-9]. However, whether liver fat accumulation is a causal factor of CAD and T2D remains to be elucidated and, more importantly, whether or not agents aimed at targeting NAFLD will ultimately decrease the risk of either T2D or CAD is unknown. In a previous investigation, we showed a strong genetic correlation between NAFLD, waist circumference, T2D, and CAD[10]. However, little is known about the directionality of these relations and whether NAFLD lies in the causal pathway linking abdominal adiposity and T2D/CAD.

In order to gain insight about the causality and directionality of these associations, causal inference methods such as Mendelian randomization (MR) have been developed[11]. MR uses genetic variants (which are randomly distributed at meiosis) such as single-nucleotide polymorphisms (SNPs), as instruments to infer causality. This method is in many ways comparable to a randomized control trial in which participants are naturally randomized based on the presence or absence of genetic variants that influence traits of interest[11]. In previous MR studies, a body fat distribution pattern consistent with low peripheral/subcutaneous fat accumulation and high intra-abdominal fat accumulation as estimated by the waist-to-hip ratio (WHR) adjusted for BMI was strongly associated with T2D and CAD[12,13]. However, we do not know if similar associations exist for NAFLD.

Extensions of the MR design, such as bidirectional MR and multivariable MR (MVMR), help in clarifying causal relations. Bidirectional MR refers to an analysis where both traits are alternately evaluated as exposure and outcome. This method has the potential to remove reverse causation bias by asserting the directionality of the relationship[14]. Multivariable MR can be used when multiple genetic variants are associated with two or more exposures. It conditions the effects of the SNPs of each exposure together to assess the effect of each exposure independently on the outcome. This method allows to test for mediation when two exposures share genetic variants as if they had been adjusted for one another[15].

Here, we used a MVMR study design to investigate the respective causal contributions of adiposity (defined using BMI) and abdominal adiposity (defined using waist circumference and the waist-to-hip ratio adjusted for BMI [WHRadjBMI]) to liver fat accumulation and NAFLD. Second, using a similar strategy, we aimed to determine if abdominal adiposity and liver fat accumulation are independent causal risk factors for T2D and CAD. Taken together, our triangulation of MR methods identify waist circumference as a strong, independent, and causal contributor to NAFLD, type 2 diabetes, and coronary artery disease.

Methods

Study populations

Information on the cohorts used in this MR framework is presented in Supplementary Data 1. Briefly, we combined data from publicly accessible GWAS summary statistics of European ancestry in a two-sample MR setting. BMI and waist circumference: The summary statistics of BMI and waist circumference were obtained from the UK Biobank from 461,460 and 462,166 individuals respectively. The GWAS was performed by the MRC IEU open GWAS project[16]. GWAS summary statistics from the GIANT consortium were also included to replicate the estimates obtained with the UK Biobank. These summary statistics for BMI were obtained from a meta-analysis of up to 125 GWAS for 339,224 European individuals[17]. Summary statistics for waist circumference were obtained from a meta-analysis of 232,101 individuals[18]. Measures of BMI and waist circumference were self-reported or measured in a laboratory or in a healthcare setting. Measures were corrected for age, age squared, sex, ancestry-based principal components, and study sites. The resulting residuals were inverse ranked normal transformed with standard deviation (SD) of 1. WHR adjusted for BMI: WHR adjusted for BMI was calculated as the ratio of waist and hip circumferences adjusted for BMI in 485,486 Europeans in the UK Biobank[19]. Measures of WHR and BMI were self-reported, measured in a laboratory or measured in a healthcare setting. Measures of WHRadjBMI were corrected for age, age squared, sex, principal components, and study site. The resulting residuals were transformed to approximate normality with SD of 1 using inverse normal scores. We also included GWAS summary statistics for WHRadjBMI from 210,088 Europeans from the GIANT consortium[18]. In that study, WHRadjBMI was adjusted for age, age-squared, study-specific covariates and then inverse ranked normal transformed prior to genome-wide analysis. NAFLD: We performed a GWAS meta-analysis for clinical diagnosis of NAFLD (8434 cases and 770,180 controls) of European ancestry from four cohorts, as previously described[10]. Briefly, we performed a fixed effect GWAS meta-analysis of The Electronic Medical Records and Genomics (eMERGE) network[20], the UK Biobank, the Estonian Biobank and FinnGen using the METAL package[21]. NAFLD was defined using electronic health record codes or hospital records. Logistic regression analysis was performed with adjustment for age, sex, genotyping site and the first three ancestries-based principal components. Liver Fat: GWAS summary statistics for liver fat were obtained from a GWAS of 32,858 white British participants from the UK Biobank[22]. Magnetic resonance scans were annotated by trained radiologists following a standard procedure. Using this training dataset, deep learning algorithms were then applied to estimate liver fat. The resulting dataset comprises 32,860 liver fat quantification. Liver fat was regressed using BOLT-LMM on gene carrier status, adjusted for genetic sex, age, age2, the first 10 principal components of genetic ancestry, scaled scan date, scaled scan time, and study center as fixed effects and genetic relatedness as a random effects term. The resulting residuals were inverse normal transformed prior to GWAS. Coronary artery disease: GWAS summary statistics for CAD were obtained from a GWAS on 122,733 cases and 424,528 controls from CARDIoGRAMplusC4D and UK Biobank[23]. Samples from CARDIoGRAMplusC4D were drawn from a mixed population (Europeans, East Asian, South Asian, Hispanic and African American), with the majority (77%) of the participants from European ancestry. Case status was defined by CAD diagnosis, including myocardial infarction, acute coronary syndrome, chronic stable angina, or coronary stenosis. We also used a different dataset GWAS summary statistics from the CARDIoGRAMplusC4D excluding UK Biobank (60,801 CAD cases and 123,504 controls)[24]. Type 2 diabetes: GWAS summary statistics for type 2 diabetes were obtained from the DIAbetes Genetics Replication and Meta-analysis (DIAGRAM) consortium and UK Biobank (74,124 cases/824,006 controls)[25]. Case status was defined by electronic health records, self-reports, or laboratory derived clinical diagnostics of T2D. We also used a different dataset from the DIAGRAM consortium excluding UK Biobank (26,676 T2D case and 132,532 controls)[26].

Some of the study samples used to derive our study exposures and outcomes included summary statistics from the UK Biobank, which lead to sample overlap. In univariable MR, sample overlap will bias the estimated results towards the null only when weak instrument is present. In MVMR, the direction of the bias is unclear but will occur only in the presence of weak instrument bias[27]. We included in our primary MR analysis the UK Biobank to increase power and included sensitivity analysis excluding the UK Biobank to remove sample overlap. All GWAS summary statistics were publicly available and accessible through URL. For all included genetic association studies, all participants provided informed consent, and study protocols were approved by their respective local ethical committee. Ethical approval was not required to conduct this study as it only used anonymized GWAS summary statistics.

Selection of genetic variants and variants harmonization

For univariable MR analysis, we selected all genome-wide significant SNPs (p-value < 5e−8). We then ensured the independence of genetic instruments by clumping all neighboring SNPs in a 10 Mb window with a linkage disequilibrium r2 < 0.001 using the European 1000-genome LD reference panel. SNPs and relevant association statistics can be found for each exposure in Supplementary Data 2. For multivariable MR analyses, we first extracted all genetic instruments that were previously selected for univariable MR analysis. We then pooled these SNPs to the lowest p-value corresponding to any of the exposures, using the same parameter setting as the univariable MR (r2 = 0.001 window = 10 Mb). We also included results of two other sensitivity analysis approaches: (1) prioritizing variants with lowest p value for BMI; (2) prioritizing SNPs with lowest p value for waist circumference. When NAFLD was used as an exposure in MVMR, we pooled the combined list of SNPs by selecting the SNP with the lowest p-value for NALFD. This procedure was implemented to select a maximum number of strong genetic instruments, as fewer genetic instruments are available for NAFLD exposure. SNPs in a 2 Mb window of the HLA, ABO, and APOE genetic regions were excluded due to their complex genetic architecture and their widespread pleiotropy (in GRCh37 6:28909037-30913661, 9:135130951-137150617, and 19:44409011-46412650, respectively). Exclusion of pleiotropic genetic regions satisfies the exclusion restriction and the exchangeability assumptions of instrumental variable analyses and strengthen inference of MR analyses. Harmonization was performed by aligning the effect sizes of different studies on the same effect allele. All GWAS summary statistics were reported on the forward strand. When a particular SNP was not present in the outcome datasets, we used a proxy SNPs (r2 > 0.6) obtained using linkage disequilibrium matrix of European samples from the 1000 Genomes Project. Instrument strength was quantified using the F-statistic[28], and the variance explained was quantified using the r2[29]. We calculated r2 for each individual SNP. For binary exposures, we calculated r2 using equation 10 in Lee et al., 2012[30] used in the get_r_from_lor function in the TwoSampleMR package. We calculated the F statistics following the formula . Where n is the sample size, k is the number of instruments used and R2 is the sum of the individual r2 of each SNP. These statistics can be found in Supplementary Data 3.

Statistical analyses

For univariable primary MR analysis, we performed the inverse variance weighted (IVW) method with multiplicative random effects with a standard error correction for under dispersion[31]. MR must respect three core assumptions (relevance, independence, and exclusion restriction) for correct causal inference. MR estimates bias occurs if the genetic instruments influence several traits on different causal pathways. This phenomenon, referred to as horizontal pleiotropy, can be balanced by using multiple genetic variants combined with robust univariable MR methods[32]. To verify if pleiotropy likely influenced the primary univariable MR results, we performed 6 different robust MR analyses: MR Egger[33], the MR-Robust Adjusted Profile Score (MR-RAPS)[34], the contamination mixture[35], the weighted median, the weighted mode and the MR-PRESSO[36], each making a different assumption about the underlying nature of the pleiotropy. Consistent estimates across methods provide further confirmation about the nature of the causal links. All continuous exposure estimates were normalized and reported on a SD scale. For dichotomous traits (i.e., diseased status on NAFLD, T2D and CAD), odds ratios were reported. Univariable MR analyses were performed using the TwoSampleMR V.0.5.6 package[37].

For multivariable primary MR analysis, we conducted the IVW method[38]. The use of MVMR is analogous to the inclusion of measured covariates in a multivariate linear regression. MVMR uses a set of overlapping genetic instrument to estimate the direct effect of an exposure on an outcome. As robust MVMR analyses, we used the multivariable MR-Egger[39], the multivariable median method, and the multivariable MR-Lasso method[40]. Similar to robust univariable MR analyses, each method makes different assumptions about the underlying nature of the pleiotropy and consistent estimates give confidence in the robustness of the causal findings. Multivariable MR analyses were performed using the MendelianRandomization V.0.5.1 package[41]. Conditionnal F-statistics were calculated with formula developed by Sanderson et al., in the MVMR V.0.2.0 package[42]. Percentage of mediation was quantified using the formula () Where is the direct effect estimated with IVW-MVMR and is the total effect estimated with univariable IVW-MR[43].