Dynamic 3D genome reorganization during development and metabolic stress of the porcine liver

Liver development is a complex process that is regulated by a series of signaling pathways. Three-dimensional (3D) chromatin architecture plays an important role in transcriptional regulation; nonetheless, its dynamics and role in the rapid transition of core liver functions during development and obesity-induced metabolic stress remain largely unexplored. To investigate the dynamic chromatin architecture during liver development and under metabolic stress, we generated high-resolution maps of chromatin architecture for porcine livers across six major developmental stages (from embryonic day 38 to the adult stage) and under a high-fat diet-induced obesity. The characteristically loose chromatin architecture supports a highly plastic genome organization during early liver development, which fundamentally contributes to the rapid functional transitions in the liver after birth. We reveal the multi-scale reorganization of chromatin architecture and its influence on transcriptional regulation of critical signaling processes during liver development, and show its close association with transition in hepatic functions (i.e., from hematopoiesis in the fetus to metabolism and immunity after birth). The limited changes in chromatin structure help explain the observed metabolic adaptation to excessive energy intake in pigs. These results provide a global overview of chromatin architecture dynamics associated with the transition of physiological liver functions between prenatal development and postnatal maturation, and a foundational resource that allows for future in-depth functional characterization.

Supplementary Tables S1-2 Table S1 Functional description of the genes mentioned in this study.

Table S2
Information of primers and genomic locations for target gene promoters and enhancers validated using the Dual-Luciferase reporter assay. Supplementary Fig. S8 Chromatin accessibility between prenatal (E80) and adult (2Y) stages. a Data summary of ATAC-seq (n = 6, three replicates for two selected stages, i.e., prenatal E80 and adult 2Y). b Scatterplots showing correlations of normalized ATAC-seq signals (reads per million, RPM) between samples. The genome-wide ATAC-seq signal was highly reproducible between biological replicates (average Pearson's r = 0.95), but were more dissimilar between stages (E80 and 2Y) (average Pearson's r = 0.83). c Heatmaps depicting the enrichment of normalized ATAC-seq signal cantered on transcription start site (TSS) and ordered by signal intensity. d, e Distribution of ATAC-seq peaks (d) and stage-specific ATAC-seq peaks (e) relative to genomic features. The peaks were classified into five distinct categories based on 1 bp overlap with features (i.e., promoter, exon, intron, transcription end site [TES], and intergenic regions). f Normalized ATAC-seq signal in each 20 kb bins for compartment A (orange) or B (blue) regions. Statistical significance was calculated by the Wilcoxon rank-sum test. g Average of normalized ATAC-seq signal (RPM) in a 4 kb region cantered on stage-specific peaks between E80 and 2Y. h Overlap of stagespecific ATAC-seq peaks or rest peaks with stage-restricted compartment A regions. Statistical significance was determined by two-proportions test. i Top ten enriched transcription factor binding motifs within stage-specific ATAC-seq peaks, as identified by JASPAR database 108 . Fig. S9 Rewiring of PEIs dynamically regulates gene expression in a developmental-dependent manner. a Evaluation of resolutions for the intrachromosomal Hi-C maps generated by the merged Hi-C data of six replicates at each stage and that of HFD-fed pigs. The curve represents the proportion of bins that contain at least 1,000 intra-chromosomal contacts in a relevant bin size. Map resolution is defined as the smallest bin size where 80% of the bins have at least 1,000 reads in order to allow reliable discerning of local features 89 . b Size distribution of PEIs identified at 5-kb resolution using the merged Hi-C data. PEIs larger than 20 kb were retained. The line represents the cumulative distribution. c Number of PEIs. d Proportion of enhancers interacting with the nearest promoter or skipping at least one promoter. e Proportion of PEIs within or across TADs. f Proportion of gene promoters from each expression category interacting with zero to more than six enhancers. g Gene expression was associated with RPS. Genes with RPS > 0 were divided equally into five different percentiles. Data are presented as mean ± SD. h Expression foldchanges for genes with differential RPS between consecutive stages. P-values were calculated using the Wilcoxon rank-sum test. i Determination of enhancer activity by analysing the distribution of H3K27ac signals using the ROSE algorithm 111 . j Proportion of genes interacting with super-enhancers (SEs), regular-enhancers (REs), and poised-enhancers (PEs). k RPS values and l expression levels of genes interacting with SEs, REs and PEs. m Expression levels of genes with active promoters (ac-P, enriched by H3K4me3 peaks) or inactive promoters (in-P). n Proportion of active or inactive promoters interacting with SEs, REs and PEs. o Validation of the enhancers identified in this study using a Dual-Luciferase Reporter Assay System in HEK-293T cells. Data shown are mean ± SD (n = 3). P-values were calculated using a two-sided Student's t-test.

Supplementary Fig. S10
Genes with the developmental-dependent changes in RPS are related to functional transitions during liver development. a Top ten significantly enriched GO-BP terms (performed by the Metascape 50 ) for five representative clusters of genes identified by the STEM algorithm based on dynamic RPS. b PEI rewiring of 26 genes that were enriched in the GO-BP term 'monocarboxylic acid metabolic process'. The changes in RPS (left), interacting enhancers (middle) and gene expression (right) are shown. c-e Examples of PEI rewiring for HCC signature genes 23 , ADH4 (c), FGL2 (d), and POSTN (e). ADH4 is a signature gene for a less aggressive HCC subtype (S-I, good prognosis) and is involved in drug metabolism. FGL2 and POSTN are signature genes for a more aggressive HCC subtype (S-III, poor prognosis) and participate in leukocyte activation. Left: Schematics of PEIs, H3K4me3 and H3K27ac signals, and transcription. Dashed boxes highlight SE peaks. Right: 3D structural models and Hi-C contact maps of the corresponding genomic region (seven 5-kb bins up-and down-stream from the bin containing the TSS). Gene promoters (blue spheres) and enhancers (red spheres) are shown in the 3D models. Fig. S11 STEM clustering of RPS profiles across developmental stages for eight a priori gene sets related to core liver functions at prenatal (hematopoiesis) and postnatal stages (metabolism of amino acid, fatty acid, glucose, bile acid, and drug; tricarboxylic acid cycle; and immunity). Genes in each set were classified into five or ten profiles based on RPS changes during development. The bold purple lines in each plot represent the mean RPS, and the grey lines represent the RPS of genes in relevant cluster during development. FDRcorrected P-values were obtained from multiple hypothesis testing. Identification (a) and functional enrichment (b) of genes with differential expression (|log 2 FC| > 1 and FDR < 0.05) between the two groups. The enrichment analysis was conducted using Metascape 50 . We show the top ten significantly enriched terms. c Differences in expression, RPS, and contacting enhancers of the nine genes exhibiting concomitant changes in expression and RPS. d-g PEI organization, histone modifications, and expression levels of four NAFLD markers, including CYP2E1 (d), IL6 (e), LEP (f), and TNF (g). Left: a schematic representation of PEIs, H3K4me3 and H3K27ac signals, and transcription. Right: 3D structural models of the corresponding genomic regions. Gene promoters and enhancers are shown as blue and red spheres, respectively.

Supplementary Fig. S14
Promoter-promoter interaction profiles responding to liver metabolic stress in pigs. a Comparison of promoter-promoter interactions in liver of normally and HFD fed pigs. Top panels: 3D spatial distance among genes from different expression categories in liver of normally fed (2Y, left) and HFD fed (right) pigs. Middle panels: Overlap of the top 20% of genes with the highest expression level in liver of normally fed (2Y) and HFD fed pigs. Bottom panels: Transcription factor (TF) enrichment for the genes with specific high expression in liver of normally fed (2Y, left) and HFD fed (right) pigs. The pig genes were converted to human orthologs. TFs enriched in the specifically expressed genes were identified using the section "TRRUST Transcription Factors 2019" from the Enrichr web server (http://amp.pharm.mssm.edu/Enrichr). The top ten statistically significant TFs were depicted in bar plots. b Gene expression of predicted target genes of HNF4α and C/EBPα. P-values were calculated using a Wilcoxon rank-sum test.

Supplementary Fig. S15
The loose chromatin architecture allows for the transcription of extensive genomic regions during early liver development. a 3D genome structures of porcine livers during development, which were inferred using the normalized intra-(20-kb resolution) and inter-chromosomal (1-Mb resolution) contact maps with the miniMDS 96 and visualized with PyMOL (v 2.5.2). The nuclear radius is determined by the average distance to the nuclear center of mass. b Comparison of genomic features between regions that have distinct compartment status (A or B) during liver development, including compartment length, number and density of genes, enhancer density, as well as H3K4me3 ChIP-seq, H3K27ac ChIP-seq and ATAC-seq signals. c Expansion of the transcriptionally active compartment A regions facilitate widespread active transcription during early liver development. By comparing the representative stages between prenatal E38 and adult 2Y stages, we observed that 246.78 Mb (10.89% of the genome) and 88.60 Mb (3.91% of the genome) specifically exhibited compartment A status in E38 and 2Y, respectively. The genes located in these stage-restricted compartment A regions (1,668 for E38, and 648 for 2Y) tended to show increased gene expression. We then compared compartment scores (i.e., the A-B index values) for the regions with the same compartment status between the two stages. This allowed us to identify regions spanning 870.12 Mb (38.40% of the genome) and 476.50 Mb (21.03% of the genome) that exhibited statistically significant elevated compartment scores (P < 0.05, unpaired Student's t-test, and |△A-B index| > 0.5) in E38 and 2Y, respectively. The genes located in these stage-specific higher compartment score regions (7,269 for E38, and 4,991 for 2Y) also tended to show increased gene expression. The genes with evidence of transcription (TPM > 0.5) in E38 or 2Y stages were used for expression comparison.

Supplementary Fig. S16
Human trait-associated noncoding SNPs were enriched in enhancers in the porcine liver. a Comparison of trait-associated, noncoding SNP enrichments in enhancers and other genomic regions in the porcine liver across six developmental stages. The identification of SEs, REs, and PEs was based on H3K27ac binding peaks. SNP enrichment score calculation can be found in Materials and methods. b Heatmap of noncoding SNP enrichment scores in enhancer regions in the porcine liver across six developmental stages. The merged traits or diseases with the top ten highest SNP enrichment scores at each stage (P < 0.05,  2 test) are displayed. Traits or diseases can be generally classified into three categories: hematopoiesis-related traits, lipid metabolism-related traits, and others. The numbers in parentheses indicate the counts of trait-associated, noncoding SNPs. c Radar plots showing differences in SNP enrichment linked to traits of core liver functions, including hematopoiesis (upper) and lipid metabolism (lower), in enhancers during six developmental stages of the porcine liver. The colored dots on the axis indicate SNP enrichment scores in SEs (red points) or the union set of all enhancers (blue points) at each stage. The numbers in parentheses indicate the counts of trait-associated, noncoding SNPs. Supplementary Table S1 Functional description of the genes mentioned in this study.

Gene name Related figure in manuscript
Gene function Reference

HBB, HBE1
Hemoglobin subunit beta, Hemoglobin subunit epsilon 1 These genes encode 146-aminoacid globin chains, which compose hemoglobin in combination with αlike globin to participate in oxygen transport. A protein kinase, which phosphorylates and inactivates pyruvate dehydrogenase complex that catalyzes a rate-limiting step of glucose oxidation. Its upregulated expression is associated with fatty acid oxidation in liver and heart. It is highly expressed in proerythroblasts and erythroblasts in the fetal liver and encodes a transcription factor that stimulates erythroid cell survival, proliferation, and terminal maturation while silences epsilon globin expression in definitive erythropoiesis. GHR is a key regulator of postnatal growth and is closely related to metabolism. GH binds to GHR, followed by the activation of the JAK-STAT pathway and subsequent increase in expression of IGF1.

CPS1
Carbamoylphosphate synthase 1  PPARα is a member of PPAR family, which is abundantly expressed in the liver and is essential for modulation of lipid transport and metabolism through transcriptionally activated genes that are involved in mitochondrial and peroxisomal fatty acid βoxidation. The absence of PPARα enhances hepatic steatosis.  G6PC is a key enzyme in the maintenance of glucose homeostasis, which catalyses the hydrolysis of glucose-6-phosphate to glucose and phosphate in the terminal step of gluconeogenesis and glycogenolysis. The mutations in G6PC are associated with Glycogen storage disease type Ia.  C3 is a member of complement system and plays a central role in the complement cascade, which is essential to the innate immune process. Thorgersen