Epigenomic map of human liver reveals principles of zonated morphogenic and metabolic control

A deeper epigenomic understanding of spatial organization of cells in human tissues is an important challenge. Here we report the first combined positional analysis of transcriptomes and methylomes across three micro-dissected zones (pericentral, intermediate and periportal) of human liver. We identify pronounced anti-correlated transcriptional and methylation gradients including a core of 271 genes controlling zonated metabolic and morphogen networks and observe a prominent porto-central gradient of DNA methylation at binding sites of 46 transcription factors. The gradient includes an epigenetic and transcriptional Wnt signature supporting the concept of a pericentral hepatocyte regeneration pathway under steady-state conditions. While donors with non-alcoholic fatty liver disease show consistent gene expression differences corresponding to the severity of the disease across all zones, the relative zonated gene expression and DNA methylation patterns remain unchanged. Overall our data provide a wealth of new positional insights into zonal networks controlled by epigenetic and transcriptional gradients in human liver.

T issues are composed of heterogeneous cell populations, each contributing to the overall performance of the biological system. Within these populations, even morphologically similar cells of the same type may perform different functions. Deciphering this cellular heterogeneity on a functional level is key to understand tissue function and its alterations in disease. Here, functional genomics may play an important role if applied in a tissue context 1 , as combined transcriptional and epigenetic analysis provides deep insights into functional state and genomic programming of cells [2][3][4] . Such regionalized genomic maps may be generated by single cell analysis of isolated cells and subsequent mapping on tissue architecture using marker genes 5,6 or by analysis of cells obtained directly from tissue context using laser capture cryo microdissection (LCM). Advances in the respective protocols allow the generation of transcriptomes 7 and reduced representation bisulfide sequencing (RRBS) methylomes 8 from the 100 to 200 cells obtained typically by LCM. Such analysis could provide a combined epigenetic dataset, preserve tissue context, avoid isolation procedures and may thus offer new opportunities for the analysis of shock frozen human tissue.
Despite their rather uniform morphological appearance, hepatocytes in mammalian liver are characterized by a heterogenous enzyme distribution along the porto-central axis 9,10 which is termed metabolic zonation 11,12 . Recently, using single cell mRNA sequencing and spatial reconstruction, a global signature of zonation in mouse liver has been established 5 . Immunohistochemical analyses have confirmed the zonal distribution of metabolic enzymes, such as cytochrome P450 13,14 , alcohol dehydrogenase 15 , fatty acid-binding protein 16 , glucose-6phosphatase and glutamine synthetase 17 in human liver. Relevant questions, such as the link of metabolic zonation to underlying morphogens, a global assessment of zonation and the underlying regulatory mechanisms and not least also the relation of rodent findings to human liver 18 are as yet not fully addressed 19 .
In an effort to contribute to a better understanding of these problems, we generated a spatially resolved transcript and DNA methylation map of human hepatocytes, demonstrate an epigenetic zonation in human liver and provide pathway and transcription factor signatures along the porto-central axis that help to understand the underlying programs of liver metabolism, regeneration and tissue structure.

Results
Laser microdissection of hepatic zones from human biopsies. Hepatocytes were isolated by LCM from surgical liver samples from 19 patients. Control liver samples (NC, N = 4) were obtained from patients undergoing liver resection for metastasis surgery or surgical liver biopsy for the exclusion of malignancy. N = 15 liver samples were obtained from patients undergoing bariatric surgery and were categorized based on standardized liver histology employing the NAS score 20 into phenotypic groups of healthy obese (HO, N = 5), bland steatosis (STEA, N = 5) and early non-alcoholic fatty liver disease (EARLY, N = 5). An overview of the phenotypic characteristics of the patients is provided in Table 1. Liver tissue was shock frozen directly in the operating room ensuring ex vivo times of less than 40 s in all cases. Hepatocytes were isolated by LCM from the pericentral, intermediate and periportal areas of each patient from cryosections stained by cresyl violet (Fig. 1a, Supplementary Figure 1). RNA quality and quantity was assessed with a Bioanalyzer (Agilent, Santa Clara, CA) yielding a mean RIN of 8.09 and a mean yield of 5.43 ng (Supplementary Figure 2A). Samples were subjected to a modified SMART-Seq RNA sequencing protocol 7,21 . Comparisons of the expression of indicator transcripts from cell types, such as immune cells, endothelial cells and fibroblasts indicated an insignificant contamination by nonparenchymal cells (Supplementary Figure 2B, C). Adjacent cryosections were used to generate matching reduced representation bisulfite sequencing (RRBS) data using the CpG methylation insensitive restriction enzymes AluI and HaeIII to cover~14 million CpGs per sample 8 .
Signatures of transcriptional and epigenetic zonation. Principal component analysis of the 1000 most variable transcripts (Fig. 1b) revealed highly significant correlations of PC1 and PC2 with hepatic zones (PC1: r = 0.47, p = 2.0 × 10 −4 ; PC2: r = 0.86, p < 2.2 × 10 −16 , Supplementary Figure 3) indicating zonation as a major driver of variation between samples. Phenotypic group, sex, age, BMI and diabetes state were used as covariates for paired differential analysis using edgeR (FDR < 0.01, |log2FC| > 1). Using the full dataset, we identified 317 genes predominantly expressed in the pericentral zone and 488 genes with predominant periportal expression (Supplementary Data 1). We did not observe genes with highest expression in the intermediate zone, but rather transcriptional gradients from pericentral to intermediate to periportal zone or vice versa (Fig. 1c). Upon analysis of differential expression of steatotic samples (STEA, EARLY) versus samples with normal liver histology (NC, HO), a significant deregulation of 467 genes was observed (Supplementary Figure 4A, Supplementary Data 2), however largely irrespective of hepatic zones (Supplementary Figures 4B, 5). If zonation analysis was restricted to patients with normal liver histology (NC, HO), reduced power led to the detection of only 470 zonated differentially expressed genes (DEGs, Supplementary Data 3). Comparison of pairwise effect sizes between the analyses revealed highly similar log 2 fold changes (all p < 2.2 × 10 −16 ) with pairwise correlations (r 2 ) ranging from 0.80 to 0.98 ( Supplementary Figure 6) indicating that the transcriptional zonation signature is persisting during early phases of NAFLD (Fig. 1b).
The transcriptional zonation of the mouse liver lobule has recently been systematically studied from single cell transcriptomic data defining 9 zones 5 . For an interspecies orthologue comparison of human and mouse data, we merged zones 1-2 of the mouse for a comparison to our human pericentral data, the mouse zones 4-5 to the intermediate zone and the mouse zones 8-9 to the periportal zone. We obtained a global transcript correlation of r 2 = 0.77 (p < 2.2 × 10 −16 ) in periportal and pericentral zones and r 2 = 0.78 (p < 2.2 × 10 −16 ) in the intermediate zone (Supplementary Figure 7A-C) indicating an overall The number of patients in each phenotypic category is provided together with demographic and histologic characteristics. All numeric traits are shown as the median with the range provided in parenthesis. BMI, body mass index; NAS, NAFLD activity score 20 comparability of the datasets between species. Our zonation comparison was restricted on 313 genes zonated in humans that featured one to one matching orthologues with high confidence expression values in the mouse dataset. We detected 111 genes with conserved zonation signatures between mouse and human (Supplementary Data 4) including highly expressed zonation landmark genes, such as GLUL, CYP2E1, SDS and APOF and morphogens like LGR5 and SOX9 with low transcript levels (Fig. 1d). However, 202 genes zonated in humans (e.g. Wnt pathway members such as RSPO2) were not detected as zonated in the mouse. Most non-corresponding genes in the mouse showed a rather low expression level as compared to our deeply sequenced human data (Supplementary Figure 7D). Comparison of available immunohistochemistry (CYP1A2, EPHX1, EPHX2 and FABP) and in situ hybridization (ALB, CYP1A2) data from strongly expressed human genes revealed a similar zonation profile 14,16,17,22 (Fig. 1e). On the epigenetic level, a principal component analysis (Fig. 2a) of the 5000 most variable CpGs (coverage >5 in all samples) showed a separation by hepatic zones with highly significant correlation of PC1 (p = 6.7 × 10 −4 , r = 0.44, Supplementary Figure 4) and zonation. Following covariate adjustment, a differential methylation analysis (see method for details) between pericentral and periportal samples (MethylKit, CpG coverage >10, 500 bp tiles, at least 3 CpGs, FDR adjusted p-value <0.01) yielded 17,862 differentially methylated regions (DMRs) with at least 5% methylation difference (Supplementary Data 5). In contrast we identified only 23 significant DMRs exhibiting the highest methylation difference between PV or CV and IZ. The vast majority of DMRs featured methylation gradients from pericentral to intermediate to periportal zone (Fig. 2b). The majority of DMRs was detected in gene context and only 18.7% of DMRs were located in intergenic regions (Fig. 2c).
To integrate DNA methylation with expression changes, DMRs were annotated to the closest gene transcriptional start site. 44.3% (357) of zone-specific differentially expressed genes were associated with zonal DMRs (1094). A total of 205 DEGs were related with multiple DMRs. DMRs associated to differentially expressed transcripts were predominantly located in introns (46.4%), exons (18.7%) and promoter regions (24.7%, Fig. 2c). DNA methylation differences of significant DEG-associated DMRs and corresponding expression changes ( Fig. 2d) were largely anti-correlated (p < 2.2 × 10 −16 , r 2 = −0.37). Ninety nine DEGs associated with multiple DMRs displayed a combination of hypo-and hypermethylated patterns of DMRs. We identified 132 upregulated genes associated with 430 hypomethylated DMRs in the pericentral zone and 139 upregulated genes with 291 hypomethylated DMRs in the periportal zone. Based on this clear anti-correlation between DNA methylation and expression gradients, we define these 271 transcripts as epigenetically marked driver genes of zonation (Supplementary Data 6). An example of such potential epigenetically marked driver genes is the Wnt signalling pathway member Palmitoleoyl-Protein Carboxylesterase (NOTUM) containing a cluster of six hypomethylated promoter DMRs in the CV zone where the gene is exclusively expressed. Publicly available DNaseI-seq and ChIPseq data from the hepatocyte-like cell line HepG2 (ENCODE) identifies potential regulatory activity for all seven DMRs (Supplementary Figure 8).
DNA methylation control at transcription factor binding sites. Tissue structure and zonation may be maintained not only by external morphogen gradients, but also by differential signalling efficacy of transcription factors, that may even be uniformly present along the porto-central axis: For instance, the hepatic nuclear factor 4a (HNF4A) is a highly expressed transcription factor (logCPM = 6.12 in our dataset) binding to promoters of 12% of liver expressed genes 23 . How HNF4A exhibits is differential activity on zonated genes, such as cytochrome P450 and aldehyde dehydrogenase genes 23 remains unclear. Our RNAseq data-in line with previous reports 24 -show that HNF4A expression shows only a marginal trend towards zonation (FDR PV_CV = 0.02, Log 2 FC PV_CV = -0.19). Still Hnf4a deletion in mice leads to a clear zonal disruption of key metabolic enzymes 25 . Our genome-wide DNA methylation data combined with the ENCODE maps of transcription factor binding sites in HepG2 cells 26 , allowed us to assess a potential epigenetic layer of transcription factor binding control: Indeed 46 of the 59 available transcription factors showed an enrichment of differential methylation at their mapped binding sites in the CV or PV zone ( The strongest effect size and a pericentral hypomethylation (Δ methylation = 7%, p = 1.2 × 10 −80 ) was observed for TCF7L2, the key downstream transcription factor of the Wnt signalling pathway 28 . Preferred binding of TCF7L2 at hypomethylated sites indicates higher pericentral signalling efficiency (Supplementary Data 8). An additional analysis using previously identified TCF7L2 binding sites 29 further supported our zonation results (Supplementary Figure 9). TCF7L2 interacts with HNF4α in multiple hepatocyte specific functions 30 . Like HNF4α, overall TCF7L2 (FDR PV_CV = 0.008, Log 2 FC PV_CV = 0.3) shows only a weak trend towards zonation. Because multiple TCF7L2 protein products with differential activity have been reported 31 , we analyzed alternative splicing of TCF7L2 using the SALMON parallel inference algorithm 32 . No zonation dependent alternative splicing was detected (Supplementary Figure 10), thus indicating that differential binding efficiency rather than transcript abundance underlies the zonation of this factor. In order to provide a broader assessment of zonated TF binding site methylation, we performed an enrichment analysis using predicted TF sites (Supplementary Data 9), which for instance also revealed a pericentral hypomethylation of HNF4A and TCF7L2 sites.
Zonated features of metabolic networks. Gene ontology overrepresentation analysis both on the DEG and DMR datasets (Supplementary Data 10, 11) recapitulated the proposed distribution of functional activities along the porto-central axis 10 . We find a clear enrichment of GO terms associated with metabolism of steroids, cytochrome P450, xenobiotics or retinoic acid in the pericentral area, glutathione and amino acid metabolism in the periportal area (Fig. 3d). Metabolic zonation is also seen in the GO analysis on the epigenetic level (Supplementary Data 4). A consistent pattern of hypomethylation and overexpression for key zonated metabolic enzymes (Supplementary Data 6) suggests that they are epigenetically marked driver genes controlling zonated regulatory pathways. In the pericentral area we find a strong expression and hypomethylation of CYP1A2 (logFC = −4.0, 12 DMRs: Δ meth 5.4-33.8%) and ALDH3A1 (logFC = −5.4, two DMRs: Δ meth 6.3-19.3%). The location of these enzymes corresponds to the pericentral location of liver damage in alcoholic liver disease and Paracetamol toxicity. Periportal candidate epigenetically marked driver genes include glutaminase 2 (GLS2, logFC = 3.0, 5 DMRs: Δ meth −15.1, −26.0%) and histidine ammonia-lyase (HAL, logFC = 5.1, 3 DMRs: Δ meth −17.2, −33.2%), which are relevant in the glutamate and amino acid metabolism.
When assessing the underlying layer of transcriptional regulation, strong zonation of DNA methylation of binding sites of transcription factors of metabolic relevance (e.g. HNF4A and HNF4G) was evident (Fig. 3a, Supplementary Data 7). Notably, methylation of binding sites for Retinoid X receptor alpha (RXRA), the heterodimeric partner for nuclear signalling of key metabolic transcription factors PPARA, LXRA and FXR 33 , is implications for ligands such as free fatty acids which play a key role in the pathogenesis of non-alcoholic fatty liver disease (NAFLD). Interestingly, FXR agonists compounds and drugs such as signalling biliary acids (e.g. obeticholc acid 34 ) might therefore elicit a zonal therapeutic response, which may allow targeting of the pharmacodynamic effect to the pericentral disease location in NAFLD.
We also observed a significant gradient of mitochondrial DNA copies deduced from the ratio of mitochondrial reads in the genomic coverage of RRBS data. There is significantly less mitochondrial DNA in the pericentral zone than in the intermediate and periportal zone (p = 4.1 × 10 −6 , Fig. 3c). This is in correspondence with a metabolic porto-central oxygen gradient 35 .
Morphogen gradients and putative stem cell niche. The transcriptional and epigenetic zonation includes a series of known morphogen factors (Fig. 3d, Supplementary Data 9, 10). In the periportal zone we find an enrichment for Notch signaling pathway related genes, while there was an enrichment of Wnt signalling among the pericentrally upregulated genes (e.g. GO:0016055, p = 5.7 × 10 −4 ). A key role of Wnt signaling for the maintenance of liver structure and regeneration is well established from murine liver 12,[27][28][29][30] . A recent analysis indicated that Wnt pathways are also involved in human liver zonation 19 . Our deep integrated datasets allowed a deeper look into the putative role of pericentral Wnt signalling networks in human liver homoeostasis: On both, the epigenetic and transcriptional level, we identify a composed network of regulatory changes affecting selected members of the pathway but not a global upregulation or activation (Fig. 4a). Prominent examples are the epigenetically marked driver loci AXIN2 (Δ methCV_PV = 5.0-29.3%), NOTUM (Δ methCV_PV = 10.7-36.5%), and LGR5 (Δ methCV_PV = 9.2%). In contrast, classical Wnt agonists such as WNT3A, WNT3 and genes of the export machinery (e.g. porcupine, PORCN1 and Wntless, WLS) are not zonally expressed, suggesting that the pericentral WNT gradient in hepatocytes might not primarily be regulated by them. Negative regulators of Wnt signaling such as NOTUM, NKD1, RNF43 and ZNRF3 are upregulated 1 to 2 log levels pericentrally creating a tightly controlled signalling environment around hepatocytes of this zone. The Wnt signalling enhancers RSPO2, RSPO3 and LGR5 are expressed 2-7 log levels higher pericentrally than periportally suggesting sensitization of pericentral hepatocytes for canonical Wnt signals 29 . Interestingly, LGR5 (logFC = 7.15, FDR = 4.88 × 10 −26 ) shows the highest fold change of all significantly DEGs. Finally, binding sites of the main Wnt-driven transcription factor (T-cell factor 7, TCF7L2, p = 1.22 × 10 −80 ) in the liver 28 are >5 fold pericentrally enriched among the pericentral DMRs (Fig. 3b, c). We conclude that morphogen gradients are controlled by multi-layered transcriptional and epigenetic Wnt signalling signatures in the zonated parenchyma.

Discussion
The systematic generation of reference epigenomes and transcriptomes through the ENCODE and IHEC projects allows deep insights into the functional state of isolated cells or tissues. The integration of cell-specific epigenomic and transcriptomic profiles with positional information appears as an essential next step to better understand the organisation of organs and tissues and to dissect regional program alterations in disease. Here, we present the first spatially resolved map of human hepatocyte genomic architecture by integrating methylome and transcriptome data. Previous human data is heterogenous and patchy and has used either measurements of enzyme activities 15,44,45 or immunohistochemistry and in situ hybridization assays targeting whole-gene families. We excluded the activity assays from our confirmatory comparisons in Fig. 1e, because enzyme activities reflect often multiple transcripts, cofactor concentrations, phosphorylation state and other confounders. While confirming previous immunohistochemical and gene expression studies 14,16,17,19,22 , we identify hundreds of additionally zonated transcripts and observe intriguing correlations between transcription and DNAmethylation.
We optimized a pipeline starting from fast sampling in the operating room, followed by a documented LCM on sectional stains, which were directly used for the generation of genome . Average CpG methylation differences between CV and PP range between 12.4% and 49.2%. Blue corresponds to low and red to high average methylation values. The annotation legend displays lean controls (NC) in light blue, healthy obese (HO) in dark blue, steatosis (STEA) in red and early NASH (EARLY) in orange. c Genomic annotation of all DMRs (upper pie plot) and of DEG-associated DMRs (lower pie plot). d Visualization of gene expression (purple) and DNA methylation (dark red) z-scores of DEG-associated DMRs and corresponding genes at pericentral (green), intermediate (red) and periportal zone (blue). In the Heatmap yellow corresponds to low z-scores (reduced expression or DNA methylation) and navy to high z-scores (increased expression or DNA methylation) wide transcriptome and RRBS data 8 from 100 to 200 captured cells per section. This approach is complementary to the alternative single cell and FACS sorting approaches as it provides greater sequencing depth and uses the pathological classification and direct zonal localisation within the tissue. Single cell approaches on the other hand while providing a cleaner view on cellular heterogeneity they have to rely on computational modelling of their spatial position using independent molecular signatures as anchor points. This highly informative single cell modelling usually comes at the expense of sequencing depth and gene coverage. We believe that our approach is an ideal complementation to single cell approaches as it provides strong  Fig. 3 Transcription factor, mitochondrial and metabolic zonation along the porto-central axis. a Analysis of methylation of transcription factor binding sites across zones. The y-axis depicts the normalized average methylation difference at binding sites overlapping with DMRs for the respective transcription factor. The x-axis depicts the log 10 -fold enrichment of DMRs among binding sites. A detailed list and the analysis method are provided in Supplementary Data 7. b Methylation of DMRs containing a transcription factor binding site for TCF7L2: Raw methylation data for all sites are shown for the CV and PV regions. c Boxplot of the ratio of mitochondrial reads in reduced representation bisulfite sequencing data (n = 19) in pericentral (green), intermediate (red) and periportal (blue) zone. Significant differences (p < 0.001) are marked by asterisks. d Visualization of significantly enriched GO terms (p < 0.05) of zonally expressed genes. A negative GO term z-score corresponds to GO enrichment in the pericentral zone (green), a positive z-score to periportal enrichment (blue). Circle size reflects the number of zonated genes relating to the term. Exemplary GO terms are labelled, while the full list of significant GO terms can be found in Supplementary Data 5 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06611-5 8 sequencing depth, positional information and allows to comprehensively cover the epigenetic dimension. Such positional datasets can be queried in multiple ways offering a new and integrated approach of epigenetic and transcriptional mechanisms controlling zonation, liver function and morphogen control: First, we demonstrate for example, that epigenetic zonation underlies the transcriptional porto-central gradient and that anti-correlated DNA methylation and transcription may define a set of potentially epigenetically marked driver genes, that form a deeper layer underlying the zonation gradients. Second, we propose that strongly zonated DNAmethylation at transcription factor binding sites may form an additional layer of epigenetic zonation control. Differential methylation of binding sites may indicate a spatial activity of otherwise uniformly expressed transcription factors, such as HNF4A and RXR. Drugs targeting RXR (e.g. obeticolic acid) or its interaction partner PPARα (i.e. fibrates) are likely to preferentially affect target genes in pericentral hepatocytes, which corresponds to the disease location of metabolic liver disease as a zone 3 disorder 20 . Third, the integrated genomic dataset also allows to assess the existing concepts on liver regeneration onto the human steady-state situation: We demonstrate a strong transcriptional and epigenetic gradient of Wnt signalling and modulating factors. The pericentral expression of LGR5 and AXIN2 and the corresponding gradient of the liver progenitor marker TBX3 46,47 indicate a pericentral source of hepatocyte regeneration 38 under steady-state conditions in humans. Conversely, in the portal area, we observed a NOTCH and EPCAM signature. This may correspond to hepatic progenitor cells at the border to the biliary tree, which may form a relevant source of regeneration under injury conditions 36,39,48 . Expansion of this compartment clearly needs to be confirmed in further studies using cirrhotic human liver. In our dataset, we see that zonation in both-the transcriptomic and epigenetic-dimensions is the by far dominant signature and remains largely preserved in steatosis and early NASH. This finding is compatible with the complete regeneration and remodelling capacity of human liver in early NAFLD.
In summary our study opens a new path for a functional dissection of parenchymal differences in the liver lobule connecting positional information to known aspects of liver biology. Many of the observations described will be helpful for a deeper understanding of liver organ function in relation to disease.

Methods
Human liver samples. Liver samples were obtained intraoperatively in patients in whom an intraoperative liver biopsy was indicated on clinical grounds such as during scheduled liver resection, exclusion of liver malignancy during major oncologic surgery, or assessment of liver histology during bariatric surgery. Standardized histology by a single pathologist (C.R.) employed the NAFLD activity score (NAS) 20 . Samples were frozen immediately in liquid nitrogen ensuring an ex vivo time of less than 40 s in all cases. Patients with evidence of viral hepatitis, hemochromatosis or alcohol consumption >20 g per day for women and 30 g per day for men were excluded. All patients provided written informed consent. The study protocol was approved by the institutional review board (Ethikkommission der Medizinischen Fakultät der Universität Kiel, D425/07, A111/99), before study commencement.
Laser capture microdissection. The frozen tissue was cut in 20 µm slices and mounted on precooled glass-slides. Immediately after mounting, the tissue was fixed by incubation in 70% EtOH for 2 min and stained by incubation in a cresyl violet 50% EtOH solution for 30 s with two short subsequent wash steps in 70% EtOH, followed by a 2 min air-drying step. All steps were performed on ice. In sum, an area of 500,000 µm 2 , from of 5 to 10 periportal fields and intermediate zones or hepatozytes surrounding 15-20 central veins per sample, was microdissected using a ZEISS PALM MicroBeam LCM in Auto-LPC mode. Tissue pieces for RRBS were frozen and maintained at −80°C. Samples for RNA-Seq were immediately processed.
RNA seq. After microdissection, RNA was isolated using the RNeasy Micro Kit (Qiagen, Hilden, Germany). Concentration and RNA Quality were obtained by using an Agilent BioAnalyzer (Agilent, Santa Clara, CA). The typical yield was between 4.5 ng and a RIN between 8 and 9. The RNA was reverse transcribed and amplified with the SMARTer Universal Low Input RNA Kit for Sequencing (Takara, Clontech) and 1 ng of RNA. Libraries for next generation sequencing were constructed with the Nextera DNA Library Prep Kit (Illumina, San Diego, Ca, USA) and sequenced on a HiSeq2500 Sequencer.
RNA-seq reads were aligned with the STAR (v2.5.2a) 55 and the per sample 2pass mapping strategy processing all reads in both passes as described in the documentation. Read counts were summarized to Gencode gene models (v19) 53 with featureCounts (v1.5.0-p3) 56 counting primary alignments only. EdgeR (v3.16.5) was used to detect differentially expressed genes with maximal q-value of 0.01 and minimal absolute logFC of 1. We discarded genes for which fewer than five samples had a counts per million value above 0.5, calculated normalization factors and robustly estimated the dispersion. In our default setting we used a model adjusting for phenotype, sex, age, BMI, donor and the prevalence of diabetes (Supplementary Figure 3).
Integrative analysis and visualization. For interspecies comparisons, reconstructed mouse zones 1 and 2 were considered as pericentral, zones 4 and 5 as intermediate and zones 8 and 9 as periportal. Spearman rank statistics of sequencing depth normalized mRNA expression levels (log(CPM + 1)) were calculated for mouse single cell data 5 and human LCM-RNAseq data with one-to-one matching orthologs. The zonation comparison was restricted to genes with high confidence expression values corresponding to a CPM > 0 in the mouse dataset after background subtraction as previously described 5 . For visualization of exemplary genes across the three zones, expression data were normalized to the zone with the highest expression for each individual gene. Potential epigenetically controlled driver genes were defined as upregulated genes that also featured at least NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06611-5 ARTICLE NATURE COMMUNICATIONS | (2018) 9:4150 | DOI: 10.1038/s41467-018-06611-5 | www.nature.com/naturecommunications one hypomethylated DMR annotated to the gene. Genomic tracks of expression and methylation data (bigwig) were generated using bedtools 54 and UCSC tools 57 and visualized in the IGV browser 58 . GO term enrichment analysis of DEGs and DMR associated genes was performed with DAVID 59 and significant DEG GO terms (p < 0.05) were visualized using GOplot 60 .
For analysis of methylation of transcription factor binding sites with available information of the hepatocyte cell line HepG2 in ENCODE (V3), a 2 × 2 table was constructed for each transcription factor as follows: The 196,621 methylation regions (500 bp intervals), were categorized in regions with and without binding sites for the respective transcription factor and by differential methylation between the CV and PV zones as analysed in the main DMR analysis (Supplementary Data 6) yielding four categories (noDMR.noTFBS, noDMR.TFBS, DMR.noTFBS and DMR.TFBS). Deviation from the null hypothesis (i.e. equal distribution of binding sites between DMR and regions without differential methylation) was tested by Fisher's exact test. The methylation gradient between CV and PV for intervals without and with binding sites for the respective transcription factors yields the zonated methylation difference Delta.TXF. Because all DMRs show a mild porto-central demethylation gradient and this gradient is reflected in DMRs with transcription factor binding sites, the zonated differential methylation is normalized to the mean of Delta.TXF (−1.27%) for the final normalized Delta. TXF.Norm.

Data availability
The data generated in this study is available at the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo) under identifier GSE105127.