{'fraction_non_alphanumeric': 0.030958283927993963, 'fraction_numerical': 0.017872156947288994, 'mean_word_length': 5.097003154574132, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
3D disorganization and rearrangement of genome provide insights into pathogenesis of NAFLD by integrated Hi-C, Nanopore, and RNA sequencing

The three-dimensional (3D) conformation of chromatin is integral to the precise regulation of gene expression. The 3D genome and genomic variations in non-alcoholic fatty liver disease (NAFLD) are largely unknown, despite their key roles in cellular function and physiological processes. High-throughput chromosome conformation capture (Hi-C), Nanopore sequencing, and RNA-sequencing (RNA-seq) assays were performed on the liver of normal and NAFLD mice. A high-resolution 3D chromatin interaction map was generated to examine different 3D genome hierarchies including A/B compartments, topologically associated domains (TADs), and chromatin loops by Hi-C, and whole genome sequencing identifying structural variations (SVs) and copy number variations (CNVs) by Nanopore sequencing. We identified variations in thousands of regions across the genome with respect to 3D chromatin organization and genomic rearrangements, between normal and NAFLD mice, and revealed gene dysregulation frequently accompanied by these variations. Candidate target genes were identified in NAFLD, impacted by genetic rearrangements and spatial organization disruption. Our data provide a high-resolution 3D genome interaction resource for NAFLD investigations, revealed the relationship among genetic rearrangements, spatial organization disruption, and gene regulation, and identified candidate genes associated with these variations implicated in the pathogenesis of NAFLD. The newly findings offer insights into novel mechanisms of NAFLD pathogenesis and can provide a new conceptual framework for NAFLD therapy.


Introduction
Non-alcoholic fatty liver disease (NAFLD) is a common chronic liver disorder associated with accumulation of fats in the liver. The prevalence of NAFLD is high in obese individuals, and increases susceptibility to liver fibrosis, cirrhosis, and hepatocellular carcinoma, which is considered a major health problem worldwide 1,2 . Genetic factors are believed to contribute to 30%e50% of the risk for high-prevalence diseases such as NAFLD 3 . Efforts in understanding the molecular mechanisms of NAFLD initiation and progression have elucidated the significant role of genetic factors in disease susceptibility, with more information still emerging. Gene variants like polymorphism rs626283 in membrane bound O-acyltransferase domain containing 7 (MBOAT7) gene may affect glucose metabolism by modulating intra-hepatic fat content in Caucasian obese children and adolescents, increasing the risk of NAFLD as observed by Umano et al. 4 , the I148M patatin like phospholipase domain containing 3 (PNPLA3) variant has been identified as the major common genetic determinant of NAFLD 5 . Moreover, genome-wide association studies have uncovered robust and reproducible associations between variations in genes such as transmembrane 6 superfamily member 2 (TM6SF2), glucokinase regulator (GCKR), hydroxysteroid 17-beta dehydrogenase 13 (HSD17B13), protein phosphatase 1 regulatory subunit 3B (PPP1R3B), immunity related GTPase M (IRGM), lipin 1 (LPIN1), and the natural history of NAFLD 6,7 . In view of the complexity of NAFLD pathogenesis, the variation of individual gene cannot play a decisive role in the occurrence of disease. Analysis of full genome and disclosure of NAFLD candidate gene genes will improve understanding of NAFLD pathophysiology.
RNA-sequencing (RNA-seq) and whole-genome sequencing (WGS) are now widely applied technologies to obtain deep insights into diseases. Major advances in sequencing technologies have enabled implementation of RNA-seq and WGS in determining the molecular pathogenesis of liver diseases 8,9 . Nanopore sequencing as a third-generation sequencing approach, promises high sequencing throughput and longer read length 10,11 , allowing a comprehensive analysis of the entire genome with identification of large numbers of structural variations (SVs) and copy number variations (CNVs). Nanopore sequencing is widely used in "-omic" sciences like genomics, epigenomics, and transcriptomics to evaluate the mutational status of single or multiple genes involved in disease processes 11 .
SVs are commonly found in the genome, however their functional impacts are often elusive. SV events may be accompanied by alternations in the chromosomal three-dimensional (3D) structures disrupting the gene structures, leading to disease 12 . Tissue-specificity and development require maintenance of precisely regulated gene expression programs in higher eukaryotes 13 . This tissue-specific gene regulation is closely associated with the highly folded and organized 3D conformation of the genome. This 3D organization of the chromatin plays critical roles in many processes like regulation of DNA replication, DNA damage repair, X chromosome inactivation 14e16 . Alternations in the 3D genome may affect numerous key biological functions and their mechanistic connections with disease development are increasingly apparent 17 . Several studies about the 3D cancer genome have been performed in the past few years. It has been uncovered that about 12% of genomic regions had A/B compartment switch in breast cancer compared with normal cells, and gene-rich chromosomes displayed a lower interaction frequency with each other than in normal cells 18 . In lung squamous carcinoma, deletion of the topologically associated domain (TAD) boundary resulted in the spread of active chromatin and subsequent insulin receptor substrate 4 (IRS4) overexpression 19 . In some tumor suppressor genes, including but not limited to phosphatase and tensin homolog (PTEN), tumor protein p53 (TP53), RB transcriptional corepressor 1 (RB1), notch receptor 1 (NOTCH1), and neurofibromin 1 (NF1), sequences or promoters can be frequently interrupted by translocation and inversion, leading to inactivation and the subsequent onset and progression of cancer 20,21 . Changes in 3D conformation of genomes have also been reported in other disease such as heart failure, diabetes, and brain and chronic kidney diseases [22][23][24][25][26] .
The development and application of imaging and chromosome conformation capture (3C) during the last decade, have advanced our capacity to explore the structures in 3D genome architecture 27 . High-throughput chromosome conformation capture (Hi-C) sequencing technology is a 3C-based method, which obtains 3D genome information of a disease by detecting changes in chromatin interactions and hierarchical structures 28 , like alterations of A/B compartment, TAD boundary, and chromatin loop. By understanding of the role of 3D genome organization in the precise regulation of gene expression, the relationship between NAFLD pathogenesis and 3D genome needs to be elucidated.
In the present study, we performed an integrated analysis combing Hi-C, Nanopore sequencing, and RNA-seq data on mice with NAFLD, to better understand the molecular mechanisms of NAFLD. We investigated the spatial organization and alterations in the 3D genome interactions to reveal the relationship between 3D genome reorganization, SVs, and gene regulation in NAFLD. The integration of these technologies has the potential to bridge the gap between genome structures and gene function, to obtain deeper insights into NAFLD.

Animals
Male C56BL/6 mice (18e22 g) obtained from Liaoning Changsheng Biotechnology Co., Ltd. [SCXK (Liao) 2015-0001, Benxi, China] were housed in the controlled environment at 23 AE 2 C under a 12 h darkelight cycle with free access to food and water. After being acclimatized to laboratory conditions for 1 week, animal maintenance and experiments were performed. Mice were randomly divided into control and model groups (n Z 10). Mice in control group were fed a normal chow diet (12% kcal fat content), and mice in model group were fed a high-fat diet (HFD, 45% kcal fat content) for 8 weeks. Then, the serum samples and liver tissues were obtained for the following pathological staining and assay. Animal experiments were performed referring to the ARRIVE guidelines and approved by the Animal Research Committee of Dalian Medical University in accordance with the Laboratory's Animal Ethics Guidelines [SYXK (Liao) 2014-0002, Dalian, China].

Biochemical analysis and pathological staining
The serum levels of alanine aminotransferase (ALT), aspartate aminotransferase (AST), total cholesterol (TC), triglyceride (TG), and the levels of free fatty acid (FFA), malondialdehyde (MDA), superoxide dismutase (SOD), and glutathione (GSH) in liver tissue (n Z 10) were examined by the detection kits according to their protocols (Nanjing Jiancheng Institute of Biotechnology, Nanjing, China). At the end of the experiment, liver tissues were removed and fixed in 4% polyoxymethylene solution for more than 24 h. Then, the livers were embedded in paraffin wax and section for hematoxylin-eosin (H&E) staining to evaluate the pathological variation. The frozen liver sections were stained with oil red O to observe lipid droplets in the liver.
Three liver samples in each group were randomly chosen and simultaneously detected by RNA-seq, Hi-C, Nanopore sequencing, and qPCR methods.

Hi-C data processing
Three liver samples of mice without NAFLD (control group, C) and with NAFLD (model group, M) were mixed equally, respectively. The Hi-C library was prepared following the standard in situ Hi-C protocol 29 , using Dpn II as the restriction enzyme. The ends of restriction fragments were labeled using biotinylated nucleotides, and then ligated in a small volume. After reversal of crosslinks, ligated DNA was purified and sheared to a target size of 300e700 bp fragments, and then ligation junctions were pulled down with streptavidin beads and prepped for high-throughput Illumina sequencing. Two Hi-C libraries were established for each group. After construction of the library, the concentration of the library and the size of the inserted fragment (Insert Size) were detected by Qubit 2.0 (Agilent Technologies, Santa Clara, USA) and Agilent 2100 (Agilent Technologies), respectively. The effective concentration of the library was quantified to ensure the quality of the library. After passing the library inspection, high-throughput sequencing was carried out on the Illumina platform (Illumina, Santiago, USA), and the sequencing read length was PE150. Clean reads were obtained by removing reads containing adapter, poly-N, and low quality reads from raw data. Clean paired-end reads were mapped to the sequences of the reference genome using Burrows-Wheeler Aligner (BWA, 0.7.10-r789 version) 30 . The mapped reads were analyzed using HiC-Pro (v2.10.0) 31 , and the valid interaction pairs and invalid interaction pairs were identified. The corresponding standardized interaction matrix was obtained by HiC-Pro v2.10.0, and the Pearson correlation coefficient was calculated. Then, the resolution of Hi-C data was calculated based on the reported method 29 , and the Hi-C resolution is 10 kb in this study. Then, the Hi-C data are standardized by HiC-Prov 2.10.0 software with 100 kb, and the relationship between interaction frequency and genome linear distance is calculated by HiCdat (HiCdatR_0.99.0 version) 32 . The slop of the curve is the interaction decay exponents (IDEs). Instructure simulation of the whole genome, the mds algorithm of PASTIS software was used to simulate the 3D location of chromatin 33 .

A/B compartment analysis
We used normalized interaction matrices at 100 kb resolution to detect chromatin compartment types by HiTC (v1.24.0) software 34 . Positive or negative values of the first principal component separate chromatin regions into two spatially segregated compartments. The compartment with higher gene density and positive principal component analysis (PCA) value was assigned as A compartment, and the other compartment and negative PCA value was assigned as B compartment.

TAD analysis
We used TadLib (hitad 0.1.1-r1) software to identify TADs at 20-kb resolution 35 . Inclusion ratio (IR) value was calculated using HOMER (v4.10.1) software 36 . The TADs were obtained by filtering TAD with IR less than 0 and length less than 5 bin.

Chromatin loops analysis
Chromatin loops were identified with improved HICCUP method 32 . For the Hi-C matrix, HICCUP calculates several enrichment values over different local neighborhoods, and each enrichment value is associated with a false discovery rate (FDR) value for assessing statistical significance. We identified loops at a matrix resolution of 10 kb, only retaining loops that have an FDR no more than 0.1 every neighborhood. Then, we combined the loop results of each sample, and the original interaction frequency was used as the input of edgeR (v3.8.6) 37 . The P value was calculated by glmLRT function and adjusted by Benjamini-Hochberg. When the obtained loops with FDR less than 0.01, we recognized it as a differential loops.
2.7. Nanopore library preparation and sequencing DNA libraries were prepared using the protocol provided by Oxford Nanopore Technologies (ONT, Oxford, UK). Briefly, 1 mg isolated DNA was sheared to w8 kb fragments with g-tube using 4000Âg and 1 min centrifuging parameters. Then, the sheared products were purified using 0.45-fold XP beads according to Agencourt protocol (Beckman Coulter, Brea, USA). DNA fragments were then subjected to FFPE DNA repair and end-repair (NEB, Ipswich, USA) steps. Then DNA was ligated to the adaptor using T4 DNA ligase (NEB). Following adaptor ligation, the products were purified by adding a 0.6-fold Agencourt XP beads according to the ONT purification protocol. Then, the final DNA libraries were added to FLO-MIN109 flowcells (Covaris, Woburn, USA). For each group, three biological reduplicated samples were loaded on PromethION-48 platform (ONT, Oxford, UK) and the data were collected. Basecalling was performed live on the compute module using Guppy (ONT) 38 . The original data was further filtered to remove adapters, short reads and low quality reads, and then the total clean reads were obtained. Sequence depth and alignment efficiency were counted by aligning the positions of clean reads at the reference genome after mapping with Minimap2 (v2.14) software 39,40 .

Detection of SVs
SVs were identified using Sniffles (v1.0.10) 41 , an ONT long-read SV mutation detection software. Four types of structural variants including deletion (DEL), insertion (INS), inversion (INV), and duplication (DUP) were annotated. The SV results of biological reduplicated samples were combined using SURVIVOR (v1.0.5) software 42 , and the SV was re-identified for each sample based on the combined SV results. Difference analysis was performed on the re-identified SVs. If the variant genotypes identified between control group and model group were different, it was considered to be differential SVs.

Detection of CNVs
CNV detection was performed using QDNA-seq (QDNA-seqR_1. 18.0) software 43 . The genome was divided into 15 kbsized bins, and the alpha value of the circular binary segmentation of the DNA copy package less than 0.05 was accepted as the changed point. Gene annotation was performed in the region of identified CNV for each sample. For each bin, identify the copy number status under different experimental groups. If the status of the model group and the control group were inconsistent, the difference was identified as the differential CNV, and the adjacent differential bins were combined.

RNA-seq transcriptome analysis
Total RNA was isolated and the RNA quality was assessed using the Agilent Bioanalyzer 2100 System (Agilent Technologies). The ribosomal RNA was removed using the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, USA), and the RNA sequencing libraries were generated using NEB Next RUltra™ Directional RNA Library Prep Kit (NEB) for Illumina following manufacturer's recommendations, and index codes were added to attribute sequences to each sample. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kitv3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina platform and reads were generated. Clean data (clean reads) were obtained by removing reads containing adapter, reads containing poly-N, and low quality reads from raw data. These clean reads were then mapped to the mouse reference GRCm38 genome (ftp://ftp. ensembl.org/pub/release-95/fasta/mus_musculus/), and all the downstream analyses were based on the clean data with high quality.

qPCR
qPCR was performed to validate gene expressions in liver tissue. RNA was extracted from the liver by using RNAiso Plus (Trans-Gen Biotech, Beijing, China). Total RNA was reverse transcribed into cDNA with oligo (dR) primers by using TransScript All-in-One First-Strand cDNA Synthesis SuperMix (TransGen Biotech). qPCR were performed using the TransStart Top Green qPCRSu-perMix (TransGen Biotech) and subjected to CFX96 PCR system (Bio-Rad Laboratories, Hercules, CA, USA). The levels of gene were normalized to those of b-actin and gene expressions were analyzed by 2 ÀDDCt method. Each experiment was performed once with three technical triplicates per sample.

Statistical analyses
Data are presented as mean AE standard deviation (SD). Statistical analysis was performed using GraphPad Prism 8.0 Software (San Diego, CA, USA). Comparisons of statistical significant (*P < 0.05, **P < 0.01) between two groups were performed using a two-tailed unpaired Student's t-test.

Reprogrammed liver transcriptome in mice with NAFLD
The mice were induced to develop NAFLD, and the results indicated significant increase in serum levels of ALT, AST, TC, TG, and FFA; with changes in levels of MDA, SOD, and GSH indicating a state of oxidative stress in the liver of model group mice (Fig. 1A). The histopathological changes and lipid accumulation in the hepatic tissues of the model group were observed by H&E and oil red O staining (Fig. 1B).
On comparing gene expression levels by RNA-seq between mice without NAFLD (control group, C) and with NAFLD (model group, M), each sample retained more than 116 million (from 116 to 205 million) high quality clean reads after quality filtering (Supporting Information Table S1). A total of 91.44% reads (from 89.03% to 93.92%) were successfully mapped to the mouse reference genome. Coverage depth maps of the mapped reads were drawn on the reference genome by counting the position distribution of reads on each chromosome (Supporting Information Fig. S1A). The high degrees of similarity between different biological replicates of the liver samples are shown in Fig. 1C and Fig. S1B. Considering all replicates, we identified 2163 differentially expressed genes by FPKM analysis [absolute value of log 2 (fold change) > 1 and adjusted P value < 0.01] in NAFLD mice (Fig. 1D). Of the differentially expressed genes, 2100 genes were used in the study and can be annotated as differentially expressed genes (DEGs), including 1251 upregulated and 849 downregulated genes in NAFLD mice (Supporting Information Table S2). Enrichment analysis on the obtained differentially expressed genes was performed using GO and KEGG database and the top 10 most significantly enriched GO terms and the top 20 KEGG pathways are shown in Fig. 1E and F and Fig. S1C and D, respectively. Among the GO terms, "negative regulation of lipid biosynthetic process", "cellular response to food", "positive regulation of lipid metabolic process", "energy reserve metabolic process", "negative regulation of gluconeogenesis", and "positive regulation of glucose metabolic process", associated with glycolipid metabolism were enriched in the biological process. Furthermore, KEGG pathways including "fat digestion and absorption", "glycerolipid metabolism", and "PPAR signaling pathway" were enriched in the model group compared with the control group (Supporting Information Table S3). These results indicate a large proportion of perturbations in the transcriptome of a NAFLD liver, with the majority of dysregulated protein-coding genes associated with glycolipid metabolism. We selected a subset of the differentially expressed genes for validation by qPCR, and the . Green dots represent significantly downregulated genes in model group, red dots represent significantly upregulated genes in model group, and black dots represent genes with on significant difference in expression between two groups. (E) The top 10 most significant biological process terms in GO enrichment within DEGs. Abscissa represents the number of DEGs enriched in the corresponding biological process, while color of the column represents the adjusted P value by hypergeometric test. (F) The top 20 enriched KEGG pathways within DEGs. Abscissa represents the rich factor, which is the ratio of the proportion of DEGs annotated to the pathway to all genes annotated to the pathway. The ordinate is Àlog 10 (Q value) that represents the enrichment significance of DEGs in this pathway, where Q is the adjusted P value after multiple hypothesis test. expression levels of the selected genes were consistent with RNA-seq results (Fig. S1E).

Mapping 3D chromatin conformation
To get maps of 3D chromatin conformation, Hi-C libraries of liver samples were sequenced and more than 477 Gb clean data were acquired for each group in total. More than 802 million pairs of valid Hi-C data were obtained for each sample after removing PCR duplicates, and the mappability of the reads was 68.56% and 68.96% at a depth of sequencing that enabled examination of interactions with a 10 kb resolution for C and M groups, respectively (Supporting Information Table S4). The standardized interaction matrix was obtained using HiC-Pro v2.10.0, and the Pearson's correlation coefficient between the sub-libraries of each group was 0.95 (Fig. S2A). Cis-interaction refers to a read with both ends aligned to the same chromosome, whereas trans-interaction indicates the ends of the read aligned to different chromosomes. There was no significant difference in the number of cis-and transinteraction in the Hi-C libraries between control and model samples ( Fig. 2A and B, and Supporting Information Table S5). The ratio of cis-to trans-interaction revealed about three-quarters of interactions occurring within chromosomes and a quarter of the interactions between chromosomes (Fig. 2A). The interaction frequency between chromosomes can be calculated by comparing the observed read number of interactions between any two chromosomes by Hi-C and the number of expected interactions, this frequency can reveal the interaction strength between chromosomes in the genome and the relative position of the chromosomes in the nucleus. The interaction frequency heatmaps between chromosomes are shown in Fig. 2C. Differences in the interaction frequency were identified between chromosome (Chr) 2 and Chr 5, by color change in the heatmap from blue in the control group to light red in the NAFLD group, suggesting relatively closer distance and stronger interaction between the two chromosomes in NAFLD group. In contrast, the observed color changed from red to blue revealing a farther spatial distance between chromosomes and weaker interaction in the NAFLD group, such as the interaction between Chr 7 and Chr 11. The two Hi-C libraries from each group were combined in the subsequent comparative analyses between control and model group. Based on the normalized interactive signal strength of Hi-C read pairs covered between any two bin, the whole-genome interaction matrixes were drawn (Fig. 2D), and the intra-chromosomal interaction matrixes are shown in Fig. 2E and Fig. S2B. The distributing proportion of cis-and trans-interaction in the Hi-C libraries obtained from control and NAFLD mice. Liver tissues from 3 mice in each group were mixed equally, and two Hi-C libraries were established for each group. Cis-close, the cis-interactions with the linear distance between two terminals was less than 20 kb; cis-far, the cis-interaction with the linear distance between two terminals was more than 20 kb; Trans, the ends of the read aligned to different chromosomes. (B) The cis/trans-interaction ratio in each chromosome in the established Hi-C libraries. (C) The interaction frequency heatmaps between chromosomes of the livers from mice in control and model groups. Red color in heatmap suggests that the distance between the two chromosomes is relatively close and there is a stronger interaction; while blue color indicates a closer distance between the two chromosomes and a weaker interaction. Some obvious differences between control and model groups are outlined in yellow. (D) Hi-C heatmap of genome-wide interaction pattern of liver samples from control and NAFLD mice. (E) The representative Hi-C interaction matrices of Chr 2, Chr 3, and Chr 4.

Identification of alternations in A/B compartment and related DEGs in NAFLD
Besides gene interactions, the high-resolution resource generated by Hi-C in this study can explore physical environment changes of any genomic locus in the liver cells with NAFLD at multiple scales, like A/B compartments, TAD, and chromatin looping. To improve the resolution of Hi-C, two libraries of each sample were combined in subsequent analysis. The chromatin was segregated into A/B compartments by PCA to understand the differences in chromatin interaction patterns in the NAFLD liver genome. As euchromatin and heterochromatin regions are spatially separated within a chromosome territory, the chromatin can be divided into A and B compartments, corresponding to higher and lower gene expressions, respectively. The distribution of A/B compartments on each chromosome of each sample are shown in Fig. 3A and Supporting Information Fig. S3. In total, 714 A (42.93%) and 722 B (49.65%) compartments were identified in liver genome of control sample, whereas 819 A (42.64%) and 821 B (49.92%) compartments were identified in NAFLD sample (Supporting Information Table S6). Very few differences in the ratio of A/B compartments in each chromosome in NAFLD were observed (Fig. 3B). On combining the genomic and transcriptome data, the gene number and expression distributed in the A compartment were significantly higher than those in the B compartment ( Fig. 3C and D), with the results consistent with the feature of A/B compartments. The compartments were highly correlated between the normal and NAFLD conditions, with switch of only 2.34% of the genome following NAFLD inducement (Fig. 3E). The distributions of A/B switching compartments in each chromosome are shown in Fig. 3F. A total of 447 A/B compartment switching events occurred in all chromosomes, with 226 A to B switches and 221 B to A switches, with the highest number of compartment switches occurring on Chr 9 and Chr 13 (Supporting Information Fig. S4AeC and Table S7). Comparison of gene expression (RNA-seq) with the A/B switching compartments defined by Hi-C identified correlation between A/B compartment changes and the gene expressions (Fig. S4D). Screening 38 DEGs associated with the A/B compartment switching between the experimental groups identified upregulation of 94.12% DEGs in B-to-A compartment and downregulation of 33.33% DEGs in A-to-B compartment in the NAFLD group. We screened 16 upregulated DEGs related to B-to-A compartment and 4 downregulated DEGs related to A-to-B compartment switch ( Fig. S4E and Table S8). SUGCT is a mitochondrial enzyme synthesizing glutaryl-CoA from glutarate in tryptophan and lysine catabolism. SUGCT dysfunction is accompanied by lipid accumulation in the liver, and is an important contributing factor in metabolic disorders 44 . Fibroblast growth factor receptor 2 (FGFR2), a tumorigenesis-related gene is upregulated in advanced NAFLD 45 , and was also upregulated in our study, and may be associated with potential tumorigenic acquirement in NAFLD.

Identification of alternations in TADs and related DEGs in NAFLD
We identified the TADs using the TadLib software at 20 kb resolution and TADs with IR > 0 and length >5 bin were retained as the final TADs. On comparison of TADs in the control and model groups, we identified 8179 TADs in the control group and 7818 TADs in the NAFLD group, and their median TAD sizes were 180 kb (Supporting Information Table S9). Analysis of TAD architecture demonstrated TADs to be stable features of chromatin structure. Directionality index (DI) measures imbalance in the directionality of contacts from a given region of interest 46 , and is related to TAD boundaries which show a signature of strong shift from upstream to downstream DI (Fig. 4A). We identified 9390 and 9246 TAD boundaries in control and NAFLD, respectively, of which 538 are special boundaries of control and 394 are of model groups (Fig. 4B-D and Supporting Information Table S10). TADs were highly correlated between control and NAFLD groups, with only approximately 5% TAD boundaries changed upon NAFLD development ( Fig. 4C and D). TADs play an important role that in the regulation of gene expression, hence, we mapped the DEGs associated with the differential TAD boundaries and analyzed their numbers. DEGs associated with the unique TAD boundary of each group and the unchanged boundary in both groups are shown in Fig. 4E. Gene expression in NAFLD was affected by changes in TAD boundary such as the chromosome region with red shadow in Fig. 4F, as well as alternation of the domains shown as the blue shadow in Fig. 4F. We identified 42 upregulated and 30 downregulated genes (Supporting Information Fig. S5 and Table S11). These results indicate altered transcriptional profiles due to formation of new and disappearance of the original TAD boundaries associated with NAFLD development.

Identification of changes in chromatin looping and related DEGs in NAFLD
In addition to short-range chromatin interactions, long-range looping in 3D genome organization brings physical proximity between two linearly separated loci. We annotated chromatin loops through the Hi-C contact matrices using an improved HICCUP method with 10 kb resolution. The results demonstrated a decrease in the total number of loops (from 39,670 to 38,469,    (Fig. 5C). On classifying the identified chromatin loops in each group, statistical analysis indicated increase in number of loops without promoter and enhancer, and decrease in number of loops with anchor located on enhance and promoter in the NAFLD group compared to the control group ( Fig. 5D and E). We confirmed 1924 specific loops in the control group and 723 specific loops in the model group with an FDR threshold of 0.1 (Fig. 5A and Supporting Information Table S12). The relationship between loop variants and gene expression was investigated combing the RNA-seq data, and 98 mapped DEGs were regulated by variants at chromatin loops, including 48 upregulated and 50 downregulated genes (Supporting Information Fig. S6B and Table  S13). GO enrichment and KEGG analysis of these loop-impacted DEGs (Fig. S6C and D) identified hallmark genes like acyl-CoA desaturase 1 (Scd1) and elongation of very long chain fatty acids protein 6 (Elovl6) involved in de novo lipogenesis, and their expression was downregulated in NAFLD. On the other hand, levels of ATP-binding cassette transporters like ATP-binding cassette sub-family G member 5 and 8 (Abcg5 and Abcg8) that form a heterodimer and facilitate the efflux of cholesterol into bile 47 , and Fgfr2 as reported in advanced NAFLD 45 were upregulated in NAFLD in our study. These finding indicate the potential role of chromatin loops in cell metabolic events in NAFLD.

SVs in NAFLD identified by Nanopore sequencing
Nanopore sequencing was used to sequence the variations in the NAFLD mice genome. After quality control, 643.97 Gb of clean data was obtained from the control and model samples (n Z 3 for each group), with an average clean data volume of 107.33 Gb. The N50 value was 9.89 kb, and the average length of the reads was more than 6.95 kb (Supporting Information Table S14), with the reads length distribution of each sample given in Supporting Information Fig. S7A. The mapping ratio in each sample was over 99.04% in each sample according to the reference genome (Supporting Information Table S15). The distribution of sequencing depth and read coverage showed acceptable randomness of sequencing ( Fig. 6A and B, Supporting Information Figs. S7B and S8). The Sniffles software was used to detect SVs in Nanopore sequencing and the SUR-VIVOR software was used to identify differential SVs between the NAFLD and control groups; which led to the discovery of 10,917 differential SVs, including 6568 DEL, 3785 INS, 405 TRA, 141 INV, and 18 DUP (Fig. 6C and Supporting Information Table S16). A majority of differential variations identified were DEL and INS, with a higher distribution number in Chr 10 to Chr 17 ( Fig. 6D and E). All annotated transcribed features related to differential SVs were retrieved, to evaluate the functional impact of the differential SVs. Between the two  Information Table S17). Based on gene function, we found 14 SV-impacted DEGs involved in lipid transport and metabolism (Acaca, Acsl5, Cbr3, Ces1f, Ces1g, Cryl1, Dgka, Dgkb, Dgkk, Idi1, Nlgn1, Nlgn3, Sugct, and Hao2), and 8 SVimpacted DEGs associated with energy production and conversion (Bckdhb, Me1, Hao2, Mt-co3, Mt-nd4, Nnt, Slc1a1, and Slc1a2). Some of the SV-impacted DEGs are associated with NAFLD, however no association with NAFLD has been reported for some SV-impacted DEGs and may present as new molecular candidates in NAFLD. These results implicate SVs in the development of NAFLD through impact on gene regulation and chromatin organization.

CNVs in NAFLD identified by Nanopore sequencing
The deletion or duplication of DNA fragments causes copy number variations (CNVs), which are an important source of both normal and pathogenic variations in the genome. In this study, approximately 150 CNVs were found in each sample including loss, gain, and amplification (Supporting Information Fig. S9). Differential CNVs between the control and NAFLD groups were identified by combining with the results of RNA-seq, and 22 DEGs with differential CNVs including 12 upregulated and 10 downregulated genes were associated with the development of NAFLD (Supporting Information Table S18). CNV variants in Cntnap2, Mrgpra2b, Hist1h4n, Hist1h2br, Hist1h2ao, Pde4b, Gsta1, Gm3776, Kcnh1, and Camk1d were simultaneously identified on comparing SV-impacted DEGs. Camk1d may have a potential role in liver gluconeogenesis 48 , and a difference in CNV in the Camk1d locus was observed between the two groups, with a downregulation of gene expression in NAFLD group. A CNV difference was found in the locus of Gsta1, with increased expression in the NAFLD group, as previously reported and is associated with oxidative stress in NAFLD 49 . It is thus evident that mechanisms of gene regulation in NAFLD are associated with a difference in chromatin microenvironment.

Pathogenic gene candidates due to reorganization of chromatin structure in NAFLD
The primary goal of the study was to identify variations in genomic and spatial features to reveal the underlying pathogenesis of NAFLD. We comprehensively assessed the Hi-C and whole genome sequencing results and extracted the structure variations data by Nanopore sequencing, the 100 highest chromatin interactions in Hi-C, and the significant changes in structure variation and chromatin interaction between the control and NAFLD group (Fig. 7A). The occurrence of NAFLD is the results of these factors. Hence, we integrated Hi-C, Nanopore data with DEGs in NAFLD, and all DEGs that impacted by 3D disorganization or rearrangement were screened by comprehensive analysis of RNA-seq, Hi-C, and Nanopore sequencing results (Supporting Information Table S19). There were totally 555 The results of contact matrix, triangular, domain, and loop reflect the changes in spatial organization generated by Hi-C. SV and CNV were obtained from Nanopore sequencing. The abundance of gene expressions was generated by RNA-seq.
DEGs related to the SVs in A/B compartment, TAD, chromatin loop, and/or CNV. Moreover, we screened DEGs that were related to the pathogenesis of NAFLD in gene ontology annotation and existed variation in genomic structure as the potential candidate genes, which may contribute to NAFLD development (Table 1). Variations in only CNV and SV in the locus of the upregulated Gsta1 gene, without change in reorganized spatial structure, were discovered (Fig. 7B). The expression of Fgfr2 was upregulated in the NAFLD group, with no variations in TAD, loop, and CNVs except for A/B compartment, the switching in which may be associated with Fgfr2 upregulation (Supporting Information Fig. S10A). Plxnb1, a suppressor of liver tumor formation 50 , was downregulated in NAFLD group, with variations observed in its spatial organization ( Table 1). The deregulation of tumor-related genes in NAFLD may increase the potential of tumorigenic acquirement. Downregulation of Elovl6 was associated with the reorganization of 3D structure in TAD and loop (Fig. S10B). Camk1d expression was downregulated and comprehensively affected by alternations in CNV and SV, chromatin loop, domains, and interaction matrix (Fig. 7C). These results indicate the functional relevance of each level of the 3D genome hierarchy, and the establishment of 3D genome provides multiple regulatory layers to alter gene expression in NAFLD development.

Discussion
NAFLD is the most common liver disease worldwide and can progress to inflammatory hepatitis, cirrhosis, and hepatocellular carcinoma. Due to an imbalance between hepatic fatty acid uptake and removal, the hallmarks of NAFLD are metabolic disorders and lipotoxic hepatocyte damage. It is imperative to alleviate the conditions of NAFLD; however, the pathogenesis of NAFLD is complex and needs to be elucidated further. In the present study, we investigated the 3D genome of NAFLD and performed integrated analyses using Hi-C, Nanopore sequencing, and RNA-seq. Mice were fed a HFD to induce NAFLD, and physiologic measures of abnormally elevated lipid levels, serum transaminases, and oxidative stress state were observed. Additionally, inflammatory injury and visible accumulation of lipid droplets in the liver of NAFLD mice was observed. These metabolic and physiological changes reflect the occurrence of hepatic damage in NAFLD mice.
The hepatic injury in NAFLD is accompanied by alternations in gene expression and an underlying remodeling of the liver transcriptome. Evidence suggests functional relevance of the different 3D genome hierarchies, with the 3D genome organization providing multiple additional regulatory layers to control gene expression 51 . The 3D genome organization contains spatial compartmentalization, organization into TADs, and chromatin looping, which modulate several processes including transcription, DNA repair, and replication. Due to the essential role of 3D genome organization in physiological conditions, their derangements may lead to potential pathogenic consequences. We generated and characterized genome-wide chromosome conformation Hi-C interaction matrices from the liver of normal and NAFLD mice to link the spatial regulatory elements and gene expressions. We identified reorganization in the interactions within chromosomes and between different chromosomes in the livers of NAFLD mice, and the relationship between the reorganizations and abnormal gene expressions were analyzed hierarchically.
Chromosomes are partitioned into either the active and open chromatin named A compartment or the inactive and closed chromatin called B compartment. The characteristics of A compartment include high GC-content and gene-rich regions, with generally high transcription levels. Conversely, B compartment is Table 1 The candidate genes that might contribute to NAFLD screened from the comprehensive analysis of the results of RNA-seq, Hi-C, and ONT in this study. TADs are the basic units of chromosome and are self-contact domains conserved among tissues and species 53 . Smaller domains within TADs called contact domains can be observed by using high resolution Hi-C, and are detected as peaks on the interaction matrices and correspond to more frequent contact regions than flanking loci 29 . TAD organization is important for genomic integrity and disease, and the disruption of TADs can cause abnormal contacts between enhancers and genes, thereby activating or inhibiting disease-relevant gene 54 . TAD disruptionmediated oncogene activation is reported to be a key mechanism in several cancer types 55 . In NAFLD, on drawing a map of the TADs, most TADs were conservative, however approximately 5% of TAD boundaries were disrupted upon NAFLD inducement. DEGs mapped to the changed TADs were screened, and downregulation of Elovl6, a hallmark gene involved in de novo lipogenesis was observed in relation to TAD variation in NAFLD. Plxnb1, a suppressor of liver tumor formation, was a TADimpacted downregulated DEG in NAFLD. The disruptions in TADs and affected DEGs provide a new insight into understanding the pathogenesis of NAFLD.
Chromatin loop is an intrinsic property of the chromatin fiber and a fundamental mechanism for building the 3D genome hierarchy 56 . Chromatin looping plays an important role in gene transcription and RNA splicing by bringing distant genomic regions between the gene promoter and other regulatory elements, physically adjacent 57 . Therefore, chromatin looping provides a mechanism to magnify the long-range transcriptional and epigenetic effects of noncoding genetic variants. We identified differences in chromatin loops between the two groups and 98 DEGs associated with changes in looping. The TAD-impacted genes, Elovl6 and Plxnb1, were affected by variations in chromatin looping as well. It may be possible that loop extrusion is a mechanistic aspect giving rise to TADs 40 .
The mechanisms of derangements in 3D genome may involve preternatural mutations of 3D genome organizers or genetic changes including SVs and CNVs in the genome 58 . It is not determined yet if 3D genome disruption is the cause or the consequence of genetic variations, however abnormal 3D genome organization and genetic variations together cause dysregulation of genome function and hence, the phenotypes of different diseases 17 . Gene dysregulation caused by TAD disruption is a novel mechanism of SV-induced congenital developmental diseases 59 . Genomic rearrangements associated with human limb malformations altered the TAD organization in the WNT6/IHH/EPHA4/ PAX3 locus, as reported by Lupiáñez et al. 60 . The application of next-generation sequencing (NGS) has obviously improved the throughput of data generation, which can make genome sequencing much faster and cheaper 61 . However, there are several drawbacks and pitfalls of NGS, most notably their relatively short reads. There are commonly a large number of repeated sequences in genomes that are longer than the NGS reads, which may cause misassemblies and gaps. Large SVs are more challenging to detect and characterize for NGS 62,63 , that is an important issue due to the relevance between SVs and many diseases. In addition, NGS methods rely on PCR causes difficulties with regions of extreme GC% 63 . Nanopore sequencing distinguishes itself from NGS, in that it directly detects the input molecule without DNA synthesis or amplification, so there is no apparent limit to the length of DNA that can be sequenced 63 . Nanopore sequencing can not only directly observe the canonical four nucleobases, but also base modifications such as methylations 64 (the methylation data of this study will be uncovered in our following research). The main drawback of long-read sequencing is the relatively high error rate compared to NGS. Nanopore sequencing improved data quality by sequencing both template and complementary strands to obtain a final more accurate consensus read. Considering the above factors, Nanopore sequencing, a long-read method, was employed in the current study. We identified all SVs and CNVs in the genome, and retrieved the DEGs related to differential SVs and CNVs in NAFLD. We identified 381 SV-impacted and 23 CNV-impacted DEGs associated with NAFLD development, and dysregulation of Camk1d, Pde4b, and Gsta1 genes was due to presence of SVs and CNVs in their locus. Additionally, presence of genomic rearrangements (SV or CNV) around the locus of DEGs including Plce1, Plxnb1, Fgfr2, Grm8, Kcnma1, and Sugct, was accompanied by disruptions of different 3D genome hierarchies. This indicates the presence of shared mechanisms that underlie variations in 3D chromatin conformation, genome, and transcription regulation, and may influence hepatocyte phenotypes in NAFLD.
There are some limitations to the current work. The regulatory role of 3D genome in the nuclear processes needs to be further elucidated in NAFLD. Genome editing studies to determine the functional consequences of perturbing the various 3D organization features are required in the future to clarify the mechanisms of 3D genome disruption causing NAFLD. Additionally, deciphering the combination of epigenetic mechanisms mediated by histone modification, DNA methylation, and noncoding RNAs is required for precise determination of genetic rearrangements and spatial organization disruption, and may likely uncover new mechanisms of NAFLD and identify novel therapeutic targets.

Conclusions
By combining transcriptome, long-read sequencing, and 3D genome analyses, we revealed occurrence of multiple levels of alternations including SV, CNV, spatial genome reorganization, and their influences on gene expression in NAFLD. Our results reveal gene regulation through genetic rearrangements and spatial organization disruption and identify candidate genes associated with the pathogenesis of NAFLD. These findings extend our knowledge of the spatial disorganization and genetic rearrangements of the liver genome of NAFLD and may have implications on clinical treatment and drug development for NAFLD.