Epigenetic profiling reveals key super-enhancer networks driving oncogenesis in HPV-positive HNSCC

Abstract

Human papillomavirus-positive (HPV+) head and neck squamous cell carcinoma (HNSCC) is a growing subset of cancer cases distinct from HPV-negative by fewer genetic mutations and prevalent epigenetic dysregulation. We mapped H3K27ac-marked super-enhancers (SEs) via ChIP-seq in HPV+ patient-derived xenografts (PDXs) and normal oropharyngeal mucosa, identifying tumor-specific SE domains (T-SEDs) enriched for transcription factors (TFs) including TP63, FOSL1, and JUND. These SE-associated TFs regulate key oncogenic pathways and are downregulated by BRD4 inhibition with JQ1, highlighting sensitivity to epigenetic modulation. RNA-seq data revealed coordinated dysregulation of enhancer RNAs and mRNAs near T-SEDs, linked to upregulated pathways including epithelial-mesenchymal transition and E2F targets. JQ1 treatment significantly repressed these tumor-specific pathways, suggesting a therapeutic potential for targeting SE-driven transcription in HPV+ HNSCC. This study underscores the critical role of SEs in epigenetic and transcriptional dysregulation in HPV+ HNSCC, revealing therapeutic targets and providing a framework for future mechanistic studies in this area.

Keywords: Bioinformatics; Cancer; Epigenetics; Molecular network.

Graphical abstract

Figure 1

Study workflow overview (I) H3K27ac-ChIP-seq was performed on two patient-derived xenograft (PDX) tumors and two non-cancerous tissues, analyzed via the LILY algorithm to identify super-enhancers (SE), enhancers (E), and promoters (P). These elements were categorized into domains specific to the tumor (T-SED), normal (N-SED), or common to both (C-SED). (II) TF enrichment in SED, ED, and PD was analyzed using the TF cistrome, complemented by TF expression analysis from primary and cell line mRNA. (III) RNA-seq of 47 tumors and 25 normal samples quantified mRNA and eRNA expressions, with differential analysis conducted using DESeq2. (IV) Functional analysis on two HPV+ HNSCC cell lines, UM-SCC-047 and UPCI-SCC-090. RNA-Seq post-treatment with DMSO or JQ1, followed by DESeq2 for differential expression. (V) Pathway analysis using Hallmark gene sets from MSigDB was based on differentially expressed genes from both primary samples and cell lines. ChIP-seq, chromatin immunoprecipitation sequencing; H3K27ac, histone H3 lysine 27 acetylation; T, tumor; N, normal; PD, promoter domains; ED, enhancer domains; SED, super-enhancer domains; T-SED, tumor-specific super-enhancer domains; N-SED, normal-specific super-enhancer domains; C-SED, common super-enhancer domains; eRNA, enhancer RNA; DESeq2, R package to analyze differential expression.

Figure 2

Transcription factor analysis in HNSCC (A) The transcription factor enrichment analysis was completed using the TF cistrome for PD (left), ED (middle), and SED (right), and, for every TF, the TFBS enrichment was calculated for the tumor (y axes) and normal (x axes) specific domains separately. TFs with statistically significant TFBS enrichment (Bonferroni adjusted p value≤ 0.05) are shown as orange dots in each graph. (B) TFs with statistically significant enrichment in tumors [Δlog(enrichment) > 0] or normals [Δlog(enrichment) < 0] are shown. Their specificity to PD, ED, and/or SED is color-coded to indicate domain types: green (promoter and enhancer), yellow (promoter, enhancer, and super-enhancer), pink (super-enhancer), blue (enhancer), and red (enhancer and super-enhancer). (C) Heatmap of the differential expression of tumor-specific TFs from B. The differential expression is shown as log2-transformed fold change (log2FC) of expression levels of TFs in 47 HPV+ HNSCC tumor vs. 25 normal samples and in two HPV+ cell lines (UM-SCC-047 [047] and UPCI-SCC-090 [090]) treated with JQ1 (treatment) vs. treated with DMSO (control). The heatmap is color-coded according to the log2FC scale. The heatmap is annotated with asterisks (∗, ∗∗, ∗∗∗) to indicate statistical significance (FDR-adjusted p value <0.05, <0.01, and <0.001, respectively). This figure shows TF-binding enrichment statistics only; no chromatin-loop or Hi-C data are depicted.

Figure 3

Cross-model integration identifies a core gene set regulated by tumor-specific super-enhancers and repressed by JQ1 in HPV+ HNSCC. Volcano plots in (A–C) show differential gene expression (log₂ fold change vs. –log₁₀ adjusted p value) for (A) HPV+ tumors vs. normal mucosa (n = 47 vs. 25), (B) UM-SCC-047 cells treated with JQ1 vs. DMSO, and (C) UPCI-SCC-090 cells treated under the same conditions. Red and green dots denote significantly upregulated and downregulated genes (p.adj <0.05), respectively; gray denotes non-significant genes Panels. (D and E) display Venn diagrams showing the intersection between genes that are upregulated in tumors, located near tumor-specific SEs (T-SEDs), and downregulated after JQ1 treatment in (D) UM-SCC-047 and (E) UPCI-SCC-090. This integrative analysis yields a 91-gene core set that is epigenetically activated in tumors and repressed by BET inhibition. (F) shows a heatmap of expression patterns for the 91-gene core set across tumors and both cell lines, revealing consistent upregulation in HPV+ tumors and suppression following JQ1 treatment. (G) Presents gene ontology enrichment analysis (biological process terms) for the core set, with top terms including mitotic nuclear division, chromosome segregation, and G2/M checkpoint regulation. Dot size reflects the number of genes per term; all terms shown have adjusted p-values <0.001.

Figure 4

Distance analysis of super-enhancer domain (SED) activity (A) Differential expression of genes by distance to the nearest SED within a 2-megabase (Mb) range. The mRNA log2 fold change in expression in tumor vs. normal samples is plotted for genes adjacent to T-SED (red) and N-SED (blue) across 100 kilobases (Kb) increments. Density ridgeline plots illustrate the distribution of log2 fold changes, with asterisks indicating statistically significant differences from two-sided Wilcoxon rank-sum tests, with p values Benjamini-Hochberg (BH) adjusted across bins (reported as “p.adj” = FDR) (∗ FDR <0.05, ∗∗ FDR <0.01, ∗∗∗ FDR <0.001) between genes adjacent to T-SED and N-SED within specific distance bins. (B) Number of genes categorized according to their distance from the nearest SED. The number of genes within each 100 Kb distance bin from the nearest SED is shown for T-SED (red) and N-SED (blue). The table inset details the gene counts and BH-adjusted Wilcoxon p values (“p.adj”, i.e., FDR) for significant distance bins, revealing a higher concentration of genes near T-SEDs and significant differences in log2 fold change distributions between T-SED and N-SED. Significant enrichment of T-SED and N-SED is observed across different distance bins, with the strongest enrichment at approximately 0.1 Mb (p.adj: 1.04E-76) and decreasing but remaining significant at distances up to 1.3 Mb (p.adj: 9.82E-03).

Figure 5

Differential expression and correlation of enhancer RNAs (eRNAs) and mRNAs in tumor and normal samples (A) Volcano plot of differential eRNA expression. The x axis represents the log2 fold change of eRNA expression, and the y axis represents the -log10 of FDR-adjusted p values (“p.adj”, Benjamini-Hochberg). eRNAs with significant changes are highlighted, with vertical dashed lines indicating log2 fold change thresholds (±1.5) and the horizontal dashed line representing an FDR (“p.adj”) threshold of 0.05 (plotted as -log10 for visualization). (B) Hexbin plot comparing log2 fold changes of eRNAs and mRNAs. The x axis represents the log2 fold change of mRNA expression, and the y axis represents the log2 fold change of eRNA expression. Points are color-coded by sample source: red for tumor, blue for normal, and black for common SEDs. The color saturation inside the hexagons indicates the count, reflecting the frequency or abundance of the data points. The correlation between eRNA and mRNA expression changes is shown, with Kendall’s Tau value of 0.5856 (p < 0.000001; two-sided, from R cor.test), indicating a statistically significant relationship.

Figure 6

Gene set enrichment analysis of primary tissues and JQ1-treated cells The dot plot illustrates the results of gene set enrichment analysis using MSigDB Hallmark gene sets, comparing differentially expressed genes across patient samples (tumors vs. normal) and cell lines UPCI-SCC-090 and UM-SCC-047 after treatment with JQ1 (JQ1 vs. DMSO). Each row represents a specific Hallmark pathway; dot color encodes the normalized enrichment score (NES; red = positive enrichment, blue = negative), and gray dots indicate no significant enrichment (FDR-adjusted p value [p.adj, Benjamini-Hochberg] > 0.05), while colored dots indicate significant enrichment (p.adj ≤0.05). Color intensity reflects |NES| (stronger color = larger magnitude of enrichment).