Nonalcoholic fatty liver disease (NAFLD) is a metabolic disease in which fat accumulates in the liver without significant alcohol intake (Cholankeril et al., 2017). NAFLD includes both nonalcoholic fatty liver or simple steatosis, which is a relatively early stage of disease and low risk, and nonalcoholic steatohepatitis (NASH), which is an advanced condition with inflammation and fibrosis (Schuster et al., 2018). Due to the complicated and multifactorial mechanisms underlying the development and progression of NAFLD, our understanding of NAFLD remains insufficient.

Genetic alterations associated with NAFLD may play a role in the pathology thereof. Single nucleotide polymorphisms of PNPLA3, GCKR, and TM6SF2 are known genetic changes associated with NAFLD progression (Sliz et al., 2018). Genomic, transcriptomic, and epigenomic analyses of tissue samples from NASH patients have recently been reported and potential biomarkers for NAFLD have also been suggested (Govaere et al., 2020; Oh et al., 2020). Another study also compared the changes in DNA methylation throughout the genome to identify epigenetic markers associated with NAFLD severity (Loomba et al., 2018). These studies led to the identification of several biomarkers specific to various stages of NAFLD progression. However, they are still insufficient to accurately diagnose the stages of NAFLD progression, which is significantly affected by several environmental factors. Thus, epigenetic analysis of the NAFLD progression is essential for understanding the development and progression of the disease.

In recent years, the assay for transposase-accessible chromatin using sequencing (ATAC-seq) has enabled the genome-wide profiling of chromatin accessibility (Corces et al., 2017). Particularly for small amounts of frozen tissues, ATAC-seq has shown wide applicability; however, it has not yet been conducted for NAFLD. Here, we report a comprehensive epigenomic analysis using ATAC-seq technology across the full histological spectrum ranging from healthy controls to fibrotic NASH in an Asian biopsy-proven NAFLD cohort (Kim et al., 2017; Yoo et al., 2021). This provides a valuable resource that can be used to elucidate transcriptional regulation. This study may also help to understand NAFLD caused by regulatory dysfunction.

Sample collection

The use of human liver samples was approved by the Institutional Review Board (IRB) of SMG-SNU Boramae Medical Center (IRB No. 16-2014-86). All participants were informed of the study protocol and provided written and signed consent. The analyzed population consisted of 19 male and 25 female study subjects, aged 19-80, who visited SMG-SNU Boramae Medical Center. Liver histology was assessed using the NASH Clinical Research Network (CRN) histological scoring system (Sanyal et al., 2021). NAFLD was defined as the presence of >5% macrovesicular steatosis. We graded steatosis, lobular inflammation, and ballooning according to the NAFLD activity score (Kleiner et al., 2005). NASH was determined based on an overall pattern of histological hepatic injury consisting of steatosis, lobular inflammation, and ballooning according to the criteria of Brunt et al. (1999; 2011). Fibrosis was assessed according to a 5-point scale proposed by Brunt and modified by Kleiner et al. (2005): F0, absence of fibrosis; F1, perisinusoidal or periportal fibrosis; F2, perisinusoidal and portal/periportal fibrosis; F3, bridging fibrosis; and F4, cirrhosis (Kleiner et al., 2005). Then, we categorized all study subjects into healthy control, steatosis, and fibrotic NASH according to the histological spectrum of NAFLD. Liver tissues were obtained from 44 individuals through percutaneous needle liver biopsy. Each sample was quickly frozen in liquid nitrogen and stored until nuclei extraction was performed.

The isolation of nuclei and the transposition for ATAC-seq

The flash-frozen tissues were manually dissociated. Isolated nuclei were quantified using a hemocytometer, and 100,000 nuclei were tagmented as previously described (Buenrostro et al., 2013), with some modifications (Corces et al., 2017) using the enzyme and buffer provided in the Nextera Library Prep Kit (Illumina, USA). The tagmented DNA was then purified using the MinElute PCR Purification Kit (Qiagen, Germany), amplified with 10 cycles of polymerase chain reaction (PCR), and purified using Agencourt AMPure SPRI beads (Beckman Coulter, USA). The resulting material was quantified using the KAPA Library Quantification Kit for Illumina platforms (KAPA Biosystems, USA), and sequenced with PE42 sequencing on the NextSeq 500 sequencer (Illumina).

ATAC-seq and computational analysis

ATAC-seq library construction was processed. All ATAC-seq libraries were sequenced on the Illumina NextSeq 500 with 42 bp paired-end reads. The reads were aligned to hg38 genome using bowtie2-2.3.4.1 (Langmead and Salzberg, 2012) with parameters “-k 4 --end-to--end”. The reads with duplicates and aligned to mitochondrial DNA were excluded using SAMtools-1.7 (Li et al., 2009) for further analysis. To correct a bias caused by tn5 transposase, all mapped reads were offset by +4 bp for the + strand and –5 bp for the – strand (Buenrostro et al., 2013). To build bigwig track, genomeCoverageBed and bedGraphToBigWig were used. MACS2-2.1.2 (Zhang et al., 2008) was used to call peaks for each sample with parameters “-q 0.05 --nomodel --shift -100 --extsize 200” and peaks on ENCODE blacklist v2 regions were excluded (Amemiya et al., 2019).

Identification of differentially accessible regions (DARs)

DiffBind (Ross-Innes et al., 2012) was used to select differentially accessible peaks with FDR (false discovery rate) < 0.05. The de novo motif analysis was performed with DREME and best matches with JASPAR 2018 CORE motifs (Khan et al., 2018) were found using TOMTOM in MEME-suite version 5.0.5 (Bailey et al., 2015). GREAT (genomic regions enrichment of annotations tool) was used to link sets of genomic regions to putative biological functions based on the functional annotations of the nearby genes (McLean et al., 2010).

Data deposition and code availability

The NGS sequence data generated in this paper have been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (https://dataview.ncbi.nlm.nih.gov/object/PRJNA725028?reviewer=535ltoopds44f3nkklr49alha4).

The code used for data analysis is available at https://github.com/sysgenlab/Steatosis.

Chromatin accessibility landscape of NAFLD progression

To determine the chromatin accessibility of NAFLD progression, ATAC-seq was performed on liver biopsy samples (Fig. 1A, mentary Tables S1 and S2). Liver tissues were obtained from 44 individuals, consisting of 4 healthy controls, 19 simple steatosis patients, and 21 fibrotic NASH patients, through percutaneous needle liver biopsy. The clinical features are presented in Table 1. The clinical information included sex, age, body mass index (27.6 ± 3.5 kg/m², mean ± SD), diabetes mellitus, grade of steatosis, steatosis score, lobular inflammation, portal inflammation, hepatocyte ballooning degeneration, NAFLD activity score, gamma-glutamyltransferase, fasting blood sugar, homeostatic model assessment of insulin resistance, cholesterol, high-density lipoprotein, triglyceride, and hemoglobin A1c.

Data quality control from the pipeline with the fraction of reads in peaks quality control

We evaluated our ATAC-seq dataset using a series of commonly used statistics, including the number of total aligned pairs, overall alignment rate, de-duplicated pairs, non-mitochondrial, non-redundant fraction, PCR bottleneck coefficient 1, PCR bottleneck coefficient 2, total reads, and the fraction of reads in peaks (FRiP) (Table 2). For each sample, we obtained an average of 50 million reads. The library complexity was measured using the non-redundant fraction, PCR bottleneck coefficient 1, and PCR bottleneck coefficient 2. The successful detection of accessible regions is also supported by the observation of strong enrichment of ATAC-seq reads around transcription start sites (+/– 5 kb, “promoters”), gene bodies (with 2 kb flanking regions), and merged peak regions (equating to all peak regions; +/– 5kb) (Figs. 1B and 1C). ATAC-seq signals are enriched at promoters and merged peak regions, indicating the acceptable quality of the ATAC-seq data. Additionally, the chromatin accessibility fragments showed mainly a mononucleosome peak (Fig. 1D). Together, these results show that high-resolution chromatin accessibility can be obtained from low amounts of frozen tissue samples. We identified accessible regions by using the MACS2-2.1.2 algorithm. We then applied the FRiP quality control and filtered out low-quality data, which was defined as FRiP scores that were less than 20%. Thus, we defined 14 profiles as valid datasets for downstream analysis.

Evaluation of the ATAC-seq data sets

The similarity within each group is shown by heatmap clustering of Pearson correlation coefficients from the comparison of 14 (>20 FRiP) ATAC-seq profiles (Fig. 2A). The correlation heatmap generated by the cross-correlation of two samples based on the read counts of all merged peaks showed a similar pattern. Steatosis and fibrotic NASH were divided into separate clusters. This result was also supported by principal component analysis (PCA) (Fig. 2B). PCA classified steatosis and fibrotic NASH along the PC1 axis. The correlation coefficient between samples within the steatosis group was up to 0.97 and the fibrotic NASH group showed a slightly higher correlation coefficient (Figs. 2A and 2C). The high correlation of the ATAC-seq signal between each of the three stages demonstrated high reproducibility. The ability to cluster samples by their disease stage also showed significantly different chromatin accessibility between these three NAFLD stages. The sequencing track of liver fibrogenesis signature genes (Col1a1, Acta2) (Elbadawy et al., 2020; Seki et al., 2007) showed a strong signal in fibrotic NASH samples (Fig. 2D). Thus, these data suggest that our ATAC-seq analyses could detect accessible chromatin regions in the NAFLD patient genome and could be used as a reference dataset for future studies.

Identification of DARs

To select differentially accessible peaks, DiffBind protocol, originally intended for ChIP-seq data, which constructs a consensus read count matrix from MACS2, was used. When 10 fibrotic NASH samples were compared with 3 steatosis samples, 272 up-regulated and 1,137 down-regulated DARs were identified (Figs. 3A-3C). Functional enrichment analysis was performed on 1,137 down-regulated regions. Significant terms of enriched Gene Ontology were listed in Fig. 3D. Functional enrichment analysis of more closed regions in fibrotic NASH confirmed that these regions were near genes related to hormone response, cell cycle, and liver development. The other DAR results are presented in Supplementary Tables S3-S8.

Transcription factor binding motif and de novo motif analyses

To further characterize these DARs, we performed the transcription factor binding site enrichment and de novo motif analyses on the 272 regions that were up-regulated in fibrotic NASH and the 1,137 regions that were down-regulated in fibrotic NASH when compared to steatosis (Fig. 4). As NAFLD progresses from steatosis to fibrotic NASH, factors such as E2F6, ELF3, RBPJ, ETV4, IKZF1, OSR2, ETV1, ZKSCAN5, GABPA, SPIC, and EWSR1-FLI1 may play an important role in the regions where chromatin is open, and factors such as ETV1, EHF, GABPA, IKZF1, ELF3, ELF5, SPI1, ELF1, ETV4, STAT2, ZNF384, and IKZF1 in the regions where chromatin is closed. These factors could be promising diagnostic biomarkers that accurately discriminate the stages of NAFLD progression.

In the present study, ATAC-seq analysis was conducted to define the changes in chromatin accessibility in liver tissues from healthy controls, steatosis, to fibrotic NASH. Our human liver ATAC-seq data generated for different stages of the NAFLD cohort also suggested open chromatin regions were specifically up- or down-regulated in fibrotic NASH tissues compared to steatosis tissues.

Analysis of high-quality ATAC-seq data obtained from biopsy samples of 44 patients with NAFLD revealed 112,303 merged peak regions (Fig. 1, Supplementary Table S1). These chromatin accessibility profiles identified stage-specific DNA regulatory elements that allow the classification of NAFLD subtypes with newly recognized prognostic importance. We also found 1,648 DARs, showing differences in chromatin openness among the three groups of patients (Fig. 3A, Supplementary Tables S3-S8). In addition, we determined 23 new biomarker candidates through transcription factor (TF) binding motif enrichment and de novo motif analysis (Fig. 4, Supplementary Tables S9 and S10).

For example, among the candidates, the chromatin accessibility of EWSR1-FLI1 motif-enriched regions was specifically upregulated in fibrotic NASH (Fig. 4B, upper). EWSR1-FLI1 has an enhancer-regulating role in Ewing sarcoma through divergent chromatin remodeling mechanisms (Riggi et al., 2014). An integrated analysis of chromatin states of Ewing sarcoma showed that EWS-FLI1 could function as a pioneer factor in generating enhancers de novo at repeat elements or inhibiting conserved enhancers by competing with endogenous erythroblast transformation specific (ETS) factors. The ETS transcription factor family might be functionally relevant for NAFLD progression since ESE3 overexpression suppressed hepatocellular carcinoma (HCC) (Lyu et al., 2020). ETV4 is positioned in the regulatory module of HCC progression (Kim et al., 2018). The ETS family belongs to candidate TFs whose binding motifs are enriched during the disease progression (Fig. 4A). In addition, SPIC could be a transcription factor that plays a role in the activation of a hepatic stellate cells, which is critical in fibrosis (Marcher et al., 2019). Another attractive candidate is GABPA since it appears in the gene list of both open and closed chromatin. In a previous report by Lake et al. (2016), it was shown that the decreased expression of GABPA may cause down-regulation of its target genes in NASH. It is suggested that these factors could be promising diagnostic and therapeutic biomarkers.

Identifying open chromatin regions is the first step in understanding the regulatory function of non-coding regions, which are not revealed in gene expression changes. However, chromatin accessibility is not always correlated with transcription because histone modifications and DNA methylation status as well as an access to DNA also affect both transcription initiation and its output. Thus, it is necessary to integrate chromatin accessibility with other epigenetic data sets, such as whole-genome bisulfite sequencing and histone modification data (Choi et al., 2020; Gelbart et al., 2009; Shia et al., 2006). The integration of RNA-seq and ATAC-seq will help to find factors whose expression is sufficient for both motif protection and nucleosome positioning. This may enable the development of a quantitative model that links the accessibility of regulatory elements to the expression of the predicted target gene.

In conclusion, the present study will help us understand the epigenetic regulation of the NAFLD progression. Future efforts to combine the genome with the epigenome will pave the way for dealing with the non-coding genome in NAFLD.

This research was supported by the Collaborative Genome Program for Fostering New Post-Genome Industry of the National Research Foundation (NRF) and funded by the Ministry of Science and ICT (MIST) (NRF-2017M3C9A6044199) (R.H.S.). This work was supported by the National Research Foundation of Korea grant funded by the Korea government (MEST) (2021R1A2C2005820, 2021M3A9E4021818) (W.K.).

B.K. (Byeonggeun Kang), R.H.S., and W.K. conceived the idea. W.K. collected samples. B.K. (Byeonggeun Kang) generated the data. B.K. (Byeonggeun Kang) and B.K. (Byunghee Kang) analyzed the data. B.K. (Byeonggeun Kang) wrote the manuscript and T.-Y.R., R.H.S., and W.K. supervised the study and revised the manuscript. All authors reviewed and approved the final manuscript.

The authors have no potential conflicts of interest to disclose.

Clinical information

Data QC from the pipeline with FRiP quality control