Transcriptional regulation of PNPLA3 and its impact on susceptibility to nonalcoholic fatty liver Disease (NAFLD) in humans

eQTL mapping

Studying genotypic data of 620,900 SNPs as well as imputed genotypic data (in total over 3.9 million SNPs using HapMap release 27 as reference), we focused our analysis on identification of cis-acting eQTLs on chromosome 22. We found that the most significant eQTLs were mapped to the PNPLA3-SAMM50 region (Fig. 1A-B). Two SNPs, rs139051 (51bp upstream of rs738409 in intron 2, and 28bp away from the intron-exon boundary) and rs2294918 (E434K in exon 9) were identified as the most significant cis-acting eQTLs for PNPLA3 transcription [p = 6.6×10^-8 and 2.4×10^-7, respectively; false discovery rate (FDR) < 0.05 for both] (Fig. 2A-B).

Detailed linkage disequilibrium (LD) analysis of variants across the locus revealed that whilst rs139051 was in low LD with all other significant eQTLs (r² ≤ 0.50 for all tests), rs2294918 was in high LD with the remaining eQTLs (r² ≥ 0.80), suggesting an independent role of rs139051. The plotted data indicated that the A allele of rs139051 and G allele of rs2294918 were associated with increased PNPLA3 gene expression (Fig. 2A, Fig. 2B). In contrast, after correcting for multiple testing, the rs738409 polymorphism was not a significant eQTL for hepatic PNPLA3 expression (FDR > 0.05) (Fig. 2C).

We further analyzed the inter-dependence of rs139051, rs2294918 and rs738409 in altering PNPLA3 gene expression. Analysis of the HapMap data (CEU, Northern Europeans in Utah) demonstrated that, although there is low LD among these three SNPs as measured by the r² values (r² ≤ 0.26 for all), the D’ values (coefficient of linkage disequilibrium) (D’=0.61, 0.88 and 1 between rs139051 and rs738409; between rs139051 and rs2294918; and between rs738409 and rs2294918, respectively) suggested alleles of these SNPs may share the same haplotype(s). Similar data (low r² but relatively high D’ values) were observed in both Asian (CHB+JPT, Han Chinese in Beijing and Japanese in Tokyo) and African (YRI, Yoruba in Ibadan, Nigeria) populations (data not shown). The phased HapMap haplotypic data (CEU) suggest that the rs139051A and rs2294918G alleles (both of which are associated with increased PNPLA3 expression) are almost always co-segregated in the same haplotype (data not shown). Further multivariate haplotypic association with gene expression in liver tissue demonstrated that the association between the rs2294918 G allele and PNPLA3 expression was actually driven by the rs139051 A allele. In rs139051 G/G homozygous individuals, rs2294918 genotypes did not significantly contribute additional variability in gene expression, which was in contrast with the association between rs139051 genotypes and gene expression among the rs2294918 A/A homozygous individuals (Fig. 2D). Similarly, the association between rs139051-rs738409 diplotypes and the PNPLA3 mRNA expression also suggested that rs139051 solely contributed to gene expression (Fig. 2E). The independent role of rs139051 in association with gene expression was further verified in a multivariate analysis in which, after controlling for rs139051 genotype, none of the other eQTLs at the locus remained significant (p > 0.05 for all tests, data not shown), while after conditioning all other SNPs, rs139051 remained significant (p = 8.18×10^-6). Taken together, these data suggest a potential causal role of rs139051 in regulating hepatic PNPLA3 transcription.

Validation of the rs139051 as a significant eQTL

To confirm the association between rs139051 and PNPLA3 mRNA expression, we used another set of liver samples among which genome-wide transcriptome data has been previously collected in a subset of samples (N = 54) [22]. These samples were collected from bariatric surgery of obese patients in Caucasian origin [22]. We found that again, rs139051 rather than rs738409 was significantly associated with PNPLA3 transcription (p = 0.037 for rs139051; p = 0.95 for rs738409). After controlling for age, gender and BMI, this association at rs139051 remained significant (p = 0.029).

Electrophoresis mobility shift assay (EMSA)

In order to further explore the potential causal role of rs139051, we performed an EMSA to examine the potential interaction between the SNP sequence and hepatic nuclear factors. We found that rs139051 exhibited an allele-specific binding to nuclear proteins extracted from the HepG2 cell line, with the rs139051 A allele possessing significantly higher (> 2 fold) affinity to nuclear factor(s) compared to the G allele (p = 0.017). This allelic binding specificity to nuclear protein was supported by a much weaker signal in the reaction with proteins extracted from the cytoplasm as compared to the nuclear extracts (Fig. 3A and Fig. 3B), suggesting that this binding is between the DNA and a nuclear-specific protein rather than a cytoplasmic protein that may non-specifically bind to DNA.

Association between SNPs and NAFLD

To test whether rs139051 or rs2294918 contributes to NAFLD risk, we performed association analyses between rs139051, rs2294918 as well as rs738409 and NAFLD phenotypes in a large cohort of European Caucasian NAFLD patients with histologically characterized disease collected by the FLIP Consortium (N = 949). Consistent with previous studies [23], rs738409 was significantly associated with histological disease severity across key quantitative phenotypes including degree of steatosis (p = 1.1×10^-13), NAFLD Activity Score (NAS, the sum of the scores for steatosis, hepatocyte ballooning degeneration and lobular inflammation and reflects disease activity) (p = 1.5×10^-8) and stage of fibrosis (p = 7.7×10^-5) after adjusting for relevant covariates (age, BMI and presence of type 2 diabetes mellitus). A less strong but statistically significant association between rs139051 and steatosis (p = 1.8x10^-4), steatohepatitis (p = 2.1x10^-4) but not fibrosis (p = 0.435) was also observed. Similarly, rs2294918 was also significantly associated with steatosis (p = 2.1x10^-5), steatohepatitis (p = 2.8x10^-4) and fibrosis (p = 0.01). However, in the multivariate analysis conditioning on the rs738409 genotype, neither rs139051 nor rs2294918 was significantly associated with these phenotypes (p > 0.05 for all tests), suggesting that neither rs139051 nor rs2294918 is independent of rs738409 as a risk factor for NAFLD.

To investigate further whether rs139051 or rs2294918 modifies the effect of rs738409 alleles in contributing to NAFLD risk, we tested associations between the rs139051-rs738409 or rs2294918-rs738409 haplotypes and aforementioned phenotypes. However, analyses examining the relative contribution of each SNP to the overall haplotype effect did not identify a statistically significant independent effect of rs139051 or rs2294918 in the FLIP population. Indeed, the carriage of the A-G haplotype which combines both the high-expression rs139051 A allele and the disease-risk rs738409 G (148M) allele, conferred similar risk for NAFLD phenotypes as compared to the G-G haplotype which contains the rs738409 risk allele, but possesses a relatively lower expression level (p = 0.51) (Table 1 for rs139051 and Table 2 for rs2294918).

To further verify these observations in an independent sample, we performed similar analysis in the aforementioned bariatric patient cohort (N = 212) in which quantitative pathological evaluation of NAFLD phenotypes has previously been described [22,24,25]. Given the dependence of rs2294918 on rs139051 as a significant eQTL, we did not include rs2294918 in the following analyses. Again, rs738409 but not rs139051 was significantly associated with NAFLD in the univariate analysis (p = 2.0×10^-4 and 0.093, respectively), and no significant independent association between rs139051 and NAFLD phenotypes was observed after conditioning on rs738409 and adjusting for age, gender and BMI. Also, no significant difference between A-G and G-G haplotypes in conferring risk for NAFLD or influencing NAFLD histological features was observed (Table 3). In a subset of 54 samples in which the PNPLA3 mRNA level was measured, we examined if the PNPLA3 mRNA expression is a direct co-factor modifying the genetic effect of rs738409 on NAFLD. We found that the PNPLA3 mRNA level was not significantly associated with any of the NAFLD phenotypes (p > 0.05 for all tests) (data not shown). The PNPLA3 mRNA level did not significantly modify the risk of NAFLD conferred by rs738409 (main effect of PNPLA3 mRNA level on NAFLD p = 0.61; interaction between rs738409 and PNPLA3 mRNA level with respect to NAFLD (p = 0.97).

In order to expand our observations to other ethnic groups, we further tested the genetic association in a Han Chinese NAFLD patient (N = 384) and healthy control (N = 384) population with matched age and gender that are reported in our previous study [26]. After adjusting for age and BMI, both rs738409 and rs139051 were individually associated with sonographically detected NAFLD (p = 2.0×10^-4 and 1.3×10^-3, respectively). Once again, after controlling for rs738409, rs139051 was no longer associated with NAFLD (p = 0.22). In testing the association between rs139051-rs738409 haplotypes and NAFLD, we again found that the A-G and G-G haplotypes confer similar risk of NAFLD (p = 0.27) (Table 3).

Given that the rs738409 G allele was consistently demonstrated to be significantly associated with increased risk for NAFLD phenotypes in all three cohorts, the rs738409 C allele was correspondingly significantly associated with reduced risk. We investigated whether rs139051 modifies this protective effect of the rs738409 C allele. Again, no statistical significance was found between the A-C and G-C haplotypes in conferring the reduced risk for NAFLD in all three cohorts (p > 0.2) (Table 1, 3-4).

Ethical statement

This study consists of samples obtained from several sources and parent studies as described below. Samples were collected with written informed consent obtained from participants with approval of institutional review boards (IRBs) at each participating institution. The Purdue University IRB has approved the use of the liver samples for the purpose of this study (approval number 1202011870). The study was carried out in compliance with the Helsinki Declaration.

Study populations

FLIP samples: A cohort of 949 unrelated patients of European-Caucasian descent with histologically characterized NAFLD representing the full spectrum of disease from simple steatosis through steatohepatitis to advanced fibrosis and cirrhosis was assembled from centers across Europe as part of a GWAS study by members of the FLIP consortium. The study had the necessary ethical approvals in all recruiting centers and all participants gave informed consent. Alternative diagnoses were excluded, including excess alcohol intake, chronic viral hepatitis (hepatitis B and hepatitis C), autoimmune liver diseases, hereditary hemochromatosis, α1-antitrypsin deficiency, Wilson’s disease and drug induced liver injury.

Bariatric NAFLD samples: The samples (N = 212) were collected under a study protocol that had been previously reviewed and approved by the Medical College of Wisconsin's Institutional Review Board. Each participating subjected gave a written Informed consent for participation. Subjects were of northern European descent, morbidly obese (BMI ≥ 40 kg/m² or > 35 kg/m² with significant comorbidities) and prior unsuccessful attempts to lose and maintain weight, who underwent bariatric surgery. A protocol intraoperative liver biopsy was performed on all patients for histological phenotyping. Patients with alcohol intake >20 g/d and those with other liver diseases based on positive disease-specific serological tests and suggestive liver histology were excluded. Patients using drugs associated with NAFLD prior to liver biopsy were excluded [22,24]. An experienced pathologist read the liver biopsy according to the NIH NASH Clinical Research Network working group system [34]. Genome-wide transcriptome analysis in the 54 samples in this study was performed previously [22].

Han Chinese samples: Blood samples for ultrasonography-defined NAFLD patients (N = 384) and healthy controls (N = 384) were collected in the health examination center of Shanghai Jiaotong University Affiliated Hospital (Shanghai, China) based on a previously study approved by the Shanghai Jiaotong University Affiliated Hospital [26]. The NAFLD diagnosis was according to the guideline defined by the Chinese National Consensus Workshop on Nonalcoholic Fatty Liver Disease [35]. Detailed information for this community-based cohort was described previously [26]. Briefly, physical examination, ultrasound screening, medical history review, questionnaire-based survey for diet, smoking and alcohol intake, as well as biochemical laboratory tests were performed. Cases with known causes of steatosis, e.g. heavy alcohol intake (> 20 g/day), the use of medications/herbals known to contribute to hepatic steatosis, as well as hepatitis B and C virus infection were excluded. Patients with high likelihood to have other known liver diseases including autoimmune disease and primary biliary cirrhosis based on aforementioned information were also excluded. The case and control groups were matched for age and gender with both consisted of 229 males and 155 females, as well as a mean age and standard deviation (SD) of 45 ± 13 years.

The study populations were summarized in Table 5.

eQTL mapping for PNPLA3

The liver sample set (N = 206) used to perform eQTL mapping has been described in our previous study [21] (Table 5). Briefly, genotyping and transcriptome profiling on these samples were performed on the Illumina Human 610 quad beadchip platform (GPL8887) and the Agilent 4x44 arrays, respectively (GEO accession: GSE28893). Genotype imputation was conducted with IMPUTE2 program [36] after prephasing the genotypes with SHAPEIT [37] using the HAPMAP release 27 as a reference to obtain genotypes for a total of over 3.9 million SNPs. Standard quality control (QC) process was performed by assessing SNP/sample call rates, minor allele frequency (removal of SNPs with frequency < 15%), deviation from Hardy-Weinberg equilibrium (p < 0.001), identity-by-decent analysis, sex imputation for each sample, etc., as previously described [21,38]. Age, gender and the first 3 components of the principal component analysis (PCA) for quantifying the ancestry of samples were included as covariates for eQTL mapping. We also employed the probabilistic estimation of expression residuals (PEER) framework [39] to quantifying unknown hidden factors. To perform the PNPLA3 gene eQTL mapping, we focused only on the cis-acting eQTLs within the ±1Mb region of PNPLA3 gene. The eQTL mapping was performed using the --assoc function for quantitative traits in PLINK [40], which correlates allele dosage with changes in the trait. False discovery rate (FDR) was used for adjusting the multiple testing [41]. Plotting of the p values of SNPs was conducted using LocusZoom [42].

Genotyping

Genotyping of rs738409 and rs139051 in DNA samples were performed using Taqman assays according to the instructions provided by the manufacturer (Life Technologies, CA, USA).

Electrophoresis Mobility Shift Assay (EMSA)

EMSA assay was performed according to our protocol published previously [43]. Briefly, HepG2 cells were cultured in standard conditions and collected for both nuclear and cytoplasmic protein extraction using a NE-PER Nuclear and Cytoplasmic Extraction Kit commercial kit (Thermal Scientific, IL, USA). EMSA was performed with 2.5µg of total nuclear or cytoplasmic extracts using a LightShift Chemiluminescent EMSA Kit (Thermal Scientific). Biotinylated probe sequences for rs139051 A and G alleles were Biotin/GGGTGCTCTCGCCTATAACTTCTCTCTCCT and Biotin/GGGTGCTCTCGCCTGTAACTTCTCTCTCCT, respectively. Non-biotin labeled probes in a higher concentration (200×) than the labeled probes were used as competitors. EMSA bands were visualized by exposed to the CL-XPosure film (Thermal Scientific) for 3 min. Films were then scanned and the intensity of the band interested was quantified using the ImageJ software [44]. EMSA assays were independently performed three times.

Haplotype estimation and diplotype assignment

Haplotype frequency between rs139051 and rs738409 among the eQTL mapping population as well as the diplotype for each individual were estimated and assigned by using PHASE [45] under default settings. To plot the data, samples with the diplotype assignment probability lower than 90% were excluded.

Data analysis

Linkage disequilibrium and haplotype analyses between rs139051, rs738409 and rs2294918 in HapMap samples (CEU, HCB+JPT and YRI) were performed using Haploview 4.2 [46]. Comparisons of gene expression level between genotypes were conducted using one-way ANOVA, and the comparison for allelic binding efficiency in EMSA assays was conducted with unpaired t-test based on the data collected from three independent experiments. Association between diabetes and NAFLD or PNPLA3 mRNA expression and interaction between polymorphisms and diabetes was examined using linear (for quantitative phenotypes) or logistic (for binary phenotypes) regression analysis within the R package.The genetic association between SNPs or haplotype and NAFLD or NAFLD-related phenotypes (including the effect of adding one SNP to a model that already included the effect of the other SNP) was analysed using the software programs R, PLINK [40] and UNPHASED [47,48]. The same programs were used for the multivariate test for inter-dependence between rs139051 and the remaining significant eQTLs. Univariate comparisons of gene expression between genotypes/diplotypes were performed using a t-test or ANOVA test. Data plotting was carried out using Graphpad Prism 6.00 (Graphpad software Inc, CA, USA).