Spatial 3D genome organization reveals intratumor heterogeneity in primary glioblastoma samples

Abstract

Glioblastoma (GBM) is the most prevalent malignant brain tumor with poor prognosis. Although chromatin intratumoral heterogeneity is a characteristic feature of GBM, most current studies are conducted at a single tumor site. To investigate the GBM-specific 3D genome organization and its heterogeneity, we conducted Hi-C experiments in 21 GBM samples from nine patients, along with three normal brain samples. We identified genome subcompartmentalization and chromatin interactions specific to GBM, as well as extensive intertumoral and intratumoral heterogeneity at these levels. We identified copy number variants (CNVs) and structural variations (SVs) and demonstrated how they disrupted 3D genome structures. SVs could not only induce enhancer hijacking but also cause the loss of enhancers to the same gene, both of which contributed to gene dysregulation. Our findings provide insights into the GBM-specific 3D genome organization and the intratumoral heterogeneity of this organization and open avenues for understanding this devastating disease.

Fig. 1.. Study design.

(A) Study design including the number of patients and samples in this study and location of each spatially mapped samples within the brain. The yellow region is the whole tumor (T2-weighted hyperintense region) and the green region is the contrast-enhancing lesion (T1-weighted post-contrast). People icon was generated from BioRender. (B) Relative distance from centroid of each sample in P530, P524, and P529. Relative distance from centroid is defined as distance to tumor centroid/(distance to tumor centroid + distance to tumor periphery). (C) Local Hi-C contact patterns and corresponding subcompartment for the 24 samples on a representative region of chr2: 130,000,000–190,000,000. Example regions with variable subcompartments were highlighted in gray.

Fig. 2.. Subcompartment analysis.

(A) Barplot of number of each subcompartment in each sample. (B and C) Boxplot of gene expression (B) and ATAC-seq read counts (C) in each subcompartment in P529_6. The center line denotes the median, and the top and bottom of the boxes denote the first and third quartiles, respectively. P values calculated using two-sided Wilcoxon rank-sum test. (D) UMAP clustering based on subcompartments. (E). Barplots of number of each subcompartment at most variable regions (left) and most stable regions (right) in each sample at 100-kb resolutions. (F) Barplot showing number of A/B switched subcompartments between each GBM sample and normal sample. (G) Sankey diagram showing subcompartments changes from NU0820 to P524_9. (H) Heatmap of recurrent A/B-switched regions between GBM samples and normal samples. (I) Boxplots showing the average gene expression in each group for the top three clusters (top boxplot) and the bottom two clusters (bottom boxplot) from the heatmap in (H). The numbers of genes were labeled on top of the boxes. (J) Subcompartment tracks showing examples of B-to-A (top) and A-to-B (bottom) compartment switches between normal samples (group 3) and GBM samples (groups 1 and 2) in (H). The average expression of CHI3L2, CAVIN2, and TMEFF2 for three groups was labeled in black, green, and purple. (K) Heatmap of A/B-switched regions between P530 temporal and frontal region samples. The numbers of A/B-switched subcompartments were labeled on the left side. (L) Boxplot of log₂ fold change (log₂FC) of gene expression in A-to-B, stable, and B-to-A regions between P530 temporal and frontal region samples. Log₂FC is defined as average gene expression in P530 temporal region samples/frontal region samples. The numbers of genes located in A-to-B– and B-to-A–switched regions were labeled at the bottom of the box. P values calculated using two-sided Wilcoxon rank-sum test.

Fig. 3.. GBM-specific chromatin interactions and intratumoral heterogeneity in chromatin interactions.

(A) Barplot showing number of predicted chromatin loops in each sample. (B) Boxplot of expression of genes located in E-P only, P-P only, E-P, and P-P overlap and P-N loop anchors in P524_9. The numbers of genes in each category were labeled on top. P values calculated using two-sided paired Wilcoxon rank-sum test. (C) Boxplot of expression of genes with different numbers of enhancers (1, 2, 3, and ≥4) in P524_9. The numbers of genes in each category were labeled on top. (D) Barplot showing number of GBM sample–specific E-P loops in each sample. (E) Aggregate peak analysis (APA) showing interaction frequency of P530 sample–specific E-P loops at three normal samples and each P530 sample at 10-kb resolution. The numbers of specific E-P loops were labeled on the left side. (F) Boxplot showing expression of genes located in P530_2 (left)– and P530_19 (right)–specific E-P loop anchors in normal sample (average of three normal samples) and P530_2 and P530_19. The number of genes was labeled on top. P values calculated using two-sided paired Wilcoxon rank-sum test. (G) Multiresolution Knight-Ruiz (KR) matrix-balance normalized Hi-C maps of example regions on chr3 of P530_temporal region combined sample. (H) APA showing interaction frequency of P530 temporal and frontal region combined specific loops at three normal samples, P530 temporal and frontal region combined sample, and each individual P530 sample, respectively, at 10-kb resolution. The numbers of specific loops were labeled on the left side. (I) KR balanced Hi-C map at chr11: 43,000,000 to 44,000,000 (left) and chr1: 169,000,000 to 170,000,000 (right) showing P530 temporal region– and frontal region–specific chromatin interactions near API5 and SELL. ATAC-seq and RNA-seq signal in temporal (purple) and frontal region samples (orange) for these regions were shown in the bottom panel.

Fig. 4.. The impact of SVs and CNVs on 3D genome.

(A) Barplot showing number of TADs in each sample. (B) Hi-C data inferred CNV profiles for P475, P498, and P524_9. (C) KR matrix-balanced Hi-C map of chr4: 53,500,000 to 56,500,000 (left) and chr2: 138,500,000 to 144,000,000 (right) in P521 (lower triangle) and NU1650 (upper triangle) from JuiceBox. TADs in normal sample were highlighted by black triangles. CNV tracks for the same region were presented below the Hi-C map. (D) Barplot showing number of SVs whose two breakpoints are in different TADs (Across_TADs) or within the same TADs (Within_TAD) in each GBM sample. (E) Percentage of SVs whose two breakpoints are in different TADs (Across_TADs) or within the same TADs (Within_TAD) in each SV category (deletion, duplication, and inversion) and all three categories combined. SV events in this figure are the combined SV events from all 21 GBM samples. (F) Barplot of number of SV breakpoints located in each subcompartment of normal sample NU0820 in each GBM sample. Numbers of SV breakpoints that were located in A/B compartment of NU0820 were labeled in red/blue. (G) Boxplot of expression of genes close to SV breakpoints and genes not near SV breakpoints in NU0820 (normal) versus P524_4 A-to-B–switched (left) and B-to-A–switched (right) regions. The center line denotes the median, and the top and bottom of the boxes denote the first and third quartiles, respectively. The numbers of genes close to SV breakpoints and genes not near SV breakpoints were labeled on top of the box. P values calculated using two-sided paired Wilcoxon rank-sum test.

Fig. 5.. Trade-off between enhancer hijacking and enhancer amputation.

(A) Phylogenetic trees based on CNVs inferred from Hi-C data. (B and C) Schematic figure showing enhancer hijacking and enhancer amputation due to deletion event (B) and translocation event (C). (D) Hi-C map (top), ATAC-seq, and RNA-seq track (bottom) of chr14: 61,360,000 to 69,780,000 region in P524_9. Deletion marker was highlighted by black circle, and deletion regions were highlighted by black dashed line and gray rectangular. (E) Reconstructed Hi-C map (top), ATAC-seq, and RNA-seq track (bottom) of chr14: 61,360,000 to 69,780,000 region after removing deletion region (chr14: 61,990,000 to 69,110,000) in P524_9. Deletion breakpoints and region were marked by black dashed line. Potential gained enhancers of ZFP36L1 were highlighted in gray rectangular. (F) Hi-C map (top), ATAC-seq, and RNA-seq track (bottom) of chr1: 203,740,000 to 204,600,000 and chr11: 131,340,000 to 132,630,000 region in frontal region sample P530_18. Translocation breakpoints were marked by black dashed line. Potential gained enhancers for ZBED6 were highlighted by gray rectangular. (G) Hi-C map (top), ATAC-seq, and RNA-seq track (bottom) of chr1: 203,200,000 to 204,600,000 region in temporal region sample P530_4. Translocation breakpoint within this region in P530 frontal region samples was marked by black arrow. Potential enhancer loss for ZBED6 in frontal region samples was highlighted in gray rectangular.