Multiomics study of nonalcoholic fatty liver disease

Nonalcoholic fatty liver (NAFL) and its sequelae are growing health problems. We performed a genome-wide association study of NAFL, cirrhosis and hepatocellular carcinoma, and integrated the findings with expression and proteomic data. For NAFL, we utilized 9,491 clinical cases and proton density fat fraction extracted from 36,116 liver magnetic resonance images. We identified 18 sequence variants associated with NAFL and 4 with cirrhosis, and found rare, protective, predicted loss-of-function variants in MTARC1 and GPAM, underscoring them as potential drug targets. We leveraged messenger RNA expression, splicing and predicted coding effects to identify 16 putative causal genes, of which many are implicated in lipid metabolism. We analyzed levels of 4,907 plasma proteins in 35,559 Icelanders and 1,459 proteins in 47,151 UK Biobank participants, identifying multiple proteins involved in disease pathogenesis. We show that proteomics can discriminate between NAFL and cirrhosis. The present study provides insights into the development of noninvasive evaluation of NAFL and new therapeutic options.

We compared the PDFF and cirrhosis effects of the 18 NAFL variants in the UKB data. The cirrhosis effects were proportional to the PDFF effects ( Fig. 2) except for p.His48Arg in ADH1B (alcohol dehydrogenase 1b) and p.Cys282Tyr in HFE (homeostatic iron regulator). These two variants are associated with cirrhosis through alcohol consumption (ADH1B) 37 and hemochromatosis (HFE) 38 , rather than solely through hepatic fat. Similarly, the effects of the 18 variants on HCC were proportional to the PDFF effects (Fig. 2).

Potential causal candidate genes
To prioritize plausible causal genes at the associated loci, we evaluated the lead affected amino acid sequence, mRNA expression (expression quantitative trait loci (eQTLs)) or splicing QTLs (sQTLs). For this analysis, we used annotation of 46.5 million sequence variants tested in the GWAS and measured mRNA levels using inhouse RNA-sequencing (RNA-seq) data from whole blood (n = 17,846) and adipose tissue (n = 770), and publicly available data in the Genotype-Tissue Expression project (GTEx, v.v8). Fourteen lead associations were with missense variants or variants in high linkage disequilibrium (LD; r 2 > 0.8) with a missense variant: in PNPLA3, APOE, GCKR, GPAM, PNPLA2, TMC4, MTARC1, APOH, ADH1B, HFE, ERLIN1, GUSB and two independent variants in TM6SF2. The missense variants in GUSB, TM6SF2, PNPLA3 and TMC4 are associated with expression levels of the corresponding gene (cis-eQTLs) in various tissues, and the missense variant in TMC4 is also associated with liver expression of MBOAT7. Loss of MBOAT7 has been associated with NAFLD [39][40][41] . These variants are top eQTLs, that is, they are either the strongest association at their loci with expression levels of these genes or highly correlated to it (r 2 > 0. 8). The variants in TOR1B, HSD17B13 and GUSB associate with splicing in whole blood as top sQTLs (Supplementary Table 3). The intronic variant rs7029757[A], located 50-bp downstream of exon 2 in TOR1B (torsin 1b), is associated with cryptic splicing (P = 1.0 × 10 −1493 ; Supplementary Table 3). This variant elongates exon 2 by 50 bp, leading to a frameshift introducing a premature stop codon in exon 3 ( Supplementary Fig. 3).
It is challenging to both diagnose and stage NAFLD. Although liver enzymes are commonly elevated, they are nonspecific and poor predictors of progression [22][23][24] . Magnetic resonance imaging (MRI)-derived proton density fat fraction (PDFF) provides accurate liver fat quantification, but liver biopsy is essential for NASH diagnosis and staging. However, biopsy involves sampling variability and a risk of complications 3,25 . Currently, no pharmacological therapy has been approved for NASH. The identification of potential drug targets and biomarkers to monitor disease progression and treatment response is therefore paramount.
We performed genome-wide association studies (GWASs) of PDFF, NAFL, cirrhosis and HCC, and a combined GWAS of PDFF and NAFL. To further characterize the risk variants, we tested them for association with clinical traits, mRNA expression in various tissues and circulating protein levels in large datasets from Iceland and the UK Biobank (UKB). We further investigated whether plasma proteins discriminate between disease stages.

GWAS analysis
We extracted liver PDFF from raw abdominal MR images of 36,116 Britons of European ancestry from the UKB (Fig. 1, Table 1 and Methods). The MR images were generated using two acquisition techniques, that is, 8,448 images using gradient multiecho (GRE) and 27,668 using iterative decomposition of water and fat with echo asymmetry and least-squares estimation (IDEAL). The results of the two methods were similar, with 20.0% of IDEAL images and 23.6% of GRE images having a PDFF > 5%. Of the 270 (IDEAL) and 97 (GRE) individuals diagnosed with NAFL, 179 (66.3%) and 76 (78.4%) had a PDFF > 5.0% ( Supplementary  Fig. 1). Our PDFF estimates correlate well (r 2 = 0.96) with 3,869 PDFF measurements calculated by others using LiverMultiScan 26 (Supplementary Fig. 2).
We performed a GWAS on the PDFF estimates (n = 36,116) using each individual's first available MRI measurement. We also meta-analyzed GWASs on clinically diagnosed NAFL (International Classification of Disease, 10th revision (ICD-10) 27 , code K76.0), including 785 cases from Iceland (deCODE genetics), 5,921 from the UK (UKB), 2,134 from the USA (Intermountain INSPIRE and HerediGene registries) and 651 from Finland (FinnGen), for a total of 9,491 NAFL cases and 876,210 controls. We combined the summary-level GWAS PDFF data and the NAFL ICD-10 code meta-analysis data to maximize power to detect associations using multitrait analysis 28 .
We identified 18 independent sequence variants at 17 loci in the combined GWAS (Table 2 and Supplementary Table 1), of which 4 have not been reported in an NAFL GWAS (in/near PNPLA2, TOR1B, APOH and GUSB). The variants associated with the combined PDFF and NAFL at genome-wide significance (GWS) were nominally significant for both phenotypes (P < 0.05). Furthermore, their effects on the two phenotypes were comparable (Fig. 2), suggesting that variants that increase PDFF also increase risk of NAFL and vice versa. We therefore refer to the combined phenotype as NAFL. The strongest association was with missense variant p.Ile148Met in PNPLA3 (P = 3.0 × 10 −217 , effect = 0.28 s.d. for PDFF and P = 9.7 × 10 −116 , odds ratio (OR) = 1.47 for NAFL). This variant has been reported as being associated with NAFLD 17,29 . Two of the eighteen variants had greater effect on PDFF in men than in women, p.Ile43Val in GPAM (effect in men = 0.09 s.d., effect in women = 0.04 s.d., P males versus females = 0.00099) and p.Glu167Lys in TM6SF2 (effect in men = 0.40 s.d., effect in women = 0.28 s.d., P males versus females = 7.3 × 10 −5 ).

Association of NAFL variants with other traits and outcomes
We tested the 20 NAFL and cirrhosis variants for association with 51 phenotypes related to liver function or NAFLD, using data from Iceland, the UK, Denmark (CHB-CVDC/DBDS) and available meta-analyses (Supplementary Table 4) [42][43][44][45] . The risk alleles of all variants were associated with increased levels of at least one liver enzyme (Figs. 2 and 3). Most variants were associated with cholesterol and sex hormone-binding globulin measures (SHBG), but with inconsistent direction of effects with regard to the NAFL risk allele (Fig. 3). The missense variants p.Thr165Ala in MTARC1 and p.Ile43Val in GPAM had the greatest similarity in their associations with the tested traits ( Figs. 3 and 4). The missense variant p.Cys325Gly in APOH was also associated with increased risk of atrial fibrillation (P = 3.5 × 10 −9 , OR = 1.12, n cases = 96,018), heart failure (P = 2.4 × 10 −9 , OR = 1.12, n cases = 99,214) and higher levels of lipoprotein (a) (Lp(a); P = 2.2 × 10 −56 , effect = 0.10 s.d.) (Supplementary Table 4). The p.His48Arg in ADH1B was associated with reduced risk of alcohol dependence (P = 2.6 × 10 −64 , OR = 0.31, n cases = 60,800). We compared variant associations with alcoholic liver disease (ALD) diagnoses and NAFL diagnoses ( Supplementary Fig. 4). The p.His48Arg was the only variant with a significantly greater effect on the risk of ALD (n cases = 3,818) than of NAFL (OR = 0.33 and OR = 0.85, respectively, P = 6.6 × 10 −10 ).

Rare loss-of-function variants in GPAM and MTARC1
Among the 20 variants with GWS associations with NAFL and cirrhosis, 14 are common or low-frequency missense variants (MAF < 3%). Although these variants implicate probable relevant genes, it is not clear whether loss or gain of function of the encoded proteins reduces or increases disease risk. To investigate this, we looked for associations with rare predicted loss-of-function (pLOF) variants in candidate causal genes at these loci, using data from Iceland and UK. We tested 47 pLOF variants for associations with the same traits that we detected with the lead variant at the locus and found two rare pLOF variant associations in Iceland: p.Arg305Ter in MTARC1 (MAF = 0.34%) and p.Thr189GlyfsTer5 in GPAM (MAF = 0.11%) (Supplementary Table 5). Both genes encode mitochondrial enzymes [46][47][48] . The rare pLOF variant p.Arg305Ter in MTARC1 was associated with lower total cholesterol levels (P = 3.2 × 10 −5 , effect = −0.18 s.d.) (Fig. 4 and Supplementary  Table 5), consistent with the protective allele of the common missense variant p.Thr165Ala in MTARC1. As p.Arg305Ter is predicted to reduce the function of MTARC1, this suggests that p.Thr165Ala also lowers cholesterol levels through reduced MTARC1 function. Therefore, its association with a reduced risk of NAFL is most probably driven by  reduced MTARC1 function. Our findings align with reports of another protein-truncating variant in MTARC1, p.Arg200Ter, that is associated nominally with reduced cirrhosis risk 20 . The minor allele of the common missense variant in GPAM, p.Ile43Val, is associated with increased risk of NAFL and higher total cholesterol levels, whereas the rare pLOF variant, p.Thr189GlyfsTer5, is associated with lower total cholesterol levels (P = 3.5 × 10 −8 , effect = −0.41 s.d.), suggesting that p.Thr189GlyfsTer5 decreases NAFL risk through decreased GPAM function. Although neither of these pLOFs is associated with NAFL or cirrhosis, they give information on whether previously reported associations are gain or loss of function. MTARC1's p.Arg305Ter is located in the gene's last exon, and we found no evidence of nonsense-mediated decay in our RNA data ( Supplementary Fig. 3). We also found a pLOF variant in GCKR, p.Arg540Ter (MAF = 0.60%), which is associated with increased triglyceride levels (P = 6.9 × 10 −7 , effect = 0.16 s.d.), similar to the NAFL risk-increasing allele of the GCKR common missense variant p.Leu446Pro.
Among participants in the UKB, 2,795 had 2 liver MR images taken 2-3 years apart. Consecutive PDFF measurement were correlated (r 2 = 0.90), as were BMI measurements (r 2 = 0.60) (Supplementary Figs. 8 and 9). Changes in PDFF and BMI associated strongly (P = 2.8 × 10 −104 ), with direction of change in the two measures that on average were the same ( Supplementary Fig. 10). We found a nominally significant interaction between p.Ile148Met in PNPLA3 and changes in BMI (P = 0.030; Supplementary Fig. 11). Furthermore, the variations in the PDFF changes were greater among p.Ile148Met carriers than noncarriers (P = 0.0055; Supplementary Fig. 12) and were not fully explained by the interaction between the genotype and change in BMI (P adjusted for change in BMI: P adj = 0.032).

Integration with circulating protein levels
To gain insight into proteins affecting the pathogenesis of NAFL or cirrhosis and search for relevant biomarkers, we analyzed protein levels measured with 4,907 aptamers by SomaScan 49 v.4 in 35,559 Icelanders 50 and 1,459 immunoassays using the Olink Explore 1536 in 47,151 European-ancestry participants from the UKB. The levels of 2,741 proteins associated (after Bonferroni's adjustment) with NAFL (n Iceland = 234 cases, n UK = 572 cases) and 948 with cirrhosis (n Iceland = 111 cases, n UK = 303 cases) in either Iceland or UKB.
We looked for associations between the 20 NAFL and cirrhosis variants and protein levels. Sixteen variants were associated with the level of one or more proteins in trans, using either Iceland (SomaScan) or UKB (Olink), as top protein QTLs (pQTLs) (Supplementary Table 7). The trans pQTLs in GCKR and SERPINA1 are nonspecific in that they are associated with 396 and 172 proteins, respectively. Focusing on the variants associated with fewer than 100 proteins, we identified 273 proteins that associated with the variants in trans (that is, variants at  Fig. 15) and eight out of these ten proteins are associated with NAFL or cirrhosis (Supplementary Tables 8-10).  Article https://doi.org/10.1038/s41588-022-01199-5  The missense variant p.Ile148Met in PNPLA3 is a top trans pQTL for 86 proteins, with aldo-keto reductase family 7 member 3 protein (AKR7A3) showing the most significant association in Iceland (SomaScan) (P = 3.1 × 10 −12 , effect = 0.07 s.d.) and keratin 18 (KRT18) in UKB (Olink) (P = 8.4 × 10 −38 , effect = 0.10 s.d.). The rs72613567 [TA] in HSD17B13 was associated most significantly with levels of the DnaJ B member homolog subfamily 9 protein (DNAJB9) in Iceland (SomaScan) and carboxylesterase 3 (CES3) in UKB (Olink). Missense variants in MTARC1 (p.Thr165Ala) and GPAM (p.Ile43Val) were both associated with levels of group XIIB secretory phospholipase A2-like protein (PLA2GXIIB) as measured in Iceland (SomaScan), and inactive Cα-formylglycine-generating enzyme (SUMF2), hypoxia-upregulated protein 1 (HYOU1) and acid sphingomyelinase-like phosphodiesterase 3a (SMPDL3A), as measured in the UKB (Olink). PLA2GXIIB is highly expressed in the liver and the NAFL protective allele of these variants was associated with lower protein levels (P = 9.3 × 10 −12 , effect = −0.062 s.d. and P = 2.4 × 10 −14 , effect = −0.072 s.d., for MTARC1 and GPAM variants, respectively). Moreover, the rare pLOF variants in Iceland in GPAM and MTARC1 were both associated with lower PLA2GXIIB protein levels (P = 7.4 × 10 −4 , effect = −0.44 s.d. and P = 0.013, effect = −0.19 s.d., respectively). We note that other NAFL variants were associated nominally with PLA2GXIIB, some with an opposite direction of effects to the LOF variants in GPAM and MTARC1 (Supplementary Fig. 13). Therefore, to explore whether any protein measured in plasma may have a causal role in disease, we performed Mendelian randomization using all-trans pQTLs for each protein. The effects of pQTLs of the transferrin receptor protein 1 (TFRC) were proportional to their effect on PDFF (P ivw-olink = 5.6 × 10 −10 ,  Table 11), suggesting that TFRC may have a causal role in NAFL. Apart from TFRC, the analysis suggests that the alterations in protein levels in plasma are not causal but rather a consequence of disease because, for many proteins, the effects of the set of NAFL variants on PDFF were proportional to their effects on protein level (Supplementary Table 12).
We performed a pathway analysis using the PANTHER v. 16.0 tool to seek understanding of the variant-liver disease associations. The 273 proteins associating with NAFL variants were enriched for multiple metabolic and catabolic processes, including the metabolism of hormones, lipids, alcohol, vitamins, steroids and monocarboxylic acid among other pathways and biological processes (Supplementary  Table 13).
To investigate whether plasma proteins can effectively discriminate between having an NAFL and cirrhosis, we trained penalized logistic regression models using liver enzymes, age, sex and BMI as a baseline, as well as using the plasma proteome and genetic risk scores (GRSs). In both Iceland (SomaScan) and UKB (Olink), the models trained with the plasma proteome outperformed other models in discriminating between NAFL and cirrhosis, NAFL and the population, and cirrhosis and the population (Fig. 5 and Supplementary  Tables 14-17).

Discussion
We identified 18 sequence variants that are associated with NAFL and two additional variants that are associated with cirrhosis and not PDFF.
The NAFL variants affect the risk of cirrhosis and HCC proportional to their effects on PDFF, supporting a causal relationship between hepatic fat accumulation and these diseases 51,52 .
The NAFLD variants are associated with other phenotypes, with variable patterns of association. One interpretation of this is that perturbations of more than one biochemical pathway lead to NAFLD. The strongest NAFL-associated variant, p.Ile148Met, in PNPLA3 interacts with rs72613567[TA] in HSD17B13 and BMI but not with age. Longitudinal PDFF measures suggest that p.Ile148Met carriers are more susceptible to change in BMI than noncarriers, in line with previous studies suggesting that the p.Ile148Met genotype increases sensitivity to the beneficial effects of dietary interventions and bariatric surgery [53][54][55][56] .
Sixteen NAFL-associated variants were coding, five were eQTLs and three sQTLs, and many implicated genes were involved in lipid metabolism, reinforcing the notion of its importance as a primary process in NAFL pathogenesis 57,58 . For example, variants at TM6SF2 exert their effects on hepatic lipid concentration by reducing TM6SF2 function, causing reduced secretion of triglyceride-rich lipoproteins 59 . Several variants affect NAFL and blood lipids in the same direction, including p.Leu446Pro in GCKR, which also decreases glycated hemoglobin and type 2 diabetes risk. This is consistent with evidence that the variant increases hepatic glucose metabolism and de novo lipogenesis 60 . The p.His48Arg in ADH1B is associated with less risk of having an NAFL diagnosis. The variant associates with less alcohol consumption and reduces the risk of ALD substantially more than the risk of NAFL. The amount of consumed alcohol (2 and 3 units of alcohol for women and men, respectively) used to distinguish between NAFL and ALD is quite arbitrary 61 . Therefore, it is likely that the association of the variant with NAFL diagnosis is driven by its effect on alcohol consumption below the ALD diagnosis threshold. None of the other NAFL variants was associated with a significantly greater effect on ALD than NAFL. The NAFL associations include missense variants in GPAM and MTARC1, which both encode mitochondrial enzymes and are highly expressed in liver and adipose tissue 31 . GPAM encodes glycerol-3-phosphate acyltransferase 1, which catalyzes the first step of triglyceride synthesis. Gpam knockout mice have lower hepatic triglyceride content and overexpression has the opposite effect 48,62 .
The results for MTARC1 are in line with previous reports of the effects of missense and LOF variants in MTARC1, but the mechanism explaining the associations was less clear 20 . The similarity in their associations with other traits indicates that mutual or similar pathways explain the MTARC1 and GPAM associations. Furthermore, both missense variants affect plasma levels of PLA2GXIIB protein. We identified rare pLOF variants in both genes affecting the same traits and reducing liver enzymes and total cholesterol levels. This suggests that GPAM and MTARC1 have an etiological role and inhibition of GPAM and MTARC1 could be therapeutic for NAFL or NASH, with a favorable effect on the metabolic profile. The lack of associations with increased risk of a large set of diseases is reassuring with regard to treatment safety.
Among associations with NAFL are missense variants in APOH and GUSB. APOH is highly and almost exclusively expressed in the liver 31 . The other trait associations of p.Cys325Gly in APOH strongly suggest a role in lipid metabolism 35,63,64 . Furthermore, the variant has been reported in GWASs of coronary artery disease 65 and Lp(a) levels 64 . We replicated the Lp(a) association and observed associations with an increased risk of atrial fibrillation and heart failure. Our results indicate a role for GUSB in the etiology of NAFL. The missense variant p.Leu649Pro in GUSB associates with both NAFL and RNA expression levels of GUSB. Furthermore, seven NAFL variants associate with GUSB plasma protein levels. GUSB encodes β-glucuronidase, a lysosomal enzyme involved in the breakdown of glycosaminoglycans 66 .
Plasma proteome analysis revealed that missense and pLOF variants in MTARC1 and GPAM are associated with levels of PLA2GXIIB in plasma. PLA2GXIIB is highly expressed in the liver and knockout mice have severe hepatosteatosis 67 . However, our proteomic analysis does not support an etiological role for PLA2GXIIB because NAFL variants affect plasma PLA2GXIIB levels with effect directions inconsistent with their NAFL effects. Evidence suggests that PLA2GXIIB is a mediator of hepatic lipoprotein secretion and its inhibition results in decreased levels of plasma lipids 67,68 . Thus, PLA2GXIIB may mediate variant effects on circulating lipids without directly affecting hepatic fat and thus could serve as a biomarker of drug target engagement.
Diagnosis and monitoring of complications in patients with NAFLD are challenging 2 . We designed models including plasma proteins that outperformed a model trained on liver enzymes and GRSs in discriminating between NAFL and cirrhosis diagnoses. Thus, levels of plasma proteins have the potential to serve as a noninvasive tool for use in the diagnosis and monitoring of disease, whereas GRSs are associated with a lifetime risk of disease. THBS2 was elevated in individuals with cirrhosis compared with NAFL and the population, and ACY1 in individuals with NAFL compared with the general population. Intrahepatic THBS2 expression levels have previously been shown to correlate with fibrosis in patients with NAFLD 69 . The association of IGFBP2 and IGFBP7 with cirrhosis and NAFL is consistent with previous studies of NASH progression 70 . Both proteins bind insulin-like growth factors (IGFs) and modulate their availability 71,72 . IGFs are mainly produced in the liver, and IGFBP2 and IGFBP7 elevation may reflect imbalances in the IGF system due to liver damage 73 . However, an etiological role has been suggested for IGFBP7, which may contribute to hepatic fibrogenesis 74 and act as a tumor suppressor in HCC 75 . Levels of SHBG are also associated with cirrhosis compared with NAFL, consistent with previous reports of a positive correlation with advanced fibrosis in NASH 70 . There are conflicting epidemiological studies about whether NAFLD is associated with increased or decreased levels of SHBG. In line with this, many NAFL variants are associated with SHBG plasma levels with inconsistent directions of effect compared with their effect on hepatic fat content 76 .
A limitation of the present study is the lack of data enabling a more detailed phenotype stratification, in particular histological data for disease staging. Furthermore, information on other causes of liver disease, such as alcohol consumption, is limited. Our approach is, however, in line with the recent opinion to not base the disease diagnosis on a criterion of exclusion of other diseases, such as ALD 77 .
In conclusion, we used multiomics data to shed light on the genetic basis of NAFLD. Analysis of pLOF variants, blood RNA expression and plasma proteomics data pointed to causative genes and whether changes in their functions contribute to the pathogenesis. We demonstrated the diverse effects of NAFL risk alleles on other diseases and traits, including blood lipids and proteins, and showed that plasma proteomics has the potential to stage NAFLD.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41588-022-01199-5.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. Article https://doi.org/10.1038/s41588-022-01199-5

Study populations
The UKB project is a large prospective cohort study of ~500,000 individuals from across the United Kingdom, aged between 40 and 69 years at recruitment 78 . The present study has collected extensive phenotypic and genotypic information on participants, including ICD-10-coded diagnoses from inpatient and outpatient hospital episodes and abdominal MRI through its imaging study 79 . Genotype imputation data were available for 487,409 participants (May 2017 release), of which 408,658 were included because they self-reported as white 80 . The UKB resource was used under application no. 56270. All phenotype and genotype data were collected following an informed consent obtained from all participants. The North West Research Ethics Committee reviewed and approved UKB's scientific protocol and operational procedures (REC reference no.: 06/MRE08/65).
The Icelandic deCODE genetics study is based on WGS data from 49,708 Icelanders participating in various research projects at deCODE genetics. Variants identified through WGS were imputed into 155,250 chip-genotyped Icelanders using long-range phasing 81 and their untyped close relatives based on genealogy. Sequencing was done using GAIIx, HiSeq, HiSeqX and NovaSeq Illumina technology. SNPs and insertions/deletions (indels) were identified, and their genotypes were called using joint calling with Graphtyper 82-84 . All participants who donated blood signed an informed consent. The personal identities of the participants and biological samples were encrypted by a third-party system. The study was approved by the Icelandic Data Protection Authority and the National Bioethics Committee of Iceland (no. VSN-20-182).
FinnGen summary statistics, including fatty liver disease and cirrhosis, were imported in December 2020 from a source available to researchers (v.4: https://www.finngen.fi/en/access_results) and methods have been documented (https://finngen.gitbook.io/documentation). The FinnGen database consists of samples collected from the Finnish biobanks and phenotype data collected at the national health registers. The Coordinating Ethics Committee of the Helsinki and Uusimaa Hospital District evaluated and approved the FinnGen research project. The project complies with existing legislation (in particular, the Biobank Law and the Personal Data Act). The official data controller of the present study is the University of Helsinki.
The Copenhagen Hospital Biobank Cardiovascular Study (CHB-CVDC) was used to acquire secondary cardiovascular phenotypes. CHB-CVDC involves a targeted selection of patients with cardiovascular disease from the CHB, a biobank based on patient blood samples drawn in Danish hospitals 43 . The CHB-CVDC has been approved by the National Committee on Health Research Ethics (1708829) and the Danish Data Protection Agency (P-2019-93). For binary phenotypes, the control group included blood donors from the Danish Blood Donor Study (DBDS) (n = 99,000), approved by the Danish Data Protection Agency (P-2019-99) and the Scientific Ethical Committee system (NVC 1700407) 44 . Chip typing and genotype imputation of CHB-CVDC and DBDS were performed at deCODE genetics using a north European sequencing panel of 15,576 individuals (including 8,429 Danes).
The samples from the USA (Intermountain dataset) were whole-genome studied using NovaSeq Illumina technology (n = 8,288) and genotyped using Illumina GSA chips (n = 28,279). Samples were filtered on a 98% variant yield. Over 245 million high-quality sequence variants and indels, to a mean depth of at least 20×, were identified using Graphtyper 82 . Quality-controlled chip genotype data were phased using Shapeit 4 (ref. 85 ). A phased haplotype reference panel was prepared from the sequence variants using the long-range phased chip genotype data. The haplotype reference panel variants were then imputed into the chip-genotyped samples using inhouse tools and methods described previously 83,84 . In the US association analysis, samples assigned <93% CEU ancestry (northern European from Utah) were excluded. We adjusted for sex, age and the first 20 principal components. Phenotypic data were based on ICD-10 code diagnoses of individuals. The Intermountain Healthcare Institutional Review Board approved the study and all participants provided written informed consent before enrollment.

Imaging protocol
The MR images used for calculating PDFF for 36,116 individuals were collected as a part of the UKB abdominal protocol, which, in turn, is part of the UKB imaging enhancement 79 . Two acquisitions were used, a single-slice GRE sequence and a single-slice IDEAL sequence 86 . The slice was captured through the center of the liver superior to the porta hepatis. The GRE sequence was captured using the following settings: repetition time (TR) = 27 ms, time to echo (TE) = [2.38, 4.76, 7.15, 9.53, 11.91, 14.29, 16.67, 19.

Imaging analysis
We used two different approaches for calculating the PDFF from the liver MR images depending on whether the acquisition was GRE (n = 8,448) or IDEAL (n = 27,668). We implemented the PDFF estimation methods using a tailored Python code. For the GRE acquisition, we used a three-point Dixon method 87 to compute a PDFF map using the second, fourth and sixth echoes. Eight 25-voxel rectangular regions of interest (ROIs) were defined within the liver and we computed the mean and s.d. of the PDFF maps over those ROIs. The reported PDFF was the ROI with the lowest s.d. By choosing the lowest s.d., we avoided ROIs with water-fat swaps. For the IDEAL acquisition, we assumed the following signal model 88 for each voxel: where ρ w and ρ f are the water and fat components, respectively, Δf is the chemical shift of fat with respect to water, φ quantifies B 0 field inhomogeneity, R * 2 is an MRI relaxation constant and t i the ith echo time. The parameters of interest, ρ w , ρ f and R * 2 , are estimated from the signal model using an iterative weighted least-squares algorithm. The PDFF map was finally constructed using |ρ f | /(|ρ w | + |ρ f |) at each voxel. The reported PDFF was calculated in the same way as for the GRE case. Iron concentration (Supplementary Fig. 1) can be estimated by using the R * 2 coefficient 89 : Iron concentration = 0.202 + 0.0254 R * 2 . To evaluate the correspondence between the PDFF scores for the IDEAL and the GRE acquisition, we investigated 1,222 PDFF scores computed for both. We also fit a linear model to quantify this relationship yielding in units of percentage: and R 2 = 0.92. There are 3,869 PDFF scores for the GRE acquisition that we computed and that are available in the UKB data showcase. Supplementary Fig. 2 shows a scatterplot of those two sets of scores demonstrating good agreement. The relationship is given in units of percentage: PDFF GRE−UKB = −0.5935 + 1.0152PDFF GRE and R 2 = 0.96.

Secondary phenotypes
NAFL-associated variants were tested for association with other phenotypes from deCODE genetics, UKB, FinnGen, CHB-CVDC/DBDS or other publicly available data sources, which contain extensive medical information on various traits. A total of 51 phenotypes was chosen for the analysis (Supplementary Table 4) based on known associations with NAFLD, liver function or the identified variants. These included other liver diseases and known blood markers of liver function, blood lipids, cardiovascular diseases, diabetes, anthropometric traits, hematological traits and hormone levels.

Association testing and meta-analysis
We used logistic regression to test for association between sequence variants and binary phenotypes assuming an additive genetic model. For deCODE genetics, the model included the following covariates: sex, county of birth, current age or age at death (first-and second-order terms included), blood sample availability for the individual and an indicator function for the overlap of the lifetime of the individual with the time span of the phenotype collection. In CHB-CVDC/DBDS, the covariates were sex, age and 20 principal components to adjust for population stratification and blood sample availability. In the UKB study, 40 principal components were used to adjust for population stratification, and age and sex were included as covariates in the logistic regression model. When analyzing PDFF, BMI was included as a covariate in the analysis to increase power of associations. We used a linear mixed model implemented in BOLT-LMM to test for association between sequence variants and quantitative traits. The measurements used were adjusted for sex, year of birth and age at measurement (when available), and these were subsequently standardized to have a normal distribution. For the meta-analysis of summary-level statistics from different populations, we used a fixed-effects inverse variance method based on effect estimates and s.e.m. With the aim of studying NAFL, we combined summary-level data from the GWAS of PDFF and the meta-analysis of GWAS using ICD-10-code-based NAFL with multitrait analysis of genome-wide association summary statistics (MTAG) 28 . To have two nonoverlapping sets, we excluded individuals from the NAFL ICD-10 code analysis who had a PDFF measurement in the UK. To account for inflation in test statistics due to cryptic relatedness and stratification, we applied the method of LD score regression 90 . For the GWS, we accounted for multiple testing with a weighted Bonferroni's adjustment, using as weights the enrichment of variant classes with predicted functional impact among association signals estimated from the Icelandic data 91 . This yielded significance thresholds of 1.8 × 10 −7 for variants with high impact (including stop-gained and loss, frameshift, splice acceptor or donor and initiator codon variants), 3.5 × 10 −8 for variants with moderate impact (missense, splice-region variants and in-frame indels), 3.2 × 10 −9 for low-impact variants (including synonymous, 3′-and 5′-UTRs and upstream and downstream variants), 1.6 × 10 −9 for other variants in DNase I hypersensitivity sites and 5.3 × 10 −10 for all other variants.

Variant-variant interaction with liver enzymes
To test whether a specific primary variant is affecting the mean effect of another secondary variant in a given quantitative trait (that is, to test if they are interacting), we split the population into three groups based on the genotype of the primary variant, denoted by g p ∈ {0, 1, 2}. We then estimated the mean effect, β 0 , β 1 and β 2 of the secondary variant in each group separately, where the quantitative trait was standardized to a normal distribution. We estimated the interaction between the primary and secondary variants by considering the following model: where b and γ are the unknown parameters to be estimated. To assess the significance of the interaction parameter γ, we applied a likelihood ratio test, comparing our model with the null model: β 0 = β 1 = β 2 (or equivalently γ = 0).

RNA-seq analysis
RNA-seq analysis was performed on whole blood (n = 17,846) and subcutaneous adipose tissue (n = 750). We isolated RNA using Chemagic Total RNA Kit special (PerkinElmer) in whole blood and RNAzol RT in adipose tissue, according to the manufacturer's protocol (Molecular Research Center, RN 190). The concentration and quality of the RNA were determined using an Agilent 2100 Bioanalyzer (Agilent Technologies). RNA was prepared and sequenced on the Illumina HiSeq 2500 and Illumina Novaseq systems according to the manufacturer's recommendation.
RNA-seq reads were aligned to personalized genomes using the STAR software package v.2.5.3 with Ensembl v.87 gene annotations 92,93 . Gene expression was computed based on personalized transcript abundances estimated using kallisto 94 . Association between sequence variants and gene expression was tested using BOLT-LMM, assuming an additive genetic model and quantile-normalized gene-expression estimates, adjusting for measurements of sequencing artefacts and demography variables. The strongest association within 1 Mb of each gene with P < 1 × 10 −7 was called a top cis-eQTL.
Quantification of alternative RNA splicing in whole blood was done using LeafCutter 95 . The cis association between sequence variants and quantified splicing (cis-sQTL) was tested using linear regression assuming an additive genetic model and quantile-normalized percentage-spliced-in values (PSI) of each splice junction, adjusting for measurements of sequencing artefacts, demography variables, and 15 leave-one-chromosome-out principal components of the quantile-normalized PSI matrix. All variants with MAF > 0.2% within 30 kb of each LeafCutter cluster were tested, and the strongest association for each splice junction with P < 1 × 10 −8 was called a top cis-sQTL.

SomaScan proteomic analysis
Plasma samples were collected from Icelanders through two main projects: the Icelandic Cancer Project (52% of participants; samples collected from 2001 to 2005) and various genetic programs at deCODE genetics, Reykjavík, Iceland (48%). The samples collected at deCODE genetics were mainly collected through the population-based deCODE Health study. The average participant age was 55 years (s.d. = 17 years) and 57% were women. In the case of repeated samples for an individual, we selected one of them at random. This left measurements for 39,155 individuals. Of these, 35,559 Icelanders were used in the protein GWASs, because they also had genotype information 50 . All plasma samples were measured with the SomaScan v.4 assay (Soma-Logic, Inc.) containing 4,907 aptamers, providing a measurement of relative binding of the plasma sample to each of the aptamers in relative fluorescence units (r.f.u.). When testing for association between protein levels and disease, logistic regression was used with age and sex as covariates. The date of diagnosis was not available and the analysis was therefore not adjusted for the time from diagnosis. A pQTL association is considered to be in cis if the associated variant is located no more than 1 Mb from the transcription start site of the gene that encodes the measured protein, and in trans otherwise. A pQTL was considered to be significant if identified in the previous study (P < 1.8 × 10 −9 ). A top pQTL is the top (most significantly associated) variant per megabase bin. Article https://doi.org/10.1038/s41588-022-01199-5

Olink proteomic analysis
For a subset of 54,265 individuals in the UKB study (47,151 of British and Irish ancestry), the levels of 1,472 proteins were measured with the Olink Explore 1536 platform as a part of the UKB-Pharma Proteomics Project (UKB application no. 65851). A large majority of the samples were randomly selected across the UKB. The UKB plasma samples were measured using the Olink Explore 1536 platform (https://www. olink.com/products-services/explore) at Olink's facilities in Uppsala, Sweden. Measurements with the Olink Explore platform use the NPX values recommended by the manufacturer, which include normalization. When testing for associations with sequence variants, the protein levels were standardized to a normal distribution.

Mendelian randomization analysis
We performed Mendelian randomization analyses with the Mendelian-Randomization R package, using both inverse variance, weighted linear regression 96 and Egger regression 97 , with default settings.

Pathway analysis
The trans pQTLs that were associated with the identified PDFF variants were grouped together and tested for enrichment in Reactome pathways, Gene ontology terms (biological process, molecular function and cellular component) and PANTHER protein classes. The analysis was performed with the PANTHER v.16.0 tool 98 using Fisher's exact test and false discovery rate (FDR) correction. The associated pQTLs were tested against a reference list that included all measured SomaAmers (n = 5,286). Pathway analysis was performed for all variants associating with more than two proteins.

Liver disease stage classification using plasma proteomics
To examine how effective circulating protein levels (SomaScan and Olink panels) are at discriminating liver disease stages, we trained a logistic regression model with an elastic net penalty on protein levels to classify between individuals diagnosed with cirrhosis and NAFL. To test how this model compares with one trained on liver enzymes, an additional logistic regression model was exclusively trained on age, sex, BMI, ALT, AST and γ-glutamyl transferase. Genetic risk scores for NAFL and cirrhosis from the identified genetic variants were constructed to compare with the protein-based prediction. We created NAFL and cirrhosis GRSs by calculating the sum-of-effect alleles of the NAFLD variants weighted by their cirrhosis and NAFL GWAS effects. In addition, logistic regression models combining the protein scores, liver enzymes and GRSs were trained. In that case, the protein scores were calculated using out-of-fold predictions (with tenfold). The SomaScan analysis was performed on 181 individuals with NAFL and 73 with cirrhosis, and the Olink analysis was performed on 610 individuals with NAFL and 262 with cirrhosis. Models trained to discriminate between the presence of disease diagnosis and population were trained, respectively, on an additional set of 20,619 individuals (SomaScan) and 38,018 individuals (Olink) without cirrhosis and NAFL. The SomaScan measurements were log(transformed) to reduce the effect of outliers. However, we found that this was not necessary for the Olink data. Boruta 99 was used to select relevant features. The selected proteins were preprocessed with Yeo-Johnson power transforms 100 and then scaled to zero mean and unit variance (estimated from the training set) before being fed to the logistic regression model. The elastic net λ and α parameters were tuned using grid search to minimize tenfold crossvalidated average precision. Model performance was evaluated by considering the mean and s.e.m. of the receiver operating characteristic (ROC) area under the curve (AUC) of 1,000 repeated, tenfold-stratified crossvalidation runs.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

March 2021
Data Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

Life sciences study design
All studies must disclose on these points even when the disclosure is negative.

Sample size
Sample sizes are reported in the article and correspond to all available data Data exclusions No available data was excluded from the study

Replication
The NAFLD GWAS analysis was performed using data from 4 populations (Iceland, UK, USA and Finland) and results across populations were compared.
Randomization Not applicatble (GWAS study, not a randomized trial)

Blinding
Not applicable (GWAS study, not a randomized trial, so no blinding is required) Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. 155,250 chip-genotyped Icelanders as well as their untyped close relatives based on genealogy. Finngen summary statistics, including fatty liver disease and cirrhosis, were imported on December 2020 from a source available to researchers (version 4; https://www.finngen.fi/en/access_results) and methods have been documented (https:// finngen.gitbook.io/documentation/). The Copenhagen Hospital Biobank Cardiovascular Study (CHB-CVDC) was used to acquire secondary cardiovascular phenotypes. CHB-CVDC involves a targeted selection of patients with cardiovascular disease from the CHB, a biobank based on patient blood samples drawn in Danish hospitals. For binary phenotypes, the control group included blood donors from The Danish Blood Donor Study (DBDS). The samples from the US (Intermountain dataset) were WGS using NovaSeq Illumina technology and genotyped using Illumina GSA chips.

Recruitment
For the deCODE Genetics study individuals were recruited through various research projects at deCODE genetics. The participants are a large fraction of the adult Icelandic population. UK Biobank holds data on half a million participants throughout the UK. All participants in UK Biobank were recruited through assessment centres, designed specifically for this purpose. The US data are individuals recruited at the Intermountain healthcare. The FinnGen database consists of samples collected from Finnish biobanks. CHB-CVDC involves a targeted selection of patients with cardiovascular disease from the CHB, a biobank based on patient blood samples drawn in Danish hospitals.

Ethics oversight
All participating subjects in the deCODE genetics study who donated blood signed informed consent. The personal identities of the participants and biological samples were encrypted by a third-party system. For the Intermountain dataset, the Intermountain Healthcare Institutional Review Board approved the study, and all participants provided written informed consent prior to enrollment.
Note that full information on the approval of the study protocol must also be provided in the manuscript.