Three-dimensional (3D) organization of mammalian genomes is now considered as an additional level of the regulation of gene expression. Molecular approaches, such as Hi-C, provided a view of genome spatial organization, and the high-resolution Hi-C contact maps in the mammalian revealed the plaid pattern of chromatin interactions^[1]. These interactions displayed the segregation of the genome in mega-base scale, namely active A, and inactive B compartments^[1]. Topologically associating domains (TADs) are structures in sub mega-base, corresponding to sequences that interact preferentially with themselves rather than with other regions of the genome^[2-3]. TADs are separated by boundaries that are enriched in sharp CTCF binding peaks^[4]. By promoting or preventing interactions between promoters and enhancers, TADs may contribute to transcriptome regulation^[5]. Chromatin loops are focal enrichments in spatial contact of point-to-point on the heatmap, and previous studies suggest that loops can bring gene promoters to distant regulatory elements within the same TAD^[5-6].

The juxtaposition of enhancer and promoter creates a transcription dependent phase-separated condensate with mediator and Pol Ⅱ, which is important for transcriptional regulations^[7-8]. The cis-elements that mediate the dysregulations of oncogenes and tumor suppressor genes are potential neoplastic cell transformation mechanisms^[9]. We recently performed several high-throughput sequencing assays, including in situ Hi-C between two representative liver cell lines: non tumorigenic human liver cell line L02 and human hepatocellular carcinoma cell line PLC/PRF/5. Hepatocellular carcinoma (HCC) is the most common type of hepatic carcinoma, accounting for approximately 80%–90% of primary liver cancer^[10-11]. In this study, we demonstrated that compared with L02, the HCC cell line PLC/PRF/5 (or PLC) has significant differences in various levels of 3D genome, including compartment, TAD and looping. By using ChIP-seq and RNA-seq, we also proved the correlation of 3D genome change with histone modification and gene transcription in the comparison between these two cell lines. Therefore, we demonstrated the molecular basis of genome spatial contact in HCC neoplastic cell transformation.

1 Materials and Methods 1.1 Cell culture

PLC/PRF/5 was obtained from the Type Culture Collection of the Chinese Academy of Sciences (Shanghai) and was authenticated by short tandem repeat PCR to be 100% identical to PLC/PRF/5 cells from ATCC. PLC/PRF/5 was treated with Mycoplasma Removal Agent (MP Biochemicals, 3050044) to remove the mycoplasma. L02 was obtained from Liver Cancer Institute of Zhongshan Hospital (Fudan University). Both of the cell lines were cultured in the high glucose DMEM (Gibco, 11965084) with 10% FBS (BI, 04-001-1ACS) and 1% Penicillin-Streptomycin (Gibco, 15140122), in a humidified incubator at 37 ℃ containing 5% CO₂. Cells were harvested by 0.25% Trypsin-EDTA (Gibco, 25200072) for further experiments.

1.2 Generation of in situ Hi-C libraries

For Hi-C assay, 1×10⁶ cells were harvested and fixed by formaldehyde. As described, Hi-C was performed by using the Dpn Ⅱ restriction endonuclease^[12]. Libraries were prepared by using the VAHTS^® Universal DNA Library Prep Kit for Illumina V3 (Vazyme, ND607) and purified through double sized selection process by using VAHTS DNA Clean Beads (Vazyme, N411).

1.3 ChIP-seq and RNA-seq

Chromatin immunoprecipitation (ChIP) of CTCF, SMC3 and H3K27ac were performed in L02 and PLC/PRF/5 cell lines as described^[13]. The following antibodies were used for ChIP: rabbit anti-SMC3 antibody (Abcam, ab9263), rabbit anti-CTCF antibody (Millipore, 07-729), and rabbit anti-H3K27ac antibody (Active Motif, 39133). ChIP-seq libraries were prepared in the same way as in situ Hi-C. RNA extraction was performed by using TRIzol (Invitrogen, 15596026). Reverse transcription for cDNA and transcriptome libraries were performed by Shenzhen HaploX Biotechnology Co., LTD. Each of the RNA-seq, ChIP-seq and Hi-C libraries were sequenced by using paired-end sequencing in an IlluminaHiseq X10 System.

1.4 Hi-C data analysis

Each of the in situ Hi-C reads was mapped to the human genome hg19 by using HiC-pro^[14]. Reads mapped in close proximity to the DpnⅡ restriction sites (5 bp), read pairs from invalid ligation products such as dangling end and self-circle ligation, and possible PCR duplicates were filtered out. The contact maps were iteratively corrected by using cooler, and normalized by the total number of sequencing reads^[15]. The contact matrices were subjected to the principal components analysis (PCA) in 25 kb resolution, and thus to define the compartment A/B^[1]. TADs were annotated by using HOMER and scored by the inclusion ratio (IR)^[16]. We took the IR value for 1.5 as the minimum inclusion score to call TAD. The annotation of overlapping TADs was calculated by BEDTools^[17]. Loops were also scored by HOMER, through the spatial interaction density which was normalized to the read depth^[16]. And the loops were called by 1.5 of fold threshold for local loop enrichment and 2 of fold threshold over avg. interactions of similar distance. The Hi-C heatmaps were visualized by Juicebox^[18].

1.5 ChIP-seq and RNA-seq data analysis

ChIP-seq data was aligned by using bowtie2^[19] to against hg19 version of the human genome, and multiple mapped and PCR duplicate reads have been sorted and removed by SAMtools^[20]. The bigwig files were converted by usingdeepTools and normalized with RPKM^[21]. Peak calling was performed with MACS2^[22]. The matrixes and heatmaps were generated by using deepTools ^[21]. RNA-seq data was remove low quality bases and adapter sequences from sequenced reads by using Trim_galore software. The filtered reads were mapped against hg19 human genome with bowtie2^[19]. Transcriptome quantification was performed by salmon software^[23]. Differentially expressed genes were analyzed by using the edgeR package based on R software and screened with P value < 0.05 and |fold change| ≥2 to obtain a list of differentially expressed genes between groups ^[24]. GO and KEGG enrichment analysis of differentially expressed genes are performed using the clusterProfiler package^[25].

2 Results and discussion 2.1 RNA-seq and Chip-seq profiling shows different patterns between L02 and PLC/PRF/5

To investigate the transcriptome variations across the cell lines, we performed the RNA-seq of the L02 and PLC/PRF/5. The KEGG analysis shows that the top 500 up-regulated genes in PLC/PRF/5 is significantly correlatedwith various lipid metabolismincluding cholesterolmetabolism (Fig. 1A). According to previous studies, high cholesterol and low ceramide levels can protect the hepatoma cell from oxidative stress and apoptosis, and thus to enhance cell proliferation^[26].

Whereas various cell signal transduction pathways and drug metabolism were enriched by the top 500 down-regulatedgenes (Fig. 1B). NOD-like receptor NLRP12 has been reported to suppress hepatocellular carcinoma via the regulation of JNK signaling^[27], and the silencing of AGE-RAGE signal pathway reduced collagen deposition and the tumorigenesis of HCC^[28]. In summary, the transcriptome alterations between PLC/PRF/5 and L02 suggest that it's a suitable cell model to study the underlyingmolecular basis in HCC.

To better understand the potential molecular mechanisms of the spatial genomic organization, we preformed Chip-seq of CTCF, SMC3 and H3K27ac in two cell lines. The only somatic insulator protein CTCF can block communication between enhancers and upstream promoters so as to regulate imprinted expression. SMC3 is a key component of the cohesin complex that holds together the non-adjacent regions of the chromatin^[29]. Notably, both of the CTCF and SMC3 binding levels around TSS significantly decrease in PLC/PRF/5 compared with L02, which may refer to more opening chromatin statues around promoters (Fig. 1C).

2.2 A/B compartment switching correlation with transcriptome alterations

We next performed in situ Hi-C in L02 and PLC/PRF/5 cell lines and calculated the PC1 values of the whole genome in 25 kb resolution. PC1 value is positively correlated with chromatin activity, and thus divides the whole genome into two parts: the region with positive PC1 value is the "active" compartment A, which usually represents the euchromatin, and the negative score region is the "inactive" compartment B, mostly representing the heterochromatin^[1]. In our experiments, L02's whole genome takes 50.80% of compartment A and 49.20% of compartment B, while the PLC/PRF/5's whole genome takes 50.07% of compartment A and 49.93% of compartment B. The compartment A/B ratio is similar between the two cell lines, but notably, within the same regions significant compartment switching (~13%) have occurred (Fig. 2A).

The histograms of PC1 values around the TSS indicate more opening chromatin statues around promoters in PLC/PRF/5, which are consistent with the SMC3 and CTCF binding levels around TSS (Fig. 1C & 2B). Among the studies of individual loci, we found that genes that change from compartment A to B tend to show reduced expression, whereas genes that change from B to A tend to show a higher expression level (Fig. 2C) which is consistent with previous studies^[1]. Next we picked two typical genes in the top 500 up/down list and checked their compartment statues: GLI3 is a transcriptional factor which is only expressed in L02, while APOB only expressed in PLC/PRF/5, which was the main apolipoprotein of chylomicrons and low density lipoproteins (LDL) (Fig. 2D). GLI3 locates in the A of the L02 cell line, and the same locus had switched to B in PLC/PRF/5. APOB locates in the B in L02 cell line, and this locus had switched to B in PLC/PRF/5 which is consist with the gene expression. (Fig. 2E).

2.3 TADs and loops partially change in PLC/PRF/5

TADs and loops are smaller scale structures of 3D genome that is relative to compartments^[2]. We called 3 393 TADs in L02 and 2429 TADs in PLC/PRF/5 by using HOMER package, and 1920 TADs across the cell lines are overlapped (Fig. 3A). In a typical genomic region, we observed a strong TAD in L02, which has vanished in PLC/PRF/5 with significant reduction of CTCF and SMC3 at its boundaries (Fig. 3B). In addition, TAD boundary shifting can also be observed in PLC/PRF/5, which means that TAD boundaries change their locations on the genome (Fig. 3C). We next calculated the inclusion ratio (IR) of TADs: high IR indicates a dense triangle that makes few interactions to adjacent regions. The results indicate that the TADs in PLC/PRF/5 are less robust and the boundaries are weaker compared with L02 (Fig. 3D & 3E). Previous studies suggest that the major function of TADs in gene expression regulation is to restrict promoter-enhancer interactions as microenvironments^[30-31]. The alteration of TADs pattern in HCC cell line may lead to aberrantly tumorigenic gene expressions.

To further investigate the spatial interactions of the liver cell lines, we called the loops of the genome. Notably, certain pairs of loci has a greater number of Hi-C reads that is expected by chance, which is referred to as a 'significant interaction'. The histogram showing the distribution of interaction lengths indicates that the significant interaction levels are almost same between L02 and PLC/PRF/5 (Fig. 4A). The results may indicate that the significant interactions are robust, and thus remain stable across the cell lines.

2.4 Spatial interaction change may lead to differential gene expression

By deeply focusing on significant interactions that contribute to gene expression, it can reveal the physiological implications of genomic structure. Despite of the similar number of significant interactions (L02: 2902, PLC/PRF/5: 2775), the loop patterns are different between the two cell lines, and thus may lead to different transcriptome (Fig. 4B & 4C). In the example of OAT, a mitochondrial ornithine aminotransferase, expressed only in PLC/PRF/5: the ChIP-seq data reveals that the promoter region of OAT gains high levels of H3K27ac binding, which marks the active gene transcription, the CTCF and SMC3 binding levels at OAT promoter and the upstream enhancer (the broad region with massively H3K27ac binding) increased, and correspondingly, the interaction between OAT promoter and enhancer significantly heightened (Fig. 4D & 4E). Notably, a new TAD emerged around the region, suggesting that TADs and loops may have synergistic regulation to genes expression (Fig. 4E). In mammalian, the enhancer- promoter looping spatially tightens the cis elements together and initiates gene transcription^[32]. According to previous studies, OAT is a β-catenin target gene that is overexpressed in HCC^[33], and the abnormal activation of Wnt/β-catenin signaling pathway leads to highly expressed target genes including OAT, which can promote cell cycle progression and liver neoplasms^[34-35]. Inhibition of OAT expression in HCC with small chemical molecules would suppress the HCC growth in mice^[36]. Furthermore, we picked others in the top 500 up/down regulated genes to show the annotated contacts around the TSS, and there are significant loop alterations that are correlated with the gene expressions (Fig. 4C).

3 Conclusion

In recent years, the progress in next-generation sequencing has substantialy advacned the understanding of the genes, transcriptomes and epigenomes in human cancer. Here, we provide joint profiling of transcriptomes, epigenomes, and 3D chromosomal interactions in HCC cell line PLC/PRF/5 in comparison with normal human liver cell line L02, thus exploring the underlying epigenetic basis of neoplastic transformation of hepatocytes. In the transcriptional level, the KEGG analysis reveals a typical gene expression pattern of HCC in PLC/PRF/5 by the significantly overactivatedcholesterol metabolism and the various down-regulated cell signal pathway, including the NOD-like receptors and the AGE-RAGE signal pathway. ChIP-seq analysis shows that the CTCF and SMC3 binding levels around the promoters of whole genome decreased in the HCC cell line (Fig. 1C). In addition, the genome regions around TSS became more opening due to the reduction of cohesin and CTCF, consisting with the result of PCA analysis (Fig. 2B). CTCF is the genome insulator protein which binds to DNA specific sequence sites directionally^[37]. CTCF and cohesin together have two apparently opposite functions: to help enhancers find their cognate promoters inside the TAD and CTCF loops and to restrict enhancer-promoter interactions between sequences located outside CTCF loops and TAD boundaries^[6-38]. Therefore, the reduction of CTCF and cohesin binding around TSS in HCC cell lines may lead to both gene activationand deactivation, depending on the location of the enhancer.

Principal component analysis reveals about 13% of compartments A and B switching in the HCC cell line, and distinct subsets of genes have concordant alterations of A/B compartments status and expression levels. Two loci (APOB and GLI3) show a similar trend of gene expression. These results indicate that changes in compartment status may influence the accessibility of genomic regions to regulatory proteins such as transcription factors, which may be important to certain subsets of genes.

In PLC/PRF/5, the inclusion ratio calculation shows a decreased number of TADs and less robust TAD boundaries (Fig. 3A & 3D & 3E). Notably, one major function of TADs is to restrict promoter-enhancer interactions^[30-39], thus the weakened TAD structures in PLC/PRF/5 may partially lead to the dysregulation of gene expression. TAD boundaries are defined by point-to-point interactions between two sequences bound by CTCF. These CTCF loops are relatively stable and stay conserved in a subpopulation of cells, and they displayed as strong punctate signals in Hi-C heatmaps^[12]. In our study, the "significant interaction" shows similar features with CTCF loops and the "number vs. distance" patterns are very close between these two cell lines (Fig. 4A). In contrast, other types of loops are not flanked by CTCF, and characterized by the presence of specific histone modifications^[12]. Notably, only a minority of CTCF binding peak coincides with active genes and enhancers, possibly indicating that the promoter-enhancer loops are more dynamic than CTCF loops^[32]. According to our observation, most of the altered enhancer-promoter contacts in the HCC cell line are non-CTCF loops (Fig. 4C), which consist with previous researches. But we also demonstrate a typical enhancer-promoter looping flanked by CTCF on the loci of OAT, which only exist in PLC/PRF/5 (Fig. 4E). In previous studies, deletion or inversion of the CFCT binding motif within the Sox2 super enhancer had negligible effects on the enhancer-promoter contact frequency^[40]. This result may suggest that the loop structure of enhancer-promoter contact is maintained in a CTCF-independent manner, even flanked by CTCF binding peaks.

Notably, we only tested one HCC cell line in this study, and PLC/PRF/5 contains the hepatitis B virus (HBV). Physical integration of HBV DNA into the genome of this cell line has been demonstrated by Southern blot analysis, and it is estimated that at least four complete and two partial copies of the viral genome are incorporated^[41]. The HBV genome is too small (~3 200 bp) to reach the resolution of Hi-C sequencing. Furthermore, the genome sequences were mapped to the hg19, which means that the exogenous HBV DNA was filtered by the software, and therefore HBV DNA rarely influence the results. However, the integration of HBV DNA may also lead to DNA fragment deletion and chromosomal rearrangement, and further change the 3D status of the genome. Up to now, there is no relevant research showing that the integration of HBV DNA in PLC/PRF/5 affects the 3D genome, it will be further verified in our subsequent project.

In summary, we show that the HCC cell line PLC/PRF/5 displays significant compartment A/B switching, weaker TAD structures, loop pattern alterations and more opening promoter genomic regions in comparison with the normal hepatocyte cell line L02. These results may help us to better understand the epigenetics of HCC tumorigenesis, together with transcriptome and genomic sequencing.