Abstract
Nonalcoholic fatty liver disease (NAFLD) is a prevalent, heritable trait that can progress to cancer and liver failure. Using our recently developed proxy definition for NAFLD based on chronic liver enzyme elevation without other causes of liver disease or alcohol misuse, we performed a multi-ancestry genome-wide association study in the Million Veteran Program with 90,408 NAFLD cases and 128,187 controls. Seventy-seven loci exceeded genome-wide significance of which 70 were novel, with an additional European-American specific and two African-American specific loci. Twelve of these loci were also significantly associated with quantitative hepatic fat on radiological imaging (n=44,289). Gene prioritization based on coding annotations, gene expression from GTEx, and functional genomic annotation identified candidate genes at 97% of loci. At eight loci, the allele associated with lower gene expression in liver was also associated with reduced risk of NAFLD, suggesting potential therapeutic relevance. Functional genomic annotation and gene-set enrichment demonstrated that associated loci were relevant to liver biology. We expand the catalog of genes influencing NAFLD, and provide a novel resource to understand its disease initiation, progression and therapy.
Introduction
Chronic liver disease is a major contributor to global morbidity and mortality, with complications of cirrhosis and hepatocellular carcinoma1. In particular, nonalcoholic fatty liver disease (NAFLD) is an increasingly common cause of chronic liver disease with an estimated world prevalence of 25% among adults and associated metabolic risk factors1-5. In the United States (US), NAFLD prevalence is projected to reach 33.5% among adult population by 2030, due in large part to the rising obesity and associated metabolic disorders6. NAFLD is defined by ≥5% fat accumulation in the liver (hepatic steatosis) in the absence of other known causes for liver disease, based on liver biopsy and/or non-invasive radiological imaging3,4. The clinical spectrum of NAFLD ranges from benign steatosis to nonalcoholic steatohepatitis (NASH) involving inflammation and hepatocellular injury with fibrosis progression. At least 20% of patients with NAFLD develop NASH with increased risk of consequent cirrhosis and liver cancer5,6. To date, there is no licensed drug approved to treat NAFLD and prevent its progression.
Individual susceptibility to NAFLD involves both genetic and environmental factors. Risk factors for NAFLD include obesity (in particular, abdominal adiposity), insulin resistance and features of metabolic syndrome2,5-7, with current estimates of NAFLD heritability ranging between 20% to 50%8. Several genetic variants that promote the full spectrum of fatty liver disease have been identified in genome-wide association studies (GWAS) utilizing cohorts based on liver biopsy, imaging, and/or isolated liver enzyme values9-22. The most prominent of these include p.I148M in PNPLA3 and p.E167K in TM6SF2, which increase NAFLD risk, and a loss-of-function variant in HSD17B13 that confers protection against NASH16. However, the limited number of genetic associations in NAFLD contrasts with other cardiometabolic disorders where hundreds of loci have been mapped to date, traits that include obesity23,24, type 2 diabetes25 and plasma lipids26. This also highlights the need for expanded discovery based on larger sample size and population diversity, with further integration with existing functional genomics data sets to identify candidate genes from leading, non-coding associations27.
The Million Veteran Program (MVP) is among the world’s largest and ancestrally diverse biobanks28. The availability of comprehensive, longitudinally collected Veterans Health Administration (VA) electronic health records for US Veteran participants in the MVP also makes it a promising resource for precision medicine. As NAFLD is markedly underdiagnosed clinically due to limited access to liver biopsy and variable use of imaging modalities4, we recently developed and validated a proxy phenotype for NAFLD to facilitate case identification in MVP21. The proxy NAFLD phenotype is based on chronically elevated serum alanine aminotransferase (cALT) levels while excluding other conditions that are known to increase liver enzymes (e.g. viral hepatitis, alcohol dependence, autoimmune liver disease and known hereditary liver disease). We applied this cALT-based proxy NAFLD phenotype to the current build of 430,400 genotyped MVP participants with defined ancestry classification29, and identified 90,408 NAFLD cases and 128,187 controls (inclusion/exclusion criteria for the remaining samples and study design described in Supplementary Figure 1 and Figure 1). We performed a primary GWAS and identified 77 trans-ancestry loci that reached genome-wide significance. We used additional approaches to define NAFLD heritability and genetic correlations with various traits including quantitative hepatic fat measured by liver imaging with computed tomography (CT) and magnetic resonance imaging (MRI), in addition to identifying coding variants in putative causal genes.
Left side of the flow diagram shows our study design with initial inclusion of 430,000 Million Veteran Program participants with genotyping and ancestry classification by HARE, exclusion of individuals with known liver disease or alcohol dependence and inclusion of subjects based on chronic ALT elevation (case) or normal ALT (control). This resulted in 90,408 NAFLD cases and 128,187 controls with EA, AA, HISP and ASN ancestries that were examined in primary trans-ancestry and ancestry-specific genome-wide association scans in discovery (stage 1) and internal replication stages (stage 2) with further meta-analysis. Right side of the flow diagram highlights our results of trans-ancestry and ancestry-specific meta-analyses identifying 77 trans-ancestry loci + 1 EA-specific + 2 AA-specific loci that met genome-wide significance, with additional results of external replications, locus fine-mapping via GCTA, signal fine-mapping (95% credible sets), heritability estimation and genetic correlations by LDSC, physiological categorization of discovered loci based on pleotropic trait associations (mainly liver, metabolic and inflammation), candidate gene nomination and polygenic risk score.
Results
Diverse NAFLD case and control subjects enriched for metabolic disorders in MVP
Our study consisted of 90,408 NAFLD cases and 128,187 controls across two stages and comprising four ancestral groups, namely European Americans (EA, 75.1%), African Americans (AA, 17.1%), Hispanic-Americans (HISP, 6.9%), and Asian-Americans (ASN, 0.9%, Supplemental Table 1) with the overall sample sizes and study design shown in Figure 1 and Supplemental Figure 1. Consistent with the US Veteran population, MVP cases and controls (n = 218,595) were predominantly male (92.3%) with an average age of 64 years at study enrollment (Supplemental Table 1). With the exclusion of other known causes of liver disease in our phenotype definition21, our cohort was enriched for metabolic disorders, with higher prevalence in cases as compared to controls for type 2 diabetes (71% vs. 47%, P<1×10−5), hypertension (73% vs. 60%, P<1×10−5 and dyslipidemia (82% vs. 70%, P<1×10−5).
Identification of novel trans-ancestry and ancestry-specific NAFLD-associated loci in the diverse MVP population
To identify loci associated with NAFLD, we performed ancestry-specific genome-wide scans by meta-analyzing summary statistics derived from each ancestry followed by trans-ancestry meta-analysis combining data across all ancestries and stages (Methods and Figure 1). In the trans-ancestry scan across stages, 77 independent sentinel SNPs exceeded genome-wide significance (P < 5×10−8), of which 70 were novel whereas 7 (namely PNPLA3, TM6SF2, ERLIN1, TNKS [PPP1R3B], MARC1, HSD17B13, and LYPLAL1) had previously reported genome-wide significant associations with NAFLD (within 500kb and/or CEU r2 LD > 0.05; Figure 2 and Supplemental Table 2)9-22. In addition, 55 loci in EAs, 8 loci in AAs, and 3 loci in HISPs, exceeded genome-wide significance (Supplemental Tables 3-5 and Supplemental Figures 2-4). One SNP (rs4940689) reached genome-wide significance in an ancestry-specific analysis of EAs only (Supplemental Table 3), whereas two SNPs (rs144127357; rs2666559) reached genome-wide significance among AAs only (Supplemental Table 4). No loci in ASNs achieved genome-wide significance, likely due to limited sample size in this group (Supplemental Figure 5).
Nominated genes are indicated for 77 loci reaching genome-wide significance (P<5×10−8). Previously reported NAFLD-loci with genome-wide significant association are indicated in green font.
Among the eight AA-specific lead SNPs, three were intronic: rs115038698 in the ABCB4 locus, rs144127357 in TJP2, and rs2666559 in NRXN2. Two of these variants were nearly monomorphic in EA but polymorphic in AA (rs115038698 MAF AA: 1.2%, MAF EA: 0%; rs144127357 MAF AA: 3.14%, MAF EA: <0.001%). In contrast, the tagged variant rs2666559 was common in both AA (MAF = 19.1% in AA, gnomAD AF = 17.2% in Africans) and EA (AF = 69.1% in EA, gnomAD AF = 68.4% in non-Finnish Europeans).
Internal and external replication of NAFLD-associated loci
We next compared the extent of association across both Stage 1 (primary analysis) and Stage 2 replication stage internally in MVP and externally in the Penn Medicine Biobank (PMBB, Methods). Of the 77 associated SNPs from the trans-ancestry meta-analysis, 56 reached genome-wide significance in Stage 1 subset, of which 32 passed Bonferroni significance (0.05/56) in Stage 2 replication in MVP (Supplemental Table 2). All 77 SNPs showed directional concordance in effect estimates between the two stages. External replication for our trans-ancestry lead SNPs in PMBB (n=72 of our loci were genotyped) with 2,570 cases and 3,802 controls demonstrated that 8 out of 72 variants were directionally consistent and nominally associated (signed binomial-test P=4.4×10−4). Furthermore, 4 out of 8 loci discovered in the AA-specific scan (signed bionomial-test P=2.5×10−5) and 1 of 3 loci discovered in the HISP-specific scan (signed bionomial test P=0.07) were also directionally consistent and nominally associated in PMBB (Supplemental Tables 6-9). In summary, we found 73 novel loci associated with NAFLD that were identified by trans and single-ancestry association studies and supported by replication in multiple stages and studies.
Concordance of cALT-based NAFLD loci with CT/MRI-based quantitative hepatic fat
To place our discoveries into physiological context, we next investigated the extent to which the 77 trans-ancestry SNPs from our NAFLD GWAS associated with quantitative measures of hepatic fat, derived from CT/MRI imaging studies (Methods). We performed a trans-ancestry meta-analysis among 44,289 participants in the UK Biobank, PMBB, Framingham Heart Study, University of Maryland Old Order Amish Study, and Multi-Ethnic Study of Atherosclerosis (Supplemental Table 10). We found that 24 SNPs were nominally associated with quantitative hepatic fat (P < 0.05), of which 12 (15.6% of 77 loci) exceeded Bonferroni multi-test correction (P < 6.5×10−4, including PNPLA3, TM6SF2, APOC1;APOE, GPAM, MARC1, KIAA0196 [TRIB1], MTTP, APOH, PIK3R2;IFI30;MPV17L2, TNKS [PPP1R3B], COBLL1;SCN2A and PPARG). Notably, PNPLA3, TM6SF2, and TNKS [PPP1R3B] were previously identified using imaging data11,12,14, and the direction of effect for all significant SNPs was concordant between chronic ALT elevation and hepatic fat, with the known exception of the variant at the PPP1R3B locus12.
Identification of additional independent NAFLD-associated variants by conditional analysis
To discover additional variants independent of our lead NAFLD signals, we next performed approximate conditional analysis using the leading sentinel variants at our 77 trans-ancestry associated loci. We detected a total of 41 conditionally independent SNPs flanking three known (PNPLA3, HSD17B13, and ERLIN1) and 17 novel NAFLD loci in EA (Supplemental Table 11). Nine conditionally independent SNPs were observed at the PNPLA3 locus in MVP, indicative of the complexity of this locus. For one of the novel loci, located on chromosome 12 between 121-122Mb, the trans-ancestry lead variant (rs1626329) was located in P2RX7, whereas the lead peak for EA mapped to HNF1A (rs1169292, Figure 3). Both are strongly linked to distinct coding variants (P2RX7: rs1718119, Ala348Thr and HNF1A: rs1169288, Ile27Leu) and are compelling candidate genes for metabolic liver disease. In AA, we observed eight conditionally independent variants at four genomic loci, one at PNPLA3 and three novel loci (Supplemental Table 12), including four in GPT, two in AKNA, one in ABCB4. In HISP, two conditionally independent variants in the PNPLA3 locus were identified (Supplemental Table 13). Collectively, 51 additional variants were identified at 24 loci across ancestries based on conditional analysis.
The lead variants in each analysis are highlighted. The orange arrow refers to the proxy SNP of rs1626329 in the European-only analysis.
Fine mapping to define potential causal variants in 95% credible sets
To leverage increased sample size and population diversity to improve fine-mapping resolution, we computed 95% credible sets using Wakefield’s approximate Bayes’ factors30 derived from the trans-ancestry meta-regression, EA, AA, and HISP scans (Supplemental Table 14-17, Methods). In a comparison of the trans-ancestry and EA-only scans, the trans-ancestry meta-regression reduced the median 95% credible set size from 9 (IQR 3 - 17) to 7.5 variants (IQR 2 - 13). A total of 11 distinct NAFLD associations were resolved to a single SNP in the trans-ancestry meta-regression, with 4 additional loci suggestive a single SNP in the EA (n=2) and AA (n=2) ancestry-specific meta-analyses that were not already resolved to a single SNP via the trans-ancestry analysis.
Heritability of NAFLD and genetic correlations with other phenotypes
To tabulate trait heritability and genetic correlation with others, we utilized LD score regression31-33 (Methods). Consistent with our discovery of novel genetic associations, we estimated the SNP-based liability-scaled heritability at 16% (95% CI: 12-19, P < 1.0×10−6) in EA. Genome-wide genetic correlations of NAFLD were calculated with a total of 774 complex traits and diseases by comparing allelic effects using LD score regression with the EA-specific NAFLD summary statistics. A total of 116 significant associations were observed (Bonferroni correction for 774 traits P < 6.5×10−5, Supplemental Table 18). Consistent with the previous epidemiological associations with metabolic syndrome traits, we observed strong correlations with cardiometabolic risk factors including measures of obesity and adiposity, type 2 diabetes, hypertension, dyslipidemia, coronary artery disease, family history of metabolic risk factors and general health conditions in addition to educational attainment.
Liver-specific enrichment of NAFLD heritability
To ascertain the tissues contributing to the disease-association underlying NAFLD heritability, we performed tissue-specific analysis using stratified LD score regression. The strongest associations were observed in genomic annotations surveyed in liver, hepatocytes, adipose, and immune cell types among others (e.g., liver histone H3K36me3 and H3K4me1, adipose nuclei H3K27ac, spleen TCRγd, eosinophils in visceral fat; Supplemental Table 19). Medical subject heading (MeSH)-based analysis showed enrichment mainly in hepatocytes and liver (False Discovery Rate (FDR) < 5%, Supplemental Table 20). Gene set analysis showed strongest associations for liver and lipid-related traits (P-value < 1×10−6, Supplemental Table 21). Enrichment analyses using publicly-available epigenomic data (implemented in GREGOR enrichment analysis, Methods) showed that most significant enrichments were observed for active enhancer chromatin state in liver, epigenetic modification of histone H3 in hepatocytes or liver-derived HepG2 cells (e.g. H3K27Ac, H3K9ac, H3K4me1, H3K4me3; Supplemental Table 22 and 23). These analysis support the hypothesis that our GWAS captures multiple physiological mechanisms that contribute heritability to NAFLD. Finally, DEPICT-based predicted gene function nominated gene candidates for 28 genes, including the known genes PNPLA3 and ERLIN1 (FDR <5%, Supplemental Table 24), as well as well-known cardiometabolic disease genes (e.g., PPARG).
Coding variants in putative causal genes driving NAFLD associations
There were six novel trans-ancestry loci for which the lead SNP itself is a coding missense variant (Supplementary Table 25): Thr1412Asn in CPS1 (rs1047891, β=0.037, P=2.8×10−8), Glu430Gln in GPT (rs141505249, β=-2.023, P=9.0×10−62), Val112Phe in TRIM5 (rs11601507, β=0.099, P=1.5×10−14), Ala163Thr in DNAJC22 (rs146774114, β=-0.157, P=2.5×10−8), Glu366Lys in SERPINA1 (rs28929474, β=0.492, P=9.01×10−73) and Cys325Gly in APOH (rs1801689, β=0.17, P=1.5×10−18). To identify additional coding variants that may drive the association between the lead SNPs and NAFLD risk, we investigated predicted loss of function (pLoF) and missense variants strongly linked to the identified NAFLD lead variants (r2 > 0.7, Supplemental Table 25-28). Four previously described missense variants were replicated in the current study, including Thr165Ala in MARC1, Ile292Val in ERLIN1, Glu167Lys in TM6SF2 and Ile148Met in PNPLA3. Among novel loci, missense variants linked (r2 > 0.7) with lead variants included the genes CCDC18, MERTK, APOL3, PPARG, MTTP, MLXIPL, ABCB4, AKNA, GPAM, SH2B3, P2RX7, NYNRIN, ANPEP,IFI30 and MPV17L2. Among the trans-ancestry coding missense variants, ten (CCDC18, MLXIPL, ABCB4, AKNA, DNAJC22, SERPINA1, ANPEP, APOH, IFI30, MPV17L2, and PNPLA3) were predicted based on two methods (SIFT, PolyPhen-2) to have potentially deleterious and/or damaging effects in protein function34,35. An AA-specific locus on chromosome (rs115038698, chr7:87024718) was strongly linked to a nearby missense variant Ala934Thr in ABCB4 (rs61730509, AFR r2=0.92) with predicted deleterious effect, where the T-allele confers an increased risk of NAFLD (β=0.617, P=1.8×10−20). In summary, 24 of our 77 trans-ancestry loci prioritized a candidate gene based on a missense variant in tight linkage with the lead SNP (Supplemental Table 25).
Additional approaches to nominate putative causal genes
We performed colocalization analyses with gene expression and splicing across 48 tissues measured by the GTEx project, and overlapped our lead SNPs with histone quantitative trait locus (QTL) data from livers to identify NAFLD-associated variants that are also associated with change in gene expression (eQTLs), splice isoforms (sQTLs), or histone modifications (hQTLs, Methods, Supplemental Table 29). Across all tissues, a total of 123 genes were prioritized with 20 genes in liver tissue (Methods). In liver tissue alone, a total of eight variant-gene pairs were identified where the allele associated with protection against NAFLD was also associated with reduced gene expression (i.e., the direction of effect was concordant between the GTEx eQTL and GWAS sentinel variant): AC091114.1, PANX1, FADS2, SHROOM3, U2SURP, NYNRIN, CD276 and EFHD1. Furthermore, sQTL analysis in GTEx v8 identified two genes in liver, HSD17B13 and ANPEP, and 12 genes (MARC1, HSD17B13, ABO, FADS1-FADS2, TMEM258, MLXIP, ANPEP, KAT7, STRADA, DDX42, TRC4AP, and APOL3) that were affected in at least two tissues (Supplemental Table 30). Finally, two of our lead SNPs were in high LD (r2 > 0.8) with variants that regulated H3K27ac levels in liver tissue (hQTLs), namely EFHD1 (hQTL SNPs rs2140773, rs7604422 in EFHD1) and FADS2 loci (hQTL rs174566 in FADS2)36.
We next mapped our NAFLD loci to regions of open chromatin using ATAC-seq in three biologically-relevant liver-derived tissues (human liver, liver cancer cell line [HepG2], and hepatocyte-like cells [HLC] derived from pluripotent stem cells)37. Additionally, we used promoter-focused Capture-C data to identify those credible sets that physically interact with genes in two relevant cell types (HepG2 and liver) (Supplemental Table 31). These datasets are useful entry points for deciphering regulatory mechanisms involved in the pathophysiology of NAFLD. Most notably, the genes DHODH, H2AFZ, PAQR9, FTO, MIR644A, BCL7B and KRT82 showed interactions with NAFLD loci that were also in open chromatin in both HLC and HepG2 cells.
Based on DEPICT gene prediction, coding variant linkage analysis, and QTL colocalization (Supplemental Tables 24-31), 215 potentially relevant genes for NAFLD were identified for the 77 loci. A protein-protein interaction (PPI) analysis revealed that among the 192 available proteins, 86 nodes were observed, with a PPI enrichment (P < 9.0×10−8) indicating that the network has substantially more interactions than expected by chance (Supplemental Table 32 and Supplemental Figure 6).
For each gene identified from all of the above described analyses, we counted the number of times that the gene was identified for each of the analyses (DEPICT gene prediction, coding variant linkage, QTL colocalization, promotor Capture-C and/or ATAC-Seq peak overlap, and PPI network analysis) and divided this by the number of analyses (e.g., 8). We labeled this measure as the gene nomination score, which reflects the cumulative evidence supporting the respective gene as causal for the observed association. Based on our gene nomination scheme, we found evidence for a single gene nomination at 52 genomic loci, two genes at 14 loci, and three genes at six loci. Six loci had more than three genes nominated (one of which was HLA), and only two loci lacked any data to support a nomination (Supplemental Table 33). We further prioritized those loci which were prioritized by at least 3 sources of evidence (or 4 sources of evidence for coding variants). This resulted in a total of 27 loci supported by multiple lines of evidence (Table 1), which included 6 loci with co-localizing eQTLs in liver or adipose tissues and connection to the predicted gene via Promoter CaptureC data (i.e., EPHA2, IL1RN, SHROOM3, HKDC1, PANX1, DHODH;HP).
Interestingly, 14 of the nominated genes are transcription factors (TF) (Supplemental Table 34). Of particular interest, two of these TFs have several downstream target genes identified using the DoRothEA data in OmniPath (Methods). Notably, the CEBPA TF targets the downstream genes PPARG, TRIB1, GPAM, FTO, IRS1, CRIM1, HP, TBC1D8, and CPS1, but also NCEH1, a gene in the vicinity of one of our associations that lacked a nominated candidate gene. Similarly, HNF1A, the lead gene in EA scan (and corresponding to the trans-ancestry P2RX7 locus) targets SLC2A2, MTTP, and APOH.
Polygenic Risk Score analyses
We calculated a candidate SNP polygenic risk score (PRS) based on Stage 1 350K dataset (primary set) to perform a Phewas in an independent sample in MVP (Stage 2 replication set). We observed that an increased NAFLD PRS was associated with abnormal results of function study of liver (Bonferroni P < 3.1 x 10−5), and showed suggestive significance with bacterial pneumonia, otalgia, gout and other crystal arthropathies and non-infectious gastroenteritis (P < 0.001, Supplemental Table 35). Furthermore, a NAFLD PRS based on the Stage 1 set was positively associated with NAFLD prediction in the Stage 2 replication set (P=3.8 x 10−5, Supplemental Table 36).
Investigation of pleiotropy of lead NAFLD SNPs
We next sought to identify additional traits that were also associated with our 77 trans-ancestry lead SNPs. First, we performed a LabWAS of distinct clinical laboratory test results38 in MVP (Methods), yielding 304 significant SNP-trait associations (Supplemental Table 37, Supplemental Figure 7). Second, we performed a PheWAS Analysis in UK Biobank data using SAIGE (Methods), which identified various SNP-trait associations that mapped to loci previously associated with liver traits, cardiometabolic traits, as well as additional enriched association for gallstones, gout, arthritis, and hernias (Supplemental Tables 38 and 39). In particular, we examined all associations for PheCode 571.5, “Other chronic nonalcoholic liver disease” which comprised 1,664 cases and 400,055 controls. Of the n=73 variants found, we noted that 14/73 were both nominally associated and directionally consistent with our scan (signed binomial test P=3.4×10−9), providing additional validation for our scan (Supplementary Table 40). Third, we performed a SNP lookup using the curated data in the IEU OpenGWAS project (Supplemental Tables 41 and 42), which identified 2,891 genome-wide significant SNP-trait associations for trans-ancestry SNPs, and additional 283 SNP-trait associations for the ancestry-specific lead SNPs. Finally, we performed cross-trait colocalization analyses using COLOC of EA, AA, and HISP lead loci with 36 other GWAS statistics of cardiometabolic and blood cell related traits (Methods). This resulted in significant regional colocalization for 64 SNP-trait pairs in EA, 32 SNP-trait pairs in AA, and 12 SNP trait pairs in HISP (Supplemental Table 43-45).
Based on the four analyses described above, we categorized relevant phenotypes observed as liver traits, metabolic traits, or inflammatory traits based on all significant SNP-trait associations and their nominated candidate genes (Supplemental Tables 37-44, Figure 4). Across the trans-ancestry lead variants (n=77), ancestry-specific (n=3), and secondary proximal associations (HNF1A, n=1), 22 SNPs showed association with only liver traits (such as ALT, ALP, AST, and GGT) (Figure 4). By contrast, 23 loci showed associations with both liver and cardiometabolic traits (such as HDL, LDL, and total cholesterol, triglycerides, BMI, glucose, and HbA1c) whereas 3 loci (IL1RN, TMEM147;ATP4A and RORA) showed association with both liver traits and inflammatory traits (e.g., C-reactive protein, white blood cell count, lymphocyte count). Finally, 25 loci showed association with all three traits: liver, cardiometabolic, and inflammation. Notably, among 12 loci that showed significant association with hepatic fat (color-coded in red in Figure 4), 11 were associated with both liver and metabolic traits, including five that were also associated with inflammatory traits. Collectively, our findings identify novel NAFLD-associated genetic loci with pleotropic effects that may impact hepatic, metabolic and inflammatory traits.
Overlapping liver (light blue), metabolic (pink) and/or inflammatory (green) traits are shown in association with 77 trans-ancestry and additional ancestry-specific lead SNPs. The trait categorizations reflect significant SNP-trait associations identified by: 1) LabWAS of clinical laboratory results in MVP; 2) PheWAS with UKBB data using SAIGE; 3) SNP lookup using the curated data in the IEU OpenGWAS projects; and 4) cross-trait colocalization analyses using COLOC of EA, AA and HISP lead loci with 36 other GWAS statistics of cardiometabolic and blood cell related traits. Red/bold font refers to the loci also associated with quantitative hepatic fat on imaging analyses.
Discussion
In this study, the largest and most diverse GWAS of NAFLD to date, we report a total of 77 trans-ancestry (of which 70 are novel) and 3 additional ancestry-specific loci that show significant genome-wide association with NAFLD. While our NAFLD definition is a proxy for chronic hepatocellular injury in the absence of other known causes of liver disease, we further used CT/MRI imaging data to compare to what extent these SNPs also associated with hepatic fat accumulation. Overall, 24 (∼30%) of these loci were nominally associated with hepatic fat based on CT or MRI, and the majority of these overlapping SNPs were associated with metabolic and/or inflammatory traits. Thus, SNPs that are associated with liver enzymes, metabolic risk factors, and inflammatory biomarkers may be the most likely to be associated with liver steatosis and should be prioritized for further follow-up. Furthermore, detailed genetic correlation analyses showed significant enrichment of these SNPs for cardiometabolic traits, metabolic pathways, and genomic annotations relevant for NAFLD. We found that most of our index NAFLD-associated SNPs were associated with metabolic and/or inflammatory traits - the most common being lipid-related, followed by glycemic traits, hypertension, and cardiovascular disease, as well as cholelithiasis (gallstones), cholecystitis, osteoarthritis, hypothyroidism, and thrombophlebitis. Collectively, our findings offer a comprehensive and refined view of the genetic contribution to NAFLD with potential clinical, pathogenic, and therapeutic relevance. Integration of these with extant phenotypic association data sets allowed us to further characterize the functional mechanisms through which our identified loci may mediate NAFLD risk.
Previous studies for liver enzyme levels, particularly serum ALT activity, have been performed 10,11,16. While there is overlap in the discoveries made by studies of natural variation in circulating levels of this biomarker, our cohort and approach to phenotyping make our results and interpretation unique. First, the diversity of our cohort provides both additional power and potential for discovery, as the bulk of studies to date have been performed in predominantly European-ancestry cohorts. Second, our approach ascertains individuals with chronic elevation of this enzyme, consistent with genuine chronic liver disease. At the same time, we excluded individuals with known causes of liver disease outside of NAFLD via ICD code definition, which served to further enrich for metabolic disorders in our cohort. We further excluded control individuals who maybe have intermittent ALT elevation, focusing on a healthier, ‘super-control’ subset of the population. The result is that our approach should have higher specificity to ascertained risk alleles that predispose to metabolic-induced fatty liver disease. In contrast, a standard-ALT scan would be powered to discover the full spectrum of causes of liver disease (and perhaps many more loci), many of which will not be specific to NAFLD and may be due to other causes. As we have shown in validation studies using quantitative measures of hepatic fat as well as ICD-code definitions of NAFLD, our results are highly directionally concordant, demonstrating the relevance of our proxy phenotype to liver disease and physiology. Genetic correlation analysis demonstrated strong correlation with cardiometabolic traits and disease, again consistent with the relevance of our trait relative to simple enzyme measures.
There are several aspects of our study that are worth highlighting. We demonstrate the strength of trans-ancestry GWAS for the discovery and interrogation of NAFLD susceptibility loci, discoveries made possible by the diversity and sample size of the Million Veteran Program cohort, of which 25% of participants are of non-European ancestry. Utilizing this data allows us to narrow down putatively causal variants through trans-ancestry fine-mapping and construction of credible sets likely to harbor the likely culprit variant(s). Construction of credible sets using trans-ancestry data has been shown to facilitate fine-mapping by producing smaller credible sets compared to sets based on single ancestries39, an effect we also observed at our loci. Moreover, we identified eight NAFLD-associated loci in AAs. In particular, the lead SNP at the ABCB4 locus (rs115038698) was in high LD with the missense variant rs61730509 (Ala934Thr, AFR r2=0.92) and segregated a very potent effect (OR=1.87, CI=1.64-2.14, P=1.8×10−20). This variant is of low frequency in AA (MAF=1.2%) but virtually absent in EA and ASN. ABCB4, also known as multidrug resistance protein 3 (MDR3), is a compelling candidate gene, as it has been previously implicated in cholestasis, gallbladder disease, and adult biliary fibrosis/cirrhosis40-42. Finally, for a number of variant gene-pairs, the observed effect on NAFLD risk and the impact of gene expression in the liver was consistent with our understanding of the expected effect given what is known about gene function, suggesting possible relevance as therapeutic targets. Among those, genetic deletion of Pannexin 1 (encoded by PANX1) was reported to have a protective effect in mouse model of acute and chronic liver disease43,44, and is consistent with the data we report here.
Twelve of our loci were associated with quantitative measures of hepatic fat after multiple-test correction. These included loci previously associated with NAFLD or all-cause cirrhosis (e.g., PNPLA3, TM6SF2, TNKS [PPP1R3B], KIAA0196 [TRIB1], and MARC1), but also included novel loci reported here (e.g., GPAM, APOE;APOC1, MTTP, APOH, IFI30;MPV17L2, SCN2A;COBLL1, and PPARG). In all cases except TNKS [PPP1R3B], the directional effect on hepatic fat was consistent with cALT levels. A discordance between measures of hepatic fat based on radiological and histological evaluation has been noted12 and may be explained by the role of the PPP1R3B-encoded protein in promoting the accumulation of hepatic glycogen45 which may influence the contrast in hepatic images46,47.
Through functional genomic and bioinformatics prioritization analyses beyond those based on coding variants or eQTLs, we were able to nominate loci that have at least one candidate gene nominations at 75 out of our 77 (97%) identified loci. We found that these genes were often highly expressed in liver and have prior biological connections to liver physiology and disease, making this list compelling for further interrogation. As an example, GPAM, tagged by the missense variant rs2792751 (Ile43Val, EA r2 = 0.99), encodes the mitochondrial glycerol-3-phosphate-acyltransferase 1, a protein used in the mitochondria to convert saturated fatty acids into glycerolipids. GPAM is highly expressed in liver48,49 and associated with metabolic disease50, consistent with our pleiotropy analyses. Mouse knockouts of GPAM had reduced weight, lower hepatic triacylglycerol content, and decreased plasma triacylglycerol51. Another example is MTTP which is tagged by the missense variant rs3816873 (Ile128Thr, EA r2=1.0) and encodes the microsomal triglyceride transfer protein, which loads lipids onto assembling VLDL particles and facilitate their secretion by hepatocytes. Liver-specific MTTP knockout mice have reduced VLDL secretion and increased hepatic steatosis52. Lomitapide, a small molecule inhibitor of MTTP, is approved as a treatment for lowering LDL cholesterol in homozygous familial hypercholesterolemia, but increases liver lipid by inhibiting VLDL secretion53. TRIM5 (Val112Phe) is a member of the tripartite motif (TRIM) family with E3 ubiquitin ligase activity with a key role in innate immune signaling and antiviral host defense54, and TRIM5 SNPs have been associated with increased risk of liver fibrosis in HIV/HCV co-infected patients55. APOH (Cys325Gly) encodes the apolipoprotein H which is exclusively expressed in liver tissue48 and which is associated with ALT, AST, triglycerides, LDL cholesterol and platelets in the MVP labWAS. Two coding variants (strongly linked) in MerTK (Arg466Lys and Ile518Val, r2=0.98) were associated with NAFLD; MerTK signaling in hepatic macrophages was recently shown to mediate hepatic stellate cell activation and promote hepatic fibrosis progression56, and variants in MERTK were associated with liver fibrosis progression in HCV-infected patients57, raising the possibility for MerTK as a novel therapeutic target against fibrosis58. We emphasize that functional studies of our nominated causal genes are needed to demosntrate casual relevance, their impact on hepatosteatosis, and ultimately to determine their underlying mechanisms.
Given the complex etiology and progression of NAFLD, we anticipated that our study would identify novel loci with putatively causal genes that span multiple molecular pathways. Indeed, our novel loci include genes that play roles in obesity (e.g., FTO, PPARG), insulin resistance (e.g., COBLL1, MIR5702 [IRS1]), and diabetes (e.g., HNF1A). Relevant for hepatic inflammation in the two-hit hypothesis of NAFLD59, our novel loci also implicate immune-mediated or inflammatory contributions to NAFLD progression, including HLA, RORA60,61, IFI3062,63, CD27664, ILRN62,65, ITCH66,67 and P2RX768-70. Among these, RORA encodes the retinoic acid receptor related orphan receptor A which may be involved in NASH pathogenesis through macrophage polarization and miRNA122, which comprises 70% of the total miRNA in liver 60,61. It is also known that loss of TRIB1 substantially decreases miR-122 levels via its impact on HNF4 and HNF1A71. IFI30 encodes gamma-interferon-inducible lysosomal thiol reductase (GILT) which is involved in antigen processing and presentation and the production of reactive oxygen species during cellular stress and autophagy. Finally, P2RX7 encodes the purinergic receptor P2×7 which is involved in inflammasome activation and IL-1β processing in liver inflammation and fibrosis68-70. Encouragingly, these and additional pathways have emerged despite the proxy nature of our phenotype, and almost certainly underestimate the true number of loci contributing to NAFLD.
In conclusion, we define 77 trans-ancestry loci (70 novel) with 3 additional ancestry-specific loci associated with NAFLD by using chronic ALT elevation in a large, ancestrally diverse cohort enriched for metabolic disorders without other known causes of liver disease. The abundance of NAFLD loci identified by our analyses constitutes a much-needed large-scale, multi-ancestry genetic resource that can be used to build prediction models, identify causal mechanisms, and understand biological pathways contributing to NAFLD initiation and disease progression.
Methods
We performed a large-scale trans-ancestry NAFLD GWAS in the Million Veteran Program. We subsequently conducted analyses to facilitate the prioritization of these individual findings, including transcriptome-wide predicted gene expression for NAFLD, secondary signal analysis, coding variant mapping, phenome-wide association analyses in various public data sources, and various forms of cardiometabolic cross-trait colocalization analyses to fine-map the genomic loci to putatively causal genes.
Discovery cohort
The Million Veteran Program (MVP) is a large cohort of fully consented veterans of the United States military forces recruited from 63 participating Department of Veterans Affairs (VA) medical facilities28. Recruitment for this ongoing sample started in 2011, and all veterans are eligible to participate. This study analyzed clinical data through July 2017 for participants who were enrolled since January 2011. All MVP study participants provided blood samples for DNA extraction and genotyping, completed surveys about their health, lifestyle, and military experiences. Consent to participate and permission to re-contact was provided after veterans received information materials by mail and met with research staff to address their questions. Study participation also includes access to the participant’s electronic health records for research purposes. Each veteran’s electronic health care record is integrated into the MVP biorepository, including inpatient International Classification of Diseases (ICD-9-CM and ICD-10-CM) diagnosis codes, Current Procedural Terminology (CPT) procedure codes, clinical laboratory measurements, and reports of diagnostic imaging modalities. Researchers are provided with de-identified data, and have neither the ability nor authorization to link these details with a participant’s identity. Blood samples are collected by phlebotomists and banked at the VA Central Biorepository in Boston, where DNA is extracted and shipped to two external centers for genotyping. The MVP received ethical and study protocol approval from the VA Central Institutional Review Board (IRB) in accordance with the principles outlined in the Declaration of Helsinki.
Genotyping
DNA extracted from buffy coat was genotyped using a custom Affymetrix Axiom biobank array. The MVP 1.0 genotyping array contains a total of 723,305 SNPs, enriched for 1) low frequency variants in AA and HISP populations, and 2) variants associated with diseases common to the VA population 28.
Genotype quality-control
The MVP genomics working group applied standard quality control and genotype calling algorithms to the data in three batches using the Affymetrix Power Tools Suite (v1.18). Excluded were duplicate samples, samples with more heterozygosity than expected, and samples with an over 2.5% missing genotype calls. We excluded related individuals (halfway between second- and third-degree relatives or closer) with KING software72. Before imputation, variants that were poorly called or that deviated from their expected allele frequency based on reference data from the 1000 Genomes Project were excluded73. After prephasing using EAGLE v2, genotypes were imputed via Minimac4 software74 from the 1000 Genomes Project phase 3, version 5 reference panel. The top 30 principal components (PCs) were computed using FlashPCA in all MVP participants and an additional 2,504 individuals from 1000 Genomes. These PCs were used to unify of self-reported race/ancestry and genetically inferred ancestry to compose ancestral groups29.
Phenotype classification
MVP NAFLD phenotype definitions were developed by combining a previously published VA CDW ALT-based approach with non-invasive clinical parameters available to practicing clinicians at the point of care. The primary NAFLD phenotype (labeled “ALT-threshold”) was defined by: (i) elevated ALT >40 U/L for men and >30 U/L for women during at least two time points at least 6 months apart within a two-year window period at any point prior to enrollment and (ii) exclusion of other causes of liver disease (e.g. presence of chronic viral hepatitis B or C [defined as positive hepatitis C RNA > 0 international units/mL or positive hepatitis B surface antigen], chronic liver diseases or systemic conditions [e.g. hemochromatosis, primary biliary cholangitis, primary sclerosing cholangitis, autoimmune hepatitis, alpha-1-antitrypsin deficiency, sarcoidosis, metastatic liver cancer, secondary biliary cirrhosis, Wilson’s disease], and/or alcohol use disorder [e.g. alcohol use disorder, alcoholic liver disease, alcoholic hepatitis and/or ascites, alcoholic fibrosis and sclerosis of liver, alcoholic cirrhosis of liver and/or ascites, alcoholic hepatic failure and/or coma, and unspecified alcoholic liver disease). The control group was defined by having a: normal ALT (≤30 U/L for men, ≤20 U/L for women) and no apparent causes of liver disease or alcohol use disorder or related conditions21. Habitual alcohol consumption was assessed with the age-adjusted Alcohol Use Disorders Identification Test (AUDIT-C) score, a validated questionnaire annually administered by VA primary care practitioners and used previously in MVP75,76.
Single-variant autosomal analyses
We tested imputed SNPs that passed quality control (i.e. HWE > 1×10−10, INFO > 0.3, call rate > 0.975) for association with NAFLD through logistic regression assuming an additive model of variants with MAF > 0.1% in European American (EA), and MAF > 1% in African Americans (AA), Hispanics (HISP), and Asians (ASN) using PLINK2a software77. Covariates included age, gender, age-adjusted AUDIT-C score, and 10 principal components of genetic ancestry. We aggregated association summary statistics from the ancestry-specific analyses and performed a trans-ancestry meta-analysis. The association summary statistics for each analysis were meta-analyzed in a fixed-effects model using METAL with inverse-variance weighting of log odds ratios78. Variants were clumped using a range of 500kb and/or CEU r2 LD > 0.05, and were considered genome-wide significant if they passed the conventional p-value threshold of 5.0×10−8.
Secondary signal analysis
GCTA79 was used to conduct conditional analyses to detect ancestry-specific distinct association signals at each of the lead SNPs utilizing the GWAS summary statistics in EA, AA, and HISP; these ancestry-stratified MVP cohorts were used to model LD patterns between variants. The reference panel of genotypes consisted of the variants with allele frequencies > 0.1% in EA, >1% in AA, and >1% in HISP that passed quality control criteria in the MVP-specific NAFLD GWAS (INFO > 0.3, HWE P > 1.0×10−10, call rate > 0.975). For each lead SNP, conditionally independent variants that reached locus-wide significance (P < 1.0×10−5) were considered secondary signals of distinct association. If the minimum distance between any distinct signals from two separate loci was less than 500kb, we performed an additional conditional analysis that included both regions and reassessed the independence of each signal.
Credible Sets
We calculated Wakefield’s approximate Bayes’ factors 30 based on the marginal summary statistics of the trans-ancestry meta-analysis and ancestry specific summary statistics using the CRAN R package corrcoverage80. For each locus, the posterior probabilities of each variant being causal were calculated and a 95% credible set was generated which contains the minimum set of variants that jointly have at least 95% posterior probability (PP) of including the causal variant.
Concordance of NAFLD with qHF
For 77 lead trans-ancestry SNPs a concordance analysis was performed to evaluate the extent to which genetic predictors of hepatocellular injury (cALT) correspond with quantitative hepatic fat derived from computed tomography (CT) / magnetic resonance imaging (MRI)-measured hepatic fat in the Penn Medicine Biobank (PMBB), UK Biobank, Multi-Ethnic Study of Atherosclerosis (MESA), Framingham Heart Study (FHS), and University of Maryland Older Order Amish study. Attenuation was measured in Hounsfield units. The difference between the spleen and liver attenuation was measured for PMBB; a ratio between liver attenuation/spleen attenuation was used for MESA and Amish; and liver attenuation/phantom attenuation ratio in FHS as previously described by Speliotes et al12. Abdominal MRI data from UK Biobank data were used to quantify liver fat using a two-stage machine learning approach with deep convolutional neural networks81. CT-measured hepatic fat was estimated using a multi-stage series of neural networks for presence of scan contrast and liver segmentation using convolutional neural networks. The PMBB included CT data on 2,979 EA and 1,250 AA participants82, the FHS included a total of 3,011 EA participants, the Amish study 754 EA participants, and MESA contributed 1,525 EA, 1,048 AA, 923 HISP, and 360 ASN participants for concordance analysis. The UK Biobank included MRI image data from 36,703 EA participants. All cohorts underwent individual-level linear regression analysis on hepatic fat, adjusted for the covariates of age, gender, first 10 principal components of genetic ancestry, and alcohol intake if available. If the lead SNP was not available in any of the studies, a proxy SNP in high LD with the lead variant was used (r2 > 0.7) or if no such variant was identified, the SNP was set to missing for that respective study. The study-specific ancestry-stratified summary statistics were first standardized to generate standard scores or normal deviates (z-scores), and then meta-analyzed using METAL in a fixed-effects model with inverse-variance weighting of regression coefficients78. In a first round of meta-analysis, ancestry-specific summary statistics were generated, which then served as input for a subsequent round of meta-analysis that represents the trans-ancestry effects of our lead SNPs on quantitative hepatic fat.
Heritability estimates and genetic correlations analysis
LD-score regression was used to estimate the heritability coefficient, and subsequently population and sample prevalence estimates were applied to estimate heritability on the liability scale83. A genome-wide genetic correlation analysis was performed to investigate possible co-regulation or a shared genetic basis between T2D and other complex traits and diseases. Pairwise genetic correlation coefficients were estimated between the meta-analyzed NAFLD GWAS summary output in EA and each of 774 precomputed and publicly available GWAS summary statistics for complex traits and diseases by using LD score regression through LD Hub v1.9.3 (http://ldsc.broadinstitute.org). Statistical significance was set to a Bonferroni-corrected level of P < 6.5 x 10−5.
Tissue- and epigenetic-specific enrichment of NAFLD heritability
We analyzed cell type-specific annotations to identify enrichments of NAFLD heritability. First, a baseline gene model was generated consisting of 53 functional categories, including UCSC gene models, ENCODE functional annotations84, Roadmap epigenomic annotations85, and FANTOM5 enhancers86. Gene expression and chromatin data were also analyzed to identify disease-relevant tissues, cell types, and tissue-specific epigenetic annotations. We used LDSC31-33 to test for enriched heritability in regions surrounding genes with the highest tissue-specific expression. Sources of data that were analyzed included 53 human tissue or cell type RNA-seq data from GTEx27; human, mouse, or rat tissue or cell type array data from the Franke lab87; 3 sets of mouse brain cell type array data from Cahoy et al88; 292 mouse immune cell type array data from ImmGen89; and 396 human epigenetic annotations from the Roadmap Epigenomics Consortium 85.
Pathway Annotation enrichment
Enrichment analyses in DEPICT90 were conducted using genome-wide significant (P < 5×10−8) NAFLD GWAS lead SNPs. DEPICT is based on predefined phenotypic gene sets from multiple databases and Affymetrix HGU133a2.0 expression microarray data from >37k subjects to build highly-expressed gene sets for Medical Subject Heading (MeSH) tissue and cell type annotations. Output includes a P-value for enrichment and a yes/no indicator of whether the FDR q-value is significant (P < 0.05). Tissue and gene-set enrichment features are considered. We tested for epigenomic enrichment of genetic variants using GREGOR software91. We selected EA-specific NAFLD lead variants with a p-value less than 5×10−8. We tested for enrichment of the resulting GWAS lead variants or their LD proxies (r2 threshold of 0.8 within 1 Mb of the GWAS lead, 1000 Genomes Phase I) in genomic features including ENCODE, Epigenome Roadmap, and manually curated data (Supplemental Table 24). Enrichment was considered significant if the enrichment p-value was less than the Bonferroni-corrected threshold of P=1.8×10−5 (0.05/2,725 tested features).
Coding variant mapping
All imputed variants in MVP were evaluated with Ensemble variant effect predictor92, and all predicted LoF and missense variants were extracted. The LD was calculated with established variants for trans-ancestry, EA, AA, and HISP lead SNPs based on 1000 Genomes reference panel73. For SNPs with low allele frequencies, the MVP dataset was used for LD calculation for the respective underlying population. For the trans-ancestry coding variant, the EA panel was used for LD calculation. Coding variants that were in strong LD (r2 > 0.7) with lead SNPs and had a strong statistical association (P-value < 1×10−5) were considered the putative causal drivers of the observed association at the respective locus.
Colocalization with gene expression
GWAS summary statistics were lifted over from GRCh37 to GRCh38 using LiftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver). Colocalization analysis was run separately for each of the 49 tissues in GTEx v827. For each tissue, we obtained an LD block for the genome with a sentinel SNP at P < 5×10−8, and then restricted analysis to the LD blocks. For each LD block with a sentinel SNP, all genes within 1Mb of the sentinel SNP (cis-Genes) were identified, and then restricted to those that were identified as eGenes in GTEx v8 at an FDR threshold of 0.05 (cis-eGenes). For each cis-eGene, we performed colocalization using all variants within 1Mb of the gene using the default prior probabilities in the ‘coloc’ function for the coloc package in R. We first assessed each coloc result for whether there was sufficient power to test for colocalization (PP3+PP4>0.8), and for the colocalization pairs that pass the power threshold, we defined the significant colocalization threshold as PP4/(PP3+PP4)>0.9.
Overlap with open chromatin
At each of the 77 NAFLD-associated loci from the trans-ancestry meta-analysis, we looked for overlaps between any variant in the credible set, and regions of open chromatin previously identified using ATAC-Seq experiments in two cell types—3 biological replicates of HepG293 and 3 biological replicates of hepatocyte-like cells (HLC)94 produced by differentiating three biological replicates of iPSCs, which in turn were generated from peripheral blood mononuclear cells using a previously published protocol36.
Overlap with Promoter Capture-C data
We used two promoter Capture-C datasets from two cell/tissue types to capture physical interactions between gene promoters and their regulatory elements and genes; three biological replicates of HepG2 liver carcinoma cells, and hepatocyte-like cells (HLC)93. The detailed protocol to prepare HepG2 or HLC cells for the promoter Capture-C experiment is previously described36. Briefly, for each dataset, 10 million cells were used for promoter Capture-C library generation. Custom capture baits were designed using an Agilent SureSelect library design targeting both ends of DpnII restriction fragments encompassing promoters (including alternative promoters) of all human coding genes, noncoding RNA, antisense RNA, snRNA, miRNA, snoRNA, and lincRNA transcripts, totaling 36,691 RNA baited fragments. Each library was then sequenced on an Illumina NovoSeq (HLC), or Illumina HiSeq 4000 (HLC), generating 1.6 billion read pairs per sample (50 base pair read length.) HiCUP95 was used to process the raw FastQ files into loop calls; we then used CHiCAGO96 to define significant looping interactions; a default score of 5 was defined as significant. We identified those NAFLD loci at which at least one variant in the credible set interacted with an annotated bait in the Capture-C data.
Protein-Protein Interaction Network Analysis
We employed the search tool for retrieval of interacting genes (STRING) v1197 (https://string-db.org) to seek potential interactions between nominated genes. STRING integrates both known and predicted PPIs and can be applied to predict functional interactions of proteins. In our study, the sources for interaction were restricted to the ‘Homo Sapiens’ species and limited to experimentally validated and curated databases. An interaction score > 0.4 were applied to construct the PPI networks, in which the nodes correspond to the proteins and the edges represent the interactions (Figure 4, Supplemental Table 32).
Gene Nomination
Based on DEPICT gene prediction, coding variant linkage analysis, QTL analysis, and annotation enrichment, and PPI networks (Supplemental Tables 24-33), a total of 215 potentially relevant genes for NAFLD were mapped to trans-ancestry 77 loci. For each locus with multiple mapped genes, we counted how many times each gene was identified through each of the analysis, and divided this by the total number of experiments (i.e., 8) to calculate an evidence burden that ranges from 0 to 100%. For each genomic locus, the gene that was most frequently identified as potentially relevant was selected as the putative causal gene. In the case of a tie break, and if the respective genes have identical nomination profiles, the gene with more eQTLs was selected as the putative causal gene. Similarly, gene nomination was preferred for loci that strongly tagged (r2 > 0.8) a coding variant. Loci that scored with 3 pieces of evidence or greater are listed for coding variant (Table 1A) and non-coding variants (Table 1B), respectively.
MVP LabWAS
A total of 21 continuous traits in the discovery MVP dataset, e.g. AST, ALP, fasting TG, HDL, LDL, TC, random glucose, HbA1c, albumin, bilirubin, platelet count, BMI, blood urea nitrogen (BUN), creatinine, eGFR, SBP, DBP, ESR, INR, and C-reactive protein were tested in 186,681 EA’s with association of 77 SNPs using linear regression of log-linear values. Covariates included age, gender and first 10 principal components of EA ancestry.
PheWAS with UK Biobank data
For the 77 lead trans-ancestry SNPs and EA and AA specific SNPs, we performed a PheWAS in a genome-wide association study of EHR-derived ICD billing codes from the White British participants of the UK Biobank using PheWeb98. In short, phenotypes were classified into 1,403 PheWAS codes excluding SNP-PheWAS code association pairs with case counts less than fifty99. All individuals were imputed using the Haplotype Reference Consortium panel 100, resulting in the availability of 28 million genetic variants for a total of 408,961 subjects. Analyses on binary outcomes were conducted using a model named SAIGE, adjusted for genetic relatedness, gender, year of birth and the first 4 principal components of white British genetic ancestry101. SAIGE stands for Scalable and Accurate Implementation of GEneralized mixed model and represents a generalized mixed-model association test that accounts for case-control imbalance and sample relatedness101.
IEU OpenGWAS project SNP lookup
An additional phenome-wide lookup was performed for 77 lead trans-ancestry SNPs and EA and AA specific SNPs in Bristol University’s MRC Integrative Epidemiology Unit (IEU) GWAS database102. This database consists of 126,114,500,026 genetic associations from 34,494 GWAS summary datasets, including UK Biobank (http://www.nealelab.is/uk-biobank), FinnGen (https://github.com/FINNGEN/pheweb), Biobank Japan (http://jenger.riken.jp/result), NHGRI-EBI GWAS catalog (https://www.ebi.ac.uk/gwas), blood metabolites GWAS103, circulating metabolites GWAS104, the MR-Base manually curated database105, and protein level GWAS106.
Regional cardiometabolic cross-trait colocalization
Bayesian colocalization tests between NAFLD-associated signals and the following trait- and disease-associated signals were performed using the COLOC R package107. To enable cross-trait associations, we compiled summary statistics of 36 cardiometabolic and blood cell-related quantitative traits and disease from GWAS studies conducted in EA ancestry individuals, and for MVP-based reports also on AA and HISP. To summarize, for total, HDL, and LDL cholesterol, triglycerides, alcohol use disorder, alcohol intake, systolic blood pressure, diastolic blood pressure, type 2 diabetes, BMI, CAD, we used the summary statistics available from various MVP-based studies26,75,108. Of these, the summary statistics for CAD and BMI GWAS have not been published or deposited as of yet. Data on WHR were derived from GIANT Consortium109, whereas summary statistics on CKD, gout, blood urea nitrogen, urate, urinary albumin-to-creatinine ratio, microalbuminuria, and eGFR were derived from CKD Genetics Consortium110-112. Finally, summary statistics of blood cell traits (e.g. platelet count, albumin, white blood cells, basophils, eosinophils, neutrophils, hemoglobin, hematocrit, immature reticulocyte fraction, lymphocytes, monocytes, reticulocytes, mean corpuscular hemoglobin, mean corpuscular volume, mean platelet volume, platelet distribution width, and red cell distribution width) were derived from a large-scale GWAS report performed in UK Biobank and INTERVAL studies113. A colocalization test was performed for all 77 NAFLD loci spanning 500kb region around the lead SNP for all 36 compiled traits. COLOC requires for each SNP data on allele frequency, sample size, beta-coefficients and variance or p values. For each association pair COLOC was run with default parameters and priors. COLOC computed posterior probabilities for the following five hypotheses: PP0, no association with trait 1 (NAFLD GWAS signal) or trait 2 (e.g., co-associated metabolic signal); PP1, association with trait 1 only (i.e., no association with trait 2); PP2, association with trait 2 only (i.e., no association with trait 1); PP3, association with trait 1 and trait 2 by two independent signals; and PP4, association with trait 1 and trait 2 by shared variants. Evidence of colocalization114 was defined by PP3 + PP4 ≥ 0.99 and PP4/PP3 ≥ 5.
NAFLD Polygenic risk score and NAFLD risk
We constructed polygenic risk score (PRS) for NAFLD in the Stage 2 replication data set containing of 73,580 MVP participants of EA ancestry by calculating a linear combination of weights derived from the discovery MVP dataset of lead 77 trans-ancestry variants. The PRS was divided into quintiles and the risk of NAFLD was assessed using a logistic regression model using the lowest decile as a reference (e.g. the 20% of participants with lowest of NAFLD PRS), together with the potential confounding factors of age, gender, age-adjusted AUDIT-C, and the first 10 principal components of EA ancestry.
NAFLD PRS Phewas
For the NAFLD PRS that was generated using the Stage 1 350K dataset, we performed a PheWAS study in the Stage 2 108K replication dataset to fully leverage full catalog of available ICD-9/ICD-10 diagnosis codes. Of genotyped veterans, participants were included in the PheWAS analysis if their respective electronic health record reflected two or more separate encounters in the VA Healthcare System in MVP up to July 2017. Using this method, a total of 73,580 veterans of EA ancestry were available for PheWAS analysis. ICD-9/ICD-10 diagnosis codes were collapsed to clinical disease groups and corresponding controls using predefined groupings 99. Phenotypes were required to have a case count over 25 in order to be included in the PheWAS analysis, and a multiple testing thresholds for statistical significance was set to P < 2.8×10−5 (Bonferroni method). The NAFLD PRS was used as a continuous exposure variable in a logistic regression adjusting for age, sex, age-adjusted AUDIT-C, and 10 principal components in an additive effects model using the PheWAS R package in R v3.2.065. The results from these analyses are reported as odds ratios, in which the estimate is the average change in odds of the PheWAS trait per NAFLD-increasing polygenic risk score.
Transcription Factor Analysis
We identified nominated genes (Supplemental Table 34) that encode for TFs based on known motifs, inferred motifs from similar proteins, and likely sequence specific TFs according to literature or domain structure115. Target genes for these TFs were extracted using DoRothEA database116 in OmniPath collection117 using the associated Bioconductor R package OmnipathR118, a gene set resource containing TF-TF target interactions curated from public literature resources, such as ChIP-seq peaks, TF binding site motifs and interactions inferred directly from gene expression.
Data Availability
The full summary level association data from the trans-ancestry, European, African American, Hispanic, and Asian meta-analysis from this report will be available through dbGAP project number phs001672.v4.p1 (study-specific accession codes will be available before publication).
https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/analysis.cgi?study_id=phs001672.v4.p1
MESA/MESA SHARe Acknowledgements
MESA and the MESA SHARe projects are conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with MESA investigators. This research was supported by R01 HL071739 and MESA was supported by contracts 75N92020D00001, HHSN268201500003I, N01-HC-95159, 75N92020D00005, N01-HC-95160, 75N92020D00002, N01-HC-95161, 75N92020D00003, N01-HC-95162, 75N92020D00006, N01-HC-95163, 75N92020D00004, N01-HC-95164, 75N92020D00007, N01-HC-95165, N01-HC-95166, N01-HC-95167, N01-HC-95168, N01-HC-95169, UL1-TR-000040, UL1-TR-001079, UL1-TR-001420. Also supported in part by the National Center for Advancing Translational Sciences, CTSI grant UL1TR001881, and the National Institute of Diabetes and Digestive and Kidney Disease Diabetes Research Center (DRC) grant DK063491 to the Southern California Diabetes Endocrinology Research Center.
Ethics statement
The Central Veterans Affairs Institutional Review Board (IRB) and site-specific Research and Development Committees approved the Million Veteran Program study. All other cohorts participating in this meta-analysis have ethical approval from their local institutions. All relevant ethical regulations were followed.
Data availability
The full summary level association data from the trans-ancestry, European, African American, Hispanic, and Asian meta-analysis from this report will be available through dbGAP (Accession codes will be available before publication).
Disclosures
H.R.K. is a member of a Dicerna scientific advisory board; a member of the American Society of Clinical Psychopharmacology’s Alcohol Clinical Trials Initiative, which during the past three years was supported by Alkermes, Amygdala Neurosciences, Arbor Pharmaceuticals, Dicerna, Ethypharm, Indivior, Lundbeck, Mitsubishi, and Otsuka; and is named as an inventor on PCT patent application #15/878,640 entitled: “Genotype-guided dosing of opioid agonists,” filed January 24, 2018.
Acknowledgements
This research is based on data from the Million Veteran Program, Office of Research and Development, Veterans Health Administration and was supported by award no. MVP000. This publication does not represent the views of the Department of Veterans Affairs, the US Food and Drug Administration, or the US Government. This research was also supported by funding from: the Department of Veterans Affairs awards I01-BX003362 (P.S.T. and K.M.C) and I01BX003341 (H.R.K. Co-Principal Investigator) and the VA Informatics and Computing Infrastructure (VINCI) VA HSR RES 130457 (S.L.D). B.F.V. acknowledges support for this work from the NIH/NIDDK (DK101478 and DK126194) and a Linda Pechenik Montague Investigator award. K.M.C, S.M.D, J.M.G, C.J.O, L.S.P, and P.S.T. are supported by the VA Cooperative Studies Program. S.M.D. is supported by the Veterans Administration [IK2 CX001780]. Funding support is also acknowledged for MS (K23 DK115897), R.M.C (R01 AA026302), D.K. (National Heart, Lung, and Blood Institute of the National Institutes of Health [T32 HL007734]), L.S.P. (VA awards I01 CX001025, and I01 CX001737, NIH awards R21 DK099716, U01 DK091958, U01 DK098246, P30 DK111024, and R03 AI133172, and a Cystic Fibrosis Foundation award PHILLI12A0). The Rader lab was supported by NIH grants HL134853 (NJH and DJR) and DK114291-01A1 (K.T.C, N.J.H, and D.J.R). We thank all study participants for their contribution. Support for imaging studies was provided by ITMAT (NIH NCATS UL1TR001878), the Penn Center for Precision Medicine Accelerator Fund and R01 HL137501. Data for external replication and hepatic fat concordance were provided by investigators using United Kingdom BioBank, Multi-Ethnic Study of Atherosclerosis (MESA), Old Order Amish Study (Amish), Framingham Heart Study (FHS) and Penn Medicine Biobank (PMBB).
References
- 1.↵
- 2.
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.
- 14.↵
- 15.
- 16.↵
- 17.
- 18.
- 19.
- 20.
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.
- 66.
- 67.
- 68.↵
- 69.
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵
- 108.↵
- 109.↵
- 110.↵
- 111.
- 112.↵
- 113.↵
- 114.↵
- 115.↵
- 116.↵
- 117.↵
- 118.↵