Topological reorganization and functional alteration of distinct genomic components in gallbladder cancer

Abstract

Altered three-dimensional architecture of chromatin influences various genomic regulators and subsequent gene expression in human cancer. However, knowledge of the topological rearrangement of genomic hierarchical layers in cancer is largely limited. Here, by taking advantage of in situ Hi-C, RNA-sequencing, and chromatin immunoprecipitation sequencing (ChIP-seq), we investigated structural reorganization and functional changes in chromosomal compartments, topologically associated domains (TADs), and CCCTC binding factor (CTCF)-mediated loops in gallbladder cancer (GBC) tissues and cell lines. We observed that the chromosomal compartment A/B switch was correlated with CTCF binding levels and gene expression changes. Increased inter-TAD interactions with weaker TAD boundaries were identified in cancer cell lines relative to normal controls. Furthermore, the chromatin short loops and cancer unique loops associated with chromatin remodeling and epithelial–mesenchymal transition activation were enriched in cancer compared with their control counterparts. Cancer-specific enhancer–promoter loops, which contain multiple transcription factor binding motifs, acted as a central element to regulate aberrant gene expression. Depletion of individual enhancers in each loop anchor that connects with promoters led to the inhibition of their corresponding gene expressions. Collectively, our data offer the landscape of hierarchical layers of cancer genome and functional alterations that contribute to the development of GBC.

Keywords

3D genome / cancer / TADs / loop / gene regulation

1 Introduction

Genomic alteration emerges as a fundamental driver of tumorigenesis. Recent compelling pieces of evidence offer new insights into chromatin structural reorganization and functional changes that regulate gene expression and activation during cancer development [1]. The three-dimensional (3D) analysis of chromatin architecture defines a functional DNA structure that orchestrates chromatin loops in which promoters (P) and enhancers (E) coordinately regulate the expressions of a variety of genes. Chromatin loops develop when CCCTC binding factor (CTCF) and cohesin located at distant loci on a linear chromosome interact to form a topological close-up structure [2]. The disruption of CTCF binding sites or their orientation in loop anchors can disrupt enhancer–promoter (E-P) interactions and impair gene expressions, which ultimately gives rise to an abnormal phenotype [3,4]. Individual chromatin CTCF loops spatially compact together to form submegabase regions, known as topologically associated domains (TADs), which are highly self-interacting and reflect genomic signatures [5]. The DNA genome is glossily designated into large-scale compartments A and B, in which compartment A involves an accessible and transcriptionally active genome, whereas compartment B is compact and relatively silent. Reorganization and conversion between these compartments frequently occur in cancer [6–9]; however, the molecular mechanisms underpinning topological reorganization and carcinogenesis remain to be established.

Gallbladder cancer (GBC) is the most common malignancy of the biliary tract system, with a median survival of less than one year [10]. The poor prognosis of this disease, which lacks notable symptoms at the early stages, is largely attributed to the delayed diagnosis [11]. A total of 25% of patients are considered acceptable for curative surgery, and less than 16% can survive for more than five years [12]. Thus, improvements in early cancer detection and prevention are eminent. GBC prevention requires the decreased levels of cancer risk factor, including cholecystitis, cholelithiasis, and gallbladder polyps [12,13]. In addition, a large body of research evidence has revealed a variety of key molecules that govern cancer development. Our previous work primarily focused on the molecular and genetic signatures of GBC and demonstrated that genomic dysfunctions promoted the malignant transformation of GBC [14–18]. However, topological changes in genome and chromatin organization in GBC are poorly understood.

In this study, we employed in situ Hi-C, RNA-seq, and ChIP-seq to systemically investigate the hierarchical layers of cancer genome reorganization and subsequent functional changes in GBC tissues and cell lines. We observed that 3D reconfiguration and functional changes in compartments, TADs, and CTCF-mediated loops that harbor genes (e.g., BMP4, KRT19, and FOXA2) driven by E-P interaction coordinately contributed to the genotypic and phenotypic changes in GBC cells, which promoted cancer cell growth. This study sheds light on the newly discovered signature of chromatin reorganization and activity that regulates the expression of cancer-related genes and drives the development of GBC.

2 Materials and methods

2.1 Patients and clinical specimens

Gallbladder cancerous and paired paracancerous tissues were obtained from patients who underwent cholecystectomy without neoadjuvant chemoradiotherapy or endocrinotherapy at the Department of General Surgery, Xinhua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, China. All the patients were subjected to accept the written informed consent before enrollment. The study was approved by the ethics committee of Xinhua hospital (No. XHEC-D-2021-071) and the study was performed in accordance with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. Informed consent was obtained from all patients for being included in the study.

2.2 Cell culture and reagents

The GBC-SD cell line was purchased from the cell bank of the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). The NOZ cell line was obtained from the Health Science Research Bank (Osaka, Japan). Both cell lines were maintained in high-glucose Dulbecco’s Modified Eagle Medium supplemented with 10% fetal bovine serum, penicillin G (100 U/mL), and streptomycin (100 g/mL). The cells were maintained as monolayer cultures at 37 °C in humidified air with 5% CO₂ and 95% air. Human biliary epithelial cells (HBECs) were purchased from ScienCell (#5100, USA) and characterized by immunostaining (CK18).

2.3 In situ Hi-C library construction

An in situ Hi-C library was constructed as described by Rao et al. with minor modifications [19]. Briefly, approximately two million cells from gallbladder cancerous and paracancerous tissues were prepared as a single-cell suspension and readied for crosslinking. The cell samples were incubated in 1% formaldehyde for 10 min at room temperature (RT) and quenched by the addition of proper volumes of glycine (final concentration: 0.2 mol/L). Then, the samples were lysed with an ice-cold Hi-C lysis buffer containing 10 mmol/L Tris-HCl (pH 8.0), 10 mmol/L NaCl, 0.2% Igepal-CA630, and 1× protease inhibitors. After centrifugation, the pelleted nuclei were resuspended in 0.5% sodium dodecyl sulphate (SDS) and incubated at 62 °C for 10 min, followed by the addition of 10% Triton X-100 to quench SDS. A restriction enzyme of 100 U MboI (R0147, NEB, USA) was used to digest the chromatin overnight at 37 °C, followed by incubation at 62 °C for 20 min to inactivate MboI. The DNA ends were filled with Klenow enzyme (M0210, NEB, USA) and simultaneously labeled with biotin-14-Datp (#19524016, Thermo, USA). Chromatin was ligated again using 2000 U T4 ligase (M0202, NEB, USA) for 4 h at RT. Then, the samples were digested with proteinase K (20 mg/mL) for 30 min at 55 °C, followed by overnight incubation at 68 °C to reverse the crosslinks. DNA was precipitated using ethanol and sodium acetate and fragmented from 300 bp to 500 bp with sonication. Biotin-labeled DNA was precipitated using streptavidin-conjugated beads (#65602, Life Technologies, USA), and the beads were washed thrice with 1× Tween wash buffer containing 5 mmol/L Tris-HCl (pH 7.5), 0.5 mmol/L ethylenediaminetetraacetic acid, 1 mol/L NaCl, and 0.05% Tween-20 at 55 °C. A DNA library was established through DNA end repair, dA-Tailing, and adaptor ligation using the NEB DNA library prep kit (E7370, NEB, USA). After washing, DNA samples were amplified via 10 cycles of polymerase chain reaction (PCR). The 300–500 bp DNA fragments were selected with Ampure XP beads (A63880, Beckman Coulter, USA) and sequenced using the Illumine NovaSeq platform.

2.4 ChIP-seq library construction

ChIP was performed using the Enzymatic Chromatin IP Kit, in accordance with the manufacturer’s instructions (#9003, CST, USA). Briefly, cells were cross-linked with 1% formaldehyde for 10 min at RT, and chromatin was fragmented by digestion with micrococcal nuclease. The lysate was then immunoprecipitated using antibodies against CTCF (#3418, CST, USA), H3K27ac (#8173, CST, USA), and H3K4me3 (#9751, CST, USA). Next, DNA was purified using DNA purification spin columns. Finally, a DNA library was constructed using the NEB DNA library prep kit and sequenced by the Illumina NovaSeq platform.

2.5 RNA-seq library construction

Total RNA was extracted from the cell lines and tissues using TRIzol reagent (#15596026, Invitrogen, USA). The resulting libraries were generated using the NEB RNA Library Prep Kit for Illumina (E7770, NEB, USA) following the manufacturer’s instructions, and barcodes were added to the attribute sequences of each sample. RNA-seq was performed with two biological replicates for GBC cell lines, normal epithelial cells, and GBC-1 tissues. No biological replicate was performed for GBC-2 tissues. The pooled libraries were sequenced on the Illumina NovaSeq platform, and paired-end reads were obtained.

2.6 Contact heatmap analysis

Adapters of raw Hi-C data were trimmed by cutadapt (version 3.2). HiC-Pro [20] was used to map the previously generated clean data to the human reference genome (hg19). Contact matrices were normalized by iterative correction and eigenvector decomposition (ICE-normalized) [21] and further normalized using 100 million contacts per sample. The Matplotlib plotting package was used to generate whole-genome heatmaps at 1 Mb resolution.

2.7 Whole-genome sequencing (WGS) and translocation analysis

Whole-genome DNA of GBC-SD and NOZ cell lines were extracted and sequenced at 80× depth. After quality control and preprocessing, the clean data were aligned to hg19 using BWA (version 0.7.17). Translocation events were identified by Delly [22] software (version: 1.0.3; command: delly call -t TRA -g hg19) and further filtered with the number of paired-end supporting the structural variant greater than or equal to 10. The resulting high-quality events were used for further analyses. We then determined the highest 100 interchromosomal interactions using the following process. Intrachromosomal bin pairs were first deleted from the matrices and then sorted based on the interaction count. Then, adjacent bin pairs were merged. When the merged adjacent bin pairs reached 100, the number “N” of bin pairs was determined. “N” was 222 in GBC-SD cells and 278 in NOZ cells. In addition, the bin pairs from Hi-C data were compared with translocations identified using the WGS data. NeoLoopFinder [23] software was used for genome reconstruction around the breakpoints of translocation events found in the WGS and Hi-C data sets (command: assemble-complex SVs –minimum-size 5000 –balance-type ICE).

2.8 Analysis of compartments A/B

Juicertools [24] was used to analyze compartmentalization through the first principal component (PC1 values, eigenvector) of Pearson’s matrix at 500 kb resolution. The whole genome was divided into two groups based on the eigenvector sign. The gene expression density with fragments per kilobase of transcript per million mapped reads per unit distance was calculated for high and low expression densities, in which the high density as a positive value was designated compartment A, and the low density as a negative value was designated compartment B.

2.9 RNA-seq analysis

Trimmomatic (version 0.39) [25] was used to trim universal Illumina adaptors and filter low-quality reads. High-quality reads were aligned to h19 using STAR (version 2.7.8a) [26]. Read counts were generated by HTseq [27] (version 0.13.5, command: -s no -a 10 -t exon -i gene_id -m intersection-nonempty -f bam). A set of R packages was used to perform downstream processing. We collected differentially expressed genes (DEGs) using DESeq2 (version 1.30.1) [28]. For samples without biological replication, DEGs were generated by edgeR [29].

2.10 ChIP-seq analysis

ChIP-seq raw data were truncated to 50 bp and aligned to hg19 using bowtie2 [30] with default parameters. Peaks were called using MACS2 [31], including H3K27ac, H3K4me3 (-B -g hs –shift 73 –broad -f BAMPE -q 0.05), and CTCF (-B -g hs -f BAMPE -q 0.05), where parameter q represents the q-value with a cutoff of 0.05. DeepTools [32] was used to calculate genome coverage with reads per kilobase of exon per million reads mapped (RPKM) normalization and generate bigWig files. For the joint analysis with A/B compartments, we merged the peak file of the tumor and control groups using the merge function of bedtools [33]. Next, we applied the multicov function to calculate the read count. Finally, edgeR package was applied in differential analysis to calculate the fold change (FC)of the Chip-seq signal.

2.11 CTCF analysis

Loops generated from juicer_tools were annotated by CTCF peaks, and the proportions of CTCF loops and non-CTCF loops were distinguished. Enrichment of the CTCF signal in loop anchors was computed and displayed using computeMatrix and plotheatmap scripts from DeepTools. The enrichment of the CTCF signal at the left and right anchors was further analyzed separately. For CTCF motif orientation analysis, we first determined the direction of each CTCF peak. The CTCF peak region sequence in hg19 was extracted using bedtools (version 2.30.0, command: getfasta). The human CTCF motif matrix (MA0139.1) was obtained from the JASPAR website, and the sequence of each peak was compared with the motif matrix using fimo software (version 5.3.3, command: –max-strand) to search for candidate binding sites with a P value cutoff equal to 1e–4. If multiple CTCF binding sites were present in one peak, the CTCF with the highest score was selected to determine the orientation of the peak. Second, we excluded loop anchors without CTCF peaks to conduct the joint analysis of loops and CTCF orientation. We also filtered out loops in which one anchor point was enriched by multiple peaks whose directions were conflicted.

2.12 Boundary analysis (identification of TADs)

Boundaries were detected by the index of insulation score as previously described [34]. Briefly, the sparse matrices generated by HiC-Pro were converted into dense ones. The total number of interactions within a square reading frame that slid along the diagonal of the matrix was assigned to the 10 kb diagonal bin. Then, the insulation score was normalized by comparison of the total number of interactions with the corresponding chromosome. Then, this score was normalized relative to the total number of interactions on the corresponding chromosome. The minima of the normalized insulation score curve represented TAD boundaries, and the bins on minima were extended by 30 kb to both sides to finally determine the range of boundaries. The insulation score of each 10 kb bin was calculated using Perl script (matrix2insulation.pl, version: 1.0.0, command: -b 10000 -im mean -bmoe 3 -nt 0.1 -v). To identify differences in the boundary strength between samples, we aligned all boundaries of each sample at the midpoint with an extension of 50 bins to both ends and compared the average insulation scores. TADs were defined as regions between two boundaries. To detect different inter/intra-TAD interactions between samples, we defined an 80 × 80 bin window with the boundary at the center. The average interactions of all windows in each sample were calculated and compared.

2.13 Loop analysis

HICCUPS [19] was used to call loops at 5, 10, and 25 kb resolutions, and different loops between samples were detected using juicer_tools. Venn diagrams were generated by R (version 3.6.1, package: eulerr). Loops were then annotated using the P and E table from the ChIP-seq analysis results. P were defined with +/−2 kb to the transcript start sites (TSS). E were defined as regions with H3K27ac peaks except P regions using bedtools (version 2.30.0, command: substract -A). The length of loops was determined by the distance between the midpoints of the left and right anchors.

2.14 Gene Ontology (GO) analysis

The loops were annotated to genes based on the genomic locations of loop ends and genes using bedtools. Genes located at cancer or normal unique loop ends were imported to the R package ClusterProfiler v4.0.5 [35] for GO analysis using the enrichGO function with the following parameters: pvalueCutoff = 0.05, qvalueCutoff = 0.05, and pAdjustMethod = BH. The default background parameters of this package were applied in which whole human reference genes (R package: org.Hs.eg.db) were used as background.

2.15 Transcription factor (TF) motif enrichment analysis

Scripts from Homer [36] were used for motif enrichment analysis. Loops with lengths greater than 100 kb were defined as long loops, and those with a length of less than 100 kb were short loops. Then, the Perl script (findMotifsGenome.pl) was used in different motif enrichment analyses with long loops as the background (command: size 2000). For E-P loop-related E and P TF enrichment analysis, E and P were classified into two groups, that is, those that were related and not related to E-P loops, with the latter used as the background.

2.16 Gene expression and E-P loop interaction analysis

Loop files of two compared cell lines were concatenated (HBEC with GBC-SD; HBCE with NOZ) and annotated in the E and P table. Interactions were evaluated as the interaction counts between bins where the left and right anchors were located. Heatmaps were generated using the R package pheatmap.

2.17 Gene expression with CTCF and H3K27ac signal

The comparisons of CTCF and H3K27ac signals at the P regions of DEGs were performed using bigWigCompare 3.5.1 (command: –skipNAs). A heatmap was plotted with computeMatrix (command: reference-point –reference Point center –sortRegions keep –missingDataAsZero –skipZeros -a 2000 -b 2000) and plotHeatmap (–sortRegions keep).

2.18 Genes selected for CTCF loop–gene expression model

We selected genes that met the following three criteria: (1) highly expressed (\left|{\mathrm{l}\mathrm{o}\mathrm{g}}_{2}FC\right|≥ 1 and P value ≤0.01) in GBC cell lines compared with HBEC; (2) highly expressed (\left|{\mathrm{l}\mathrm{o}\mathrm{g}}_{2}FC\right|≥ 1 and P value ≤ 0.05) in nine GBC tissues compared with the matched normal tissues obtained from the RNA microarray data in the Gene Expression Omnibus (GEO) database (GSE76633) [37]; (3) expression regulated via cancer-specific E-P loops. More specifically, different loops between normal cells and GBC cell lines obtained from Juicer_tools were annotated based on the P and E information. To further identify putative E, we defined genome loci with H3K27ac peaks and without H3Kme3 peaks as candidate E. Then, genes regulated by cancer-specific E-P loops and with higher expression in GBC cell lines than in HBEC cell lines were selected. Seven genes were finally identified for further analysis.

2.19 Western blot

GBC-SD and NOZ cells were lysed with radioimmunoprecipitation assay buffer (P0013C, Beyotime, China). Protein lysates were separated on a 10% polyacrylamide gel and transferred to a polyvinylidene difluoride membrane (IPVH00010, Millipore, Germany). Then, the membrane was incubated overnight at 4 °C with BMP4 (#4680, CST, USA) or GAPDH antibody (#5174, CST, USA) followed by incubation with horseradish peroxidase-conjugated IgG antibody for 1 h at RT. Finally, the chemiluminescence signal was detected by a ChemiDoc Touch System (Bio-Rad).

2.20 Quantitative real-time PCR

RNA was extracted with TRIzol reagent (#15596026, Invitrogen, USA) and reverse transcribed to cDNA using an reverse transcription reagent kit (RR047, Takara, Japan). For BMP4 mRNA expression analysis, the 2^-ΔΔCt method was applied, and GAPDH was used as an endogenous control.

2.21 Transient transfection

SiRNAs were transfected into cells using Lipofectamine 2000 transfection reagent (#11668019, Invitrogen, USA), in accordance with the manufacturer’s instructions. For BMP4 plasmid transfection, 2 µg plasmid was transfected into six-well plates using a Viafect transfection reagent (E4981, Promega, USA). After 48 h, the cells were harvested for subsequent experiments.

2.22 Cell proliferation analysis

Cell Counting Kit-8 (CCK8) assay was performed to evaluate cell proliferation. Each 96-well plate was seeded with 1000 cells. The optical density (OD) after the addition of CCK8 reagent (#40203ES60, Yeasen, China) was determined for five consecutive days.

2.23 E knockout

The templates for producing target sgRNAs were constructed in the pSpCas9(BB)-2A-Puro (PX459) V2.0 plasmid (Addgene: #62988, a gift from Feng Zhang) and confirmed by sequencing [38]. At 24 h after transfection, puromycin was added to a final concentration of 2 µg/mL for NOZ cells and 20 µg/mL for GBC-SD cells. Three days later, the DNA was extracted for PCR analysis to confirm the effectiveness of gene editing.

2.24 CTCF motif disruption

For CTCF binding site disruption, a sgRNA targeting sequence close to the CTCF motif was designed. The lentivirus was produced by cotransfection of CRISPR plasmids with VSVG and Δ8.9 plasmids into HEK293T cells using polyethylenimine. The CRISPR lentivirus was used to infect GBC-SD cells. After puromycin selection, editing efficiency was confirmed by TA cloning.

2.25 Statistical analysis

The bootstrap method was used to determine the relationship between interchromosomal interactions and translocation events identified by WGS data. Briefly, random interchromosomal interactions were selected to calculate the overlap with translocation 1000 times. The P value is equal to (the number of times random sites exceeded the top 100 sites)/1000. The same method was also used to determine the relationship between CTCF peaks and loop ends. For gene expression, CTCF, and H3K27ac signal changes within compartment switching, the Wilcoxon rank-sum test was used to compare the {\mathrm{l}\mathrm{o}\mathrm{g}}_{2}F C. The Wilcoxon rank-sum test was also used for TAD sizes, loop length, and insulation score comparison. The significance of different interactions between up/down expression genes was calculated using the Wilcoxon rank-sum test on the interaction count ratio. For cell proliferation, significance was determined using an unpaired two-sided t test on the OD values on day 5. For quantitative PCR (qPCR), an unpaired two-sided t test was used to compare delta cycle threshold (CT) values.

2.26 Data availability statement

All data can be viewed in the National Omics Data Encyclopedia datasets by pasting the accession number OEP002892.

3 Results

3.1 Global changes in genome architecture in GBC

To understand the reorganization of chromatin structure in GBC, we performed in situ Hi-C on two pairs of cancerous and benign gallbladder tissues derived from patients, HBECs, and two GBC cell lines (GBC-SD and NOZ). We combined these in situ Hi-C data sets with data from CTCF binding state (CTCF ChIP-seq), E/P activity (H3K27ac/H3K4me3 ChIP-seq), and gene expression (RNA-seq) to explore the molecular mechanisms underpinning hierarchical layers of chromatin reorganization, including compartments, TADs, and loops (Fig.1). The whole-genome contact heatmaps generated using the ICE-normalized matrices revealed stronger interchromosomal interactions in tissues than in cell lines (Fig.1 and S1A). In benign and cancer tissues, interchromosomal interactions accounted for 47%–52% of total chromatin interactions, but they decreased to 11%–17% in HBECs and cancer cells (Table S1), which suggest that tissue heterogeneity may account for strong interchromosomal interactions. Moreover, a number of strong specific interchromosomal interactions were observed in the GBC cell lines (Fig.1 and S1A). Given the previous reports that strong interchromosomal interactions are involved chromosomal translocation [9], we compared these strong interchromosomal interactions with the corresponding translocation events identified using WGS data. We found 30 and 32 interchromosomal interactions from the 100 highest interchromosomal interactions, and they were ascribed to translocation events in GBC-SD and NOZ cell lines, respectively (bootstrap P value < 0.001, Fig.1 and S1B, Table S2). To further explore the specific genes that may be influenced by aberrant chromosomal interactions, we rearranged genome translocation events using NeoLoopFinder software and found strong interactions within the PTPRD locus on chromosome 9 joined with chromosome 2 and the ZNF423 locus on chromosome 16 joined with chromosome 3 (Fig.1 and S1C, Table S3). PTPRD is a member of the protein tyrosine phosphatase family and regulates a variety of cellular processes, including cancer progression [39]. ZNF423 and its mouse ortholog Zfp423 are critical transcriptional modulators in neuroblastoma and leukemia [40]. Altogether, we integrated high-resolution chromatin interaction maps and epigenetic data on gallbladder tissues and cell lines and observed the association of strong interchromosomal interactions with chromosomal translocations.

3.2 Compartments and TADs are reorganized in GBC

To determine whether changes in large-scale compartment organization occur in GBC, we quantified the A/B compartments using an eigenvector-based method at 500 kb resolution. The comparisons revealed a significant difference between GBC and normal controls in A/B compartmentalization. Although most of the genome regions retained their compartment identities, a significant fraction of the genome switched from compartment A to B (5%–12%) or from compartment B to A (6%–15%) between normal and cancer (Fig.2 and 2B, Table S4). In the region of chromosome 4, a 30 Mb region switched compartments between normal and cancer cells (Fig. S2A). Notably, when compartments A/B switched, the expression levels of genes decreased considerably relative to the genes in the A/A compartment in GBC-SD or NOZ cells versus HBECs (Fig.2 and S2B). Conversely, when compartment B/A switched, the gene expression levels increased. CTCF and H3K27ac signals exhibited similar changes with compartment switching in these cells (Fig.2, 2E, and S2C). In addition, a majority of CTCF peaks (~70%) were located in compartment A, and the signal of CTCF peaks in compartment A was higher than that in compartment B (Fig. S2D and S2E). However, DEGs were not always associated with compartment switching, with only 11%–27% of DEGs located in the compartment-switching region (Fig. S2F). To further elucidate the relationship between chromatin structural changes and gene expression changes in GBC, we analyzed the TAD organization in normal and cancer tissues and cell lines. Larger but slightly fewer TADs were identified in cancer lines than in HBECs. However, the opposite relationship was observed in cancer and normal tissues (Fig.2 and 2G, Table S5). Notably, the insulation scores for the TAD boundaries were higher in one GBC tissue and two GBC cell lines than in benign tissue and HBECs (Fig.2 and S2G, Table S6), which indicates weakened boundaries in GBC. To further interpret the functional implications of these alterations in insulation scores, we generated pile-up heatmaps of all TADs from cancer and normal cells by aligning TAD boundaries and then analyzed the differences in their inter/intra-TAD interactions. The results revealed that inter-TAD interactions increased but intra-TAD interactions decreased in GBC cell lines versus HBECs (Fig.2). Similarly, GBC-1 of two cancerous tissues showed higher inter-TAD interactions but lower intra-TAD interactions than benign tissues (Fig. S2H). Altogether, compartment and TAD organization exhibited considerable differences between GBC and normal samples. The TAD boundaries in cancer tissues and cell lines were weakened relative to their counterparts in normal tissues and cell lines, and this result may have profound implications for gene expressions.

3.3 Altered chromatin loops are associated with genes involved in cancer progression

Given that TADs are established by large-scale chromatin loops [41,42], we subsequently analyzed whether the different loop structures contribute to the alteration of TADs. The majority of loops over 50% were shared between HBECs and GBC cell lines, and 21%–42% of loops were unique to each cell line. Similarly, more than 50% of loops were common between tumor and normal tissues (Fig.3 and S3A, Table S7). In addition, as shown in Fig.3 and S3B, cancer-specific loops were significantly shorter than normal-specific or common loops in cancer cells/tissue versus normal cells/tissue. These loop features differed from the TAD size characteristics between tissues and cell lines. By analyzing the distribution of cancer-specific, normal-specific, and common loops with different lengths in tissues and cells, we discovered that a higher proportion of loops < 100 kb was located in cancer-specific loops than in other loops (Fig.3 and S3C). These data suggest that these short loops may be primarily associated with cancer progression. To assess whether cancer-specific short loops harbor distinct transcription factor binding motifs, we analyzed the enrichment of transcription factor motifs in the anchor locus of short loops and observed that six transcription factor motifs were enriched in these short loops compared with those in long loops in cancers versus normal controls, in which Sp1 transcription factor (SP1) [43] and activating transcription factor 1 (ATF1) [44] are known as the activated transcription factors (Fig. S3D). To further explore the potential biological process possibly ascribed to cancer-specific loops, we performed GO enrichment analysis of genes associated with different groups of loops. Biological processes related to development were enriched in normal-specific loops, whereas cancer-associated biological processes, including epithelial–mesenchymal transition (EMT), chromatin remodeling, and regulation of MAP kinase activity, were identified in cancer-specific loops [45–47] (Fig.3 and S3E, Table S8). SRY-box transcription factor 9 (SOX9) was a typical gene example in the chromatin remodeling process, and it was located at the anchor of a loop unique to NOZ cells, rather than in HBECs (Fig.3). Altogether, we discovered that the distinct chromatin loop landscapes between GBC and normal controls may be associated with genotypic alterations in the occurrence and development of GBC.

3.4 CTCF-mediated loop formation

Next, we explored the underlying mechanism for the differences in loop structures. Given that CTCF blocks the cohesion-mediated loop extrusion process at chromatin loop anchors, we analyzed the CTCF ChIP-seq data to define the CTCF signature in the loops in GBC and normal cells. Consistent with previous studies [19,48], we found that the majority of loops were mediated by CTCF in HBEC and GBC cell lines, of which 83%–91% loops contained CTCF binding sites in at least one anchor, and 42%–60% loops had CTCF binding sites in both anchors (bootstrap P value < 0.001, Fig.4). In addition, compared with the random region of the chromosome, CTCF signals were enriched at the loop anchors to similar extents in normal and cancer cells (Fig.4). We further classified the loops into four types based on the polarity of CTCF binding (forward-reverse, reverse-forward, forward-forward, and reverse-reverse) [4,49]. In GBC and HBEC cells, ~79% of loops were anchored by CTCF binding sites in a convergent orientation, ~20% of loops in the same orientation, and ~1% of loops in divergent direction (Fig.4). Interestingly, CTCF was more likely to bind to either the left side of the left anchor midpoint or the right side of the right anchor midpoint (Fig.4). In an example of chromatin loops on chr10, in addition to two loops commonly located in normal cells and NOZ cells, a specific loop connected by two convergent CTCF binding sites was identified only in NOZ cells (Fig.4). This cancer-specific loop, which was associated with a newly emerged CTCF peak in one of the loop anchors, suggests that changes in CTCF genome-wide occupancy may underlie the reorganization of chromatin loop landscapes in GBC cells. Collectively, our data suggest that the cancer-specific chromatin loops associated with the convergent CTCF binding sites may provide a structural basis for cancer-associated gene regulation.

3.5 CTCF-mediated E-P interaction loop regulates cancer-related gene expression

To identify CTCF-mediated loops as gene regulatory elements, we specifically focused on the relationship between CTCF, E-P loop interactions, and gene expression. We analyzed E-P interaction levels using H3K27ac peaks (excluding the +/− 2 kb region of the TSS) as the putative E region in HBECs and GBC-SD/NOZ cells. First, we analyzed the proportion of E-P loops in the total loops and observed that ~20% of loops harbored E and P elements in all cell lines (Fig. S4A). In addition, gene expression, CTCF occupancy, and H3K27ac signal were positively correlated with gene expression levels in cancer cells (Fig.5 and S4B). Furthermore, CTCF and CTCFL motifs were enriched in the E and P loci of E-P loops compared with other E and P in GBC cells and HBECs (Fig. S4C). In addition, the gene expression levels in cancer cells were positively correlated with E-P interaction frequencies (Fig.5 and S4D), which suggests that E-P interactions may play important roles in directing cancer-specific gene expression programs. To further define the cancer-related genes potentially controlled by E-P loops, we analyzed the gene expression profiles of the GBC cohort (nine GBC samples, GEO accession number: GSE76633), GBC-SD, and NOZ cells. As a result, cross-combined assessment of these gene data sets revealed seven altered genes, including BMP4, KRT19, and others (Fig.5 and 5D, Table S9). FOXA2 gene was elevated in two GBC cell lines. BMP4 is a member of the TGF-β superfamily that mediates cell proliferation and EMT [50], and KRT19 is upregulated and correlated with the overall survival in lung cancer patients [51]. FOXA2 is elevated in patients with Barrett’s metaplasia, dysplasia, and adenocarcinoma [52]. By integrating RNA-seq, ChIP-seq (CTCF/H3K27ac/H3K4me3), loop interactions, and Hi-C heatmap, we observed that all three genes (BMP4, KRT19, and FOXA2) were associated with cancer-unique E-P loops in GBC-SD and NOZ cell lines (Fig.5 and S4E). Only the H3K27ac signal enriched at the end of the loops was annotated as the putative E of BMP4, KRT19, or FOXA2. H3K27ac and H3K4me3 signals enriched at one loop end were defined as P. Moreover, a unique CTCF peak downstream of the right end was found in GBC-SD cells. When we knocked down BMP4 levels using siRNA in two GBC cell lines, the proliferation of cancer cells was substantially inhibited (Fig.5 and 5G). These results suggest the need to develop a model in which the CTCF-mediated E-P interaction loop in cancer regulates gene expression for cancer cell development.

3.6 Deletion of the E or CTCF binding site downregulates gene expression

To further verify that BMP4, KRT19, and FOXA2 are regulated by cancer-specific E-P loop interactions, we deleted the putative E in GBC-SD and NOZ cells using the CRISPR/Cas9 system. The efficiency of E deletion via two sgRNAs was confirmed by Sanger sequencing and nucleic acid electrophoresis (Fig.6 and S5A). Compared with the wild-type (WT) GBC-SD and NOZ cells, the expressions of all three genes in edited cells decreased (Fig.6, 6C and S5B). To subsequently evaluate the cellular functional changes, we assessed cell proliferation in BMP4-associated putative E-edited and WT cells by CCK8 assay. Concordant with the knockdown of BMP4, the E-deleted cells also exhibited defects in cell proliferation. By contrast, the overexpression of BMP4 in the edited cells fully rescued the cell proliferation defects (Fig.6). Altogether, these data show that the E-P loop interaction can induce cancer-related gene expression and promote tumor cell proliferation in GBC. In addition, to further investigate whether BMP4 is regulated by CTCFs that mediate the formation of the E-P loop, we reconfirmed the specific CTCF binding downstream of the BMP4-associated E by ChIP‒qPCR and disrupted the CTCF binding site using CRISPR/Cas9 (Fig.7 and S5C). CTCF-edited cells with 91% indels in CTCF binding sites were obtained (Fig.7). The disruption of CTCF binding sites decreased the expression of BMP4 at the mRNA and protein levels (Fig.7 and 7D, respectively). Overall, our data support the notion that the CTCF-mediated E-P loop regulates the abnormal expression of the BMP4 gene and contributes to cancer progression.

4 Discussion

Structurally altered genomes can affect gene expression and functions, which leads to aberrant individual genotypes and phenotypes in a broad spectrum of human diseases [53–55]. Growing evidence has demonstrated that genomic organization is a multistep process involving CTCF-loop formation by distant locus contact in linear chromosomes [19,56–58], high self-interaction of TADs composed of multiple loops [41,59], and large-scale A/B compartments compacted with interactive chromatin [8,60,61]. Despite the previous intense focus on TAD reorganization, the structural and functional relationships among chromosomal compartments, CTCF loops, and gene expressions remain to be established in cancer. Hence, we systematically investigated the topological landscape of individual genomic components and functional changes in GBC tissues and cell lines by means of in situ Hi-C, RNA-seq, and ChIP-seq methodologies. In the comparison between GBC tissues/cancer cell lines and benign tissues/normal epithelial cells, we observed that GBC harbored unique cancer loops that were associated with chromatin remodeling and EMT activation, whereas controls consisted of loops regulating normal development. In addition, gene expression was correlated with CTCF levels and E-P interactions, which highlights the importance of CTCF-mediated E-P loops for gene expression in cancer. BMP4 gene expression was upregulated by the CTCF-mediated E-P loop in GBC cell lines, which was absent in normal controls. Furthermore, in GBC, increased inter-TAD interactions with weaker TAD boundaries may facilitate proximal TADs to create new loops. In addition, CTCF peaks were mostly distributed in compartment A, and increased CTCF levels and gene expressions were correlated with the switching of compartment B to A. Collectively, our findings suggest that reorganization of the 3D genome architecture at multiple levels contributes to genotypic and phenotypic abnormalities during the development of GBC.

In the analysis of inter- and intra-chromosomal interaction ratios, the ratios in our patient tissues were 47%–52% and 48%–53%, respectively, and the values in cell lines were 11%–17% and 83%–89%, which suggest that tissue heterogeneity may account for the strong interchromosomal interactions. In concert with our findings, a number of tissue and cell line studies exhibited the elevated intrachromosomal interaction ratio in tissues. The average inter/intrachromosomal interaction ratio is 37%/63% in cardiac tissue [62] and 48%/52% in human lung tissues [63]. However, the ratio reaches 14%/86% in pancreatic cell lines [64] and 17%/83% in prostate cancer cell lines [65]. Thus, all the studies should explore the molecular mechanisms underlying the different inter/intrachromosomal interaction ratios between tissues and cultured cells.

A recent study on colorectal cancer (CRC) reported that changes in the DNA loops between MYC and its E can promote the development of colon cancer [66]. In addition, a CRC risk variant approximately 300 kb away within an E of MYC facilitated loop formation [67]. In line with these findings, we defined the regulatory mechanisms of the cancer-related genes, namely, BMP4, KRT19, and FOXA2, driven by an E-P loop in GBC-SD and NOZ cell lines. The disruption of individual putative E via the CRISPR/Cas9 system resulted in downregulation of gene expressions and inhibition of cell proliferation in BMP4-E-edited cells. The individual anchors enriched in the binding of two convergent CTCF proteins at both ends of the loop secured loop integrity and function. This result was in line with a reported evidence indicating that 46% (4322/9948) of loops were mediated by two CTCF binding sites and 92% of motif pairs were convergent [19]. Disruption of the specific CTCF motif around the BMP4-related E in GBC-SD cells resulted in the downregulation of BMP4 expression. However, only one CTCF binding site was identified at the cancer-specific E-P loop at the BMP4 locus in NOZ cells. A recent study may explain this discrepancy, which revealed that P exhibited a higher CTCF binding affinity than E and that deletion of this CTCF peak had a role in determining finer-scale E-P choice and engagement [68]. Another report suggested that DNA methylation of the CTCF site can alter CTCF binding and loop formation [69]. In our system, gene modification may participate in this event but not the genetic mutation that has been excluded (Fig. S5D). Notably, we observed that CTCF peaks were not always at loop ends; instead, they were located on both sides of the ends. This phenomenon may be related to the limited resolution of the interaction matrix. In addition, CTCF was more likely to bind to either the left side of the left anchor midpoint or the right side of the right anchor midpoint. In summary, we have offered substantial pieces of evidence that reveal the gene expression mechanisms regulated by CTCF-mediated E-P interaction loops in GBC.

Our current study focused on GBC cell lines and two GBC patients and unveiled the topological and functional alterations of the cancer genome. Although the overall outcomes support our gene regulatory model driven by CTCF-mediated E-P loops, two human samples yielded several inconsistencies that limit our knowledge. The insulation score of boundaries changed in one GBC patient but not in another. Similarly, one patient possessed more cancer-specific short loops, but the other did not. Therefore, the enrollment of a substantially large number of GBC patients is essential to eliminate divergence and confirm our findings.

Given the rapid development of detection systems, such as in situ Hi-C and ChIP-seq, the 3D spatial genome organization with gene regulatory elements, including E, P, and insulators, in diseases has increasingly received remarkable attention, with the expectation of offering novel insights into the topological organization of individual genomic components. To date, a number of molecular mechanisms, such as the regulation of short loop formations, differences in inter/intra-TAD interactions, and compartment A/B switching, remain to be clarified. Nevertheless, the current study has demonstrated that topological changes in distinct genomic layers give rise to the oncogenic signature and drive cancer-specific gene expression in cancer development.

References

[1] Kloetgen A, Thandapani P, Ntziachristos P, Ghebrechristos Y, Nomikou S, Lazaris C, Chen X, Hu H, Bakogianni S, Wang J, Fu Y, Boccalatte F, Zhong H, Paietta E, Trimarchi T, Zhu Y, Van Vlierberghe P, Inghirami GG, Lionnet T, Aifantis I, Tsirigos A. Three-dimensional chromatin landscapes in T cell acute lymphoblastic leukemia. Nat Genet 2020; 52(4): 388–400

[2] Kaaij LJT, Mohn F, van der Weide RH, de Wit E, Bühler M. The ChAHP complex counteracts chromatin looping at CTCF sites that emerged from SINE expansions in mouse. Cell 2019; 178(6): 1437–1451.e14

[3] Zha J, Lai Q, Deng M, Shi P, Zhao H, Chen Q, Wu H, Xu B. Disruption of CTCF boundary at HOXA locus promote BET inhibitors’ therapeutic sensitivity in acute myeloid leukemia. Stem Cell Rev Rep 2020; 16(6): 1280–1291

[4] Guo Y, Xu Q, Canzio D, Shou J, Li J, Gorkin DU, Jung I, Wu H, Zhai Y, Tang Y, Lu Y, Wu Y, Jia Z, Li W, Zhang MQ, Ren B, Krainer AR, Maniatis T, Wu Q. CRISPR inversion of CTCF sites alters genome topology and enhancer/promoter function. Cell 2015; 162(4): 900–910

[5] Pope BD, Ryba T, Dileep V, Yue F, Wu W, Denas O, Vera DL, Wang Y, Hansen RS, Canfield TK, Thurman RE, Cheng Y, Gülsoy G, Dennis JH, Snyder MP, Stamatoyannopoulos JA, Taylor J, Hardison RC, Kahveci T, Ren B, Gilbert DM. Topologically associating domains are stable units of replication-timing regulation. Nature 2014; 515(7527): 402–405

[6] Yu M, Ren B. The three-dimensional organization of mammalian genomes. Annu Rev Cell Dev Biol 2017; 33(1): 265–289

[7] Misteli T. The self-organizing genome: principles of genome architecture and function. Cell 2020; 183(1): 28–45

[8] Lieberman-Aiden E, van Berkum NL, Williams L, Imakaev M, Ragoczy T, Telling A, Amit I, Lajoie BR, Sabo PJ, Dorschner MO, Sandstrom R, Bernstein B, Bender MA, Groudine M, Gnirke A, Stamatoyannopoulos J, Mirny LA, Lander ES, Dekker J. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science 2009; 326(5950): 289–293

[9] Wu P, Li T, Li R, Jia L, Zhu P, Liu Y, Chen Q, Tang D, Yu Y, Li C. 3D genome of multiple myeloma reveals spatial genome disorganization associated with copy number variations. Nat Commun 2017; 8(1): 1937

[10] Pandey A, Stawiski EW, Durinck S, Gowda H, Goldstein LD, Barbhuiya MA, Schröder MS, Sreenivasamurthy SK, Kim SW, Phalke S, Suryamohan K, Lee K, Chakraborty P, Kode V, Shi X, Chatterjee A, Datta K, Khan AA, Subbannayya T, Wang J, Chaudhuri S, Gupta S, Shrivastav BR, Jaiswal BS, Poojary SS, Bhunia S, Garcia P, Bizama C, Rosa L, Kwon W, Kim H, Han Y, Yadav TD, Ramprasad VL, Chaudhuri A, Modrusan Z, Roa JC, Tiwari PK, Jang JY, Seshagiri S. Integrated genomic analysis reveals mutated ELF3 as a potential gallbladder cancer vaccine candidate. Nat Commun 2020; 11(1): 4225

[11] Zhang L, Miao R, Zhang X, Chen W, Zhou Y, Wang R, Zhang R, Pang Q, Xu X, Liu C. Exploring the diagnosis markers for gallbladder cancer based on clinical data. Front Med 2015; 9(3): 350–355

[12] Aloia TA, Járufe N, Javle M, Maithel SK, Roa JC, Adsay V, Coimbra FJ, Jarnagin WR. Gallbladder cancer: expert consensus statement. HPB (Oxford) 2015; 17(8): 681–690

[13] Rustagi T, Dasanu CA. Risk factors for gallbladder cancer and cholangiocarcinoma: similarities, differences and updates. J Gastrointest Cancer 2012; 43(2): 137–147

[14] Li M, Liu F, Zhang F, Zhou W, Jiang X, Yang Y, Qu K, Wang Y, Ma Q, Wang T, Bai L, Wang Z, Song X, Zhu Y, Yuan R, Gao Y, Liu Y, Jin Y, Li H, Xiang S, Ye Y, Zhang Y, Jiang L, Hu Y, Hao Y, Lu W, Chen S, Gu J, Zhou J, Gong W, Zhang Y, Wang X, Liu X, Liu C, Liu H, Liu Y, Liu Y. Genomic ERBB2/ERBB3 mutations promote PD-L1-mediated immune escape in gallbladder cancer: a whole-exome sequencing analysis. Gut 2019; 68(6): 1024–1033

[15] Li M, Zhang Z, Li X, Ye J, Wu X, Tan Z, Liu C, Shen B, Wang XA, Wu W, Zhou D, Zhang D, Wang T, Liu B, Qu K, Ding Q, Weng H, Ding Q, Mu J, Shu Y, Bao R, Cao Y, Chen P, Liu T, Jiang L, Hu Y, Dong P, Gu J, Lu W, Shi W, Lu J, Gong W, Tang Z, Zhang Y, Wang X, Chin YE, Weng X, Zhang H, Tang W, Zheng Y, He L, Wang H, Liu Y, Liu Y. Whole-exome and targeted gene sequencing of gallbladder carcinoma identifies recurrent mutations in the ErbB pathway. Nat Genet 2014; 46(8): 872–876

[16] Hu YP, Wu ZB, Jiang L, Jin YP, Li HF, Zhang YJ, Ma Q, Ye YY, Wang Z, Liu YC, Chen HZ, Liu YB. STYK1 promotes cancer cell proliferation and malignant transformation by activating PI3K-AKT pathway in gallbladder carcinoma. Int J Biochem Cell Biol 2018; 97: 16–27

[17] Jin YP, Hu YP, Wu XS, Wu YS, Ye YY, Li HF, Liu YC, Jiang L, Liu FT, Zhang YJ, Hao YJ, Liu XY, Liu YB. miR-143-3p targeting of ITGA6 suppresses tumour growth and angiogenesis by downregulating PLGF expression via the PI3K/AKT pathway in gallbladder carcinoma. Cell Death Dis 2018; 9(2): 182

[18] Li H, Jin Y, Hu Y, Jiang L, Liu F, Zhang Y, Hao Y, Chen S, Wu X, Liu Y. The PLGF/c-MYC/miR-19a axis promotes metastasis and stemness in gallbladder cancer. Cancer Sci 2018; 109(5): 1532–1544

[19] Rao SS, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson JT, Sanborn AL, Machol I, Omer AD, Lander ES, Aiden EL. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell 2014; 159(7): 1665–1680

[20] Servant N, Varoquaux N, Lajoie BR, Viara E, Chen CJ, Vert JP, Heard E, Dekker J, Barillot E. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol 2015; 16(1): 259

[21] Imakaev M, Fudenberg G, McCord RP, Naumova N, Goloborodko A, Lajoie BR, Dekker J, Mirny LA. Iterative correction of Hi-C data reveals hallmarks of chromosome organization. Nat Methods 2012; 9(10): 999–1003

[22] Rausch T, Zichner T, Schlattl A, Stütz AM, Benes V, Korbel JO. DELLY: structural variant discovery by integrated paired-end and split-read analysis. Bioinformatics 2012; 28(18): i333–i339

[23] Wang X, Xu J, Zhang B, Hou Y, Song F, Lyu H, Yue F. Genome-wide detection of enhancer-hijacking events from chromatin interaction data in rearranged genomes. Nat Methods 2021; 18(6): 661–668

[24] Durand NC, Shamim MS, Machol I, Rao SS, Huntley MH, Lander ES, Aiden EL. Juicer provides a one-click system for analyzing loop-resolution Hi-C experiments. Cell Syst 2016; 3(1): 95–98

[25] Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014; 30(15): 2114–2120

[26] Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2013; 29(1): 15–21

[27] Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 2015; 31(2): 166–169

[28] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014; 15(12): 550

[29] Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010; 26(1): 139–140

[30] Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods 2012; 9(4): 357–359

[31] Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, Liu XS. Model-based analysis of ChIP-Seq (MACS). Genome Biol 2008; 9(9): R137

[32] Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, Heyne S, Dündar F, Manke T. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 2016; 44(W1): W160–5

[33] Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 2010; 26(6): 841–842

[34] Crane E, Bian Q, McCord RP, Lajoie BR, Wheeler BS, Ralston EJ, Uzawa S, Dekker J, Meyer BJ. Condensin-driven remodelling of X chromosome topology during dosage compensation. Nature 2015; 523(7559): 240–244

[35] Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012; 16(5): 284–287

[36] Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B-cell identities. Mol Cell 2010; 38(4): 576–589

[37] Wu XS, Wang F, Li HF, Hu YP, Jiang L, Zhang F, Li ML, Wang XA, Jin YP, Zhang YJ, Lu W, Wu WG, Shu YJ, Weng H, Cao Y, Bao RF, Liang HB, Wang Z, Zhang YC, Gong W, Zheng L, Sun SH, Liu YB. LncRNA-PAGBC acts as a microRNA sponge and promotes gallbladder tumorigenesis. EMBO Rep 2017; 18(10): 1837–1853

[38] Ran FA, Hsu PD, Wright J, Agarwala V, Scott DA, Zhang F. Genome engineering using the CRISPR‒Cas9 system. Nat Protoc 2013; 8(11): 2281–2308

[39] Ortiz B, Fabius AW, Wu WH, Pedraza A, Brennan CW, Schultz N, Pitter KL, Bromberg JF, Huse JT, Holland EC, Chan TA. Loss of the tyrosine phosphatase PTPRD leads to aberrant STAT3 activation and promotes gliomagenesis. Proc Natl Acad Sci USA 2014; 111(22): 8149–8154

[40] Harder L, Puller AC, Horstmann MA. ZNF423: transcriptional modulation in development and cancer. Mol Cell Oncol 2014; 1(3): e969655

[41] Hnisz D, Day DS, Young RA. Insulated neighborhoods: structural and functional units of mammalian gene control. Cell 2016; 167(5): 1188–1200

[42] Wutz G, Várnai C, Nagasaka K, Cisneros DA, Stocsits RR, Tang W, Schoenfelder S, Jessberger G, Muhar M, Hossain MJ, Walther N, Koch B, Kueblbeck M, Ellenberg J, Zuber J, Fraser P, Peters JM. Topologically associating domains and chromatin loops depend on cohesin and are regulated by CTCF, WAPL, and PDS5 proteins. EMBO J 2017; 36(24): 3573–3599

[43] Song J, Nabeel-Shah S, Pu S, Lee H, Braunschweig U, Ni Z, Ahmed N, Marcon E, Zhong G, Ray D, Ha KCH, Guo X, Zhang Z, Hughes TR, Blencowe BJ, Greenblatt JF. Regulation of alternative polyadenylation by the C2H2-zinc-finger protein Sp1. Mol Cell 2022; 82(17): 3135–3150.e9

[44] Mayr B, Montminy M. Transcriptional regulation by the phosphorylation-dependent factor CREB. Nat Rev Mol Cell Biol 2001; 2(8): 599–609

[45] Dongre A, Weinberg RA. New insights into the mechanisms of epithelial-mesenchymal transition and implications for cancer. Nat Rev Mol Cell Biol 2019; 20(2): 69–84

[46] Burotto M, Chiou VL, Lee JM, Kohn EC. The MAPK pathway across different malignancies: a new perspective. Cancer 2014; 120(22): 3446–3456

[47] Dawson MA, Kouzarides T. Cancer epigenetics: from mechanism to therapy. Cell 2012; 150(1): 12–27

[48] Rao SSP, Huang SC, Glenn St Hilaire B, Engreitz JM, Perez EM, Kieffer-Kwon KR, Sanborn AL, Johnstone SE, Bascom GD, Bochkov ID, Huang X, Shamim MS, Shin J, Turner D, Ye Z, Omer AD, Robinson JT, Schlick T, Bernstein BE, Casellas R, Lander ES, Aiden EL. Cohesin loss eliminates all loop domains. Cell 2017; 171(2): 305–320.e24

[49] Nora EP, Caccianini L, Fudenberg G, So K, Kameswaran V, Nagle A, Uebersohn A, Hajj B, Saux AL, Coulon A, Mirny LA, Pollard KS, Dahan M, Bruneau BG. Molecular basis of CTCF binding polarity in genome folding. Nat Commun 2020; 11(1): 5612

[50] Martínez VG, Rubio C, Martínez-Fernández M, Segovia C, López-Calderón F, Garín MI, Teijeira A, Munera-Maravilla E, Varas A, Sacedón R, Guerrero F, Villacampa F, de la Rosa F, Castellano D, López-Collazo E, Paramio JM, Vicente Á, Dueñas M. BMP4 induces M2 macrophage polarization and favors tumor progression in bladder cancer. Clin Cancer Res 2017; 23(23): 7388–7399

[51] Yuan X, Yi M, Dong B, Chu Q, Wu K. Prognostic significance of KRT19 in lung squamous cancer. J Cancer 2021; 12(4): 1240–1248

[52] Wang DH, Tiwari A, Kim ME, Clemons NJ, Regmi NL, Hodges WA, Berman DM, Montgomery EA, Watkins DN, Zhang X, Zhang Q, Jie C, Spechler SJ, Souza RF. Hedgehog signaling regulates FOXA2 in esophageal embryogenesis and Barrett’s metaplasia. J Clin Invest 2014; 124(9): 3767–3780

[53] PCAWG Transcriptome Core Group; Calabrese C, Davidson NR, Demircioğlu D, Fonseca NA, He Y, Kahles A, Lehmann KV, Liu F, Shiraishi Y, Soulette CM, Urban L, Greger L, Li S, Liu D, Perry MD, Xiang Q, Zhang F, Zhang J, Bailey P, Erkek S, Hoadley KA, Hou Y, Huska MR, Kilpinen H, Korbel JO, Marin MG, Markowski J, Nandi T, Pan-Hammarström Q, Pedamallu CS, Siebert R, Stark SG, Su H, Tan P, Waszak SM, Yung C, Zhu S, Awadalla P, Creighton CJ, Meyerson M, Ouellette BFF, Wu K, Yang H; PCAWG Transcriptome Working Group; Brazma A, Brooks AN, Göke J, Rätsch G, Schwarz RF, Stegle O, Zhang Z; PCAWG Consortium. Genomic basis for RNA alterations in cancer. Nature 2020; 578(7793): 129–136

[54] Quigley DA, Dang HX, Zhao SG, Lloyd P, Aggarwal R, Alumkal JJ, Foye A, Kothari V, Perry MD, Bailey AM, Playdle D, Barnard TJ, Zhang L, Zhang J, Youngren JF, Cieslik MP, Parolia A, Beer TM, Thomas G, Chi KN, Gleave M, Lack NA, Zoubeidi A, Reiter RE, Rettig MB, Witte O, Ryan CJ, Fong L, Kim W, Friedlander T, Chou J, Li H, Das R, Li H, Moussavi-Baygi R, Goodarzi H, Gilbert LA, Lara PN Jr, Evans CP, Goldstein TC, Stuart JM, Tomlins SA, Spratt DE, Cheetham RK, Cheng DT, Farh K, Gehring JS, Hakenberg J, Liao A, Febbo PG, Shon J, Sickler B, Batzoglou S, Knudsen KE, He HH, Huang J, Wyatt AW, Dehm SM, Ashworth A, Chinnaiyan AM, Maher CA, Small EJ, Feng FY. Genomic hallmarks and structural variation in metastatic prostate cancer. Cell 2018; 174(3): 758–769.e9

[55] Misteli T. Higher-order genome organization in human disease. Cold Spring Harb Perspect Biol 2010; 2(8): a000794

[56] Zuin J, Dixon JR, van der Reijden MI, Ye Z, Kolovos P, Brouwer RW, van de Corput MP, van de Werken HJ, Knoch TA, van IJcken WF, Grosveld FG, Ren B, Wendt KS. Cohesin and CTCF differentially affect chromatin architecture and gene expression in human cells. Proc Natl Acad Sci USA 2014; 111(3): 996–1001

[57] Vietri Rudan M, Barrington C, Henderson S, Ernst C, Odom DT, Tanay A, Hadjur S. Comparative Hi-C reveals that CTCF underlies evolution of chromosomal domain architecture. Cell Rep 2015; 10(8): 1297–1309

[58] Ong CT, Corces VG. CTCF: an architectural protein bridging genome topology and function. Nat Rev Genet 2014; 15(4): 234–246

[59] Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, Hu M, Liu JS, Ren B. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature 2012; 485(7398): 376–380

[60] Fortin JP, Hansen KD. Reconstructing A/B compartments as revealed by Hi-C using long-range correlations in epigenetic data. Genome Biol 2015; 16(1): 180

[61] Rowley MJ, Nichols MH, Lyu X, Ando-Kuri M, Rivera ISM, Hermetz K, Wang P, Ruan Y, Corces VG. Evolutionarily conserved principles predict 3D chromatin organization. Mol Cell 2017; 67(5): 837–852.e7

[62] Rosa-Garrido M, Chapski DJ, Schmitt AD, Kimball TH, Karbassi E, Monte E, Balderas E, Pellegrini M, Shih TT, Soehalim E, Liem D, Ping P, Galjart NJ, Ren S, Wang Y, Ren B, Vondriska TM. High-resolution mapping of chromatin conformation in cardiac myocytes reveals structural remodeling of the epigenome in heart failure. Circulation 2017; 136(17): 1613–1625

[63] Li T, Li R, Dong X, Shi L, Lin M, Peng T, Wu P, Liu Y, Li X, He X, Han X, Kang B, Wang Y, Liu Z, Chen Q, Shen Y, Feng M, Wang X, Wu D, Wang J, Li C. Integrative analysis of genome, 3D genome, and transcriptome alterations of clinical lung cancer samples. Genom Proteom Bioinfor 2021; 19(5): 741–753

[64] Ren B, Yang J, Wang C, Yang G, Wang H, Chen Y, Xu R, Fan X, You L, Zhang T, Zhao Y. High-resolution Hi-C maps highlight multiscale 3D epigenome reprogramming during pancreatic cancer metastasis. J Hematol Oncol 2021; 14(1): 120

[65] Luo Z, Rhie SK, Lay FD, Farnham PJ. A prostate cancer risk element functions as a repressive loop that regulates HOXA13. Cell Rep 2017; 21(6): 1411–1417

[66] Xiang JF, Yin QF, Chen T, Zhang Y, Zhang XO, Wu Z, Zhang S, Wang HB, Ge J, Lu X, Yang L, Chen LL. Human colorectal cancer-specific CCAT1-L lncRNA regulates long-range chromatin interactions at the MYC locus. Cell Res 2014; 24(5): 513–531

[67] Pomerantz MM, Ahmadiyeh N, Jia L, Herman P, Verzi MP, Doddapaneni H, Beckwith CA, Chan JA, Hills A, Davis M, Yao K, Kehoe SM, Lenz HJ, Haiman CA, Yan C, Henderson BE, Frenkel B, Barretina J, Bass A, Tabernero J, Baselga J, Regan MM, Manak JR, Shivdasani R, Coetzee GA, Freedman ML. The 8q24 cancer risk variant rs6983267 shows long-range interaction with MYC in colorectal cancer. Nat Genet 2009; 41(8): 882–884

[68] Oh S, Shao J, Mitra J, Xiong F, D’Antonio M, Wang R, Garcia-Bassets I, Ma Q, Zhu X, Lee JH, Nair SJ, Yang F, Ohgi K, Frazer KA, Zhang ZD, Li W, Rosenfeld MG. Enhancer release and retargeting activates disease-susceptibility genes. Nature 2021; 595(7869): 735–740

[69] Nanavaty V, Abrash EW, Hong C, Park S, Fink EE, Li Z, Sweet TJ, Bhasin JM, Singuri S, Lee BH, Hwang TH, Ting AH. DNA methylation regulates alternative polyadenylation via CTCF and the cohesin complex. Mol Cell 2020; 78(4): 752–764.e6

RIGHTS & PERMISSIONS

Higher Education Press