S. Carson Callahan1,2
Veena Kochat1,3Zhiyi Liu1,2
Ayush T. Raman1,4,5
Margarita Divenko1Jonathan Schulz1Christopher J. Terranova1
Archit K. Ghosh1
Ming Tang1,6Faye M. Johnson7,8
Jing Wang7,9
Heath D Skinner10,11Curtis R. Pickering2,7
Jeffrey N. Myers2,7*
Kunal Rai1,7*- 1Department of Genomic Medicine, University of Texas MD Anderson Cancer Center, Houston, TX, United States
- 2Department of Head and Neck Surgery, The University of Texas MD Anderson Cancer Center, Houston, TX, United States
- 3Department of Surgical Oncology, University of Texas MD Anderson Cancer Center, Houston, TX, United States
- 4Graduate Program in Quantitative Sciences, Baylor College of Medicine, Houston, TX, United States
- 5Epigenomics Program, Broad Institute of MIT and Harvard, Cambridge, MA, United States
- 6Department of Data Science, Dana-Farber Cancer Institute, Boston, MA, United States
- 7The University of Texas Graduate School of Biomedical Sciences, Houston, TX, United States
- 8Department of Thoracic/Head & Neck Medical Oncology, The University of Texas MD Anderson Cancer Center, Houston, TX, United States
- 9Department of Bioinformatics and Computational Biology, The University of Texas MD Anderson Cancer Center, Houston, TX, United States
- 10Department of Radiation Oncology, The University of Texas MD Anderson Cancer Center, Houston, TX, United States
- 11Department of Radiation Oncology, University of Pittsburgh, UPMC Hillman Cancer Center, Pittsburgh, PA, United States
Head and neck squamous cell carcinoma (HNSCC) is a heterogeneous disease with significant mortality and frequent recurrence. Prior efforts to transcriptionally classify HNSCC into groups of varying prognoses have identified four accepted molecular subtypes of the disease: Atypical (AT), Basal (BA), Classical (CL), and Mesenchymal (MS). Here, we investigate the active enhancer landscapes of these subtypes using representative HNSCC cell lines and identify samples belonging to the AT subtype as having increased enhancer activity compared to the other 3 HNSCC subtypes. Cell lines belonging to the AT subtype are more resistant to enhancer-blocking bromodomain inhibitors (BETi). Examination of nascent transcripts reveals that both AT TCGA tumors and cell lines express higher levels of enhancer RNA (eRNA) transcripts for enhancers controlling BETi resistance pathways, such as lipid metabolism and MAPK signaling. Additionally, investigation of higher-order chromatin structure suggests more enhancer-promoter (E-P) contacts in the AT subtype, including on genes identified in the eRNA analysis. Consistently, known BETi resistance pathways are upregulated upon exposure to these inhibitors. Together, our results identify that the AT subtype of HNSCC is associated with higher enhancer activity, resistance to enhancer blockade, and increased signaling through pathways that could serve as future targets for sensitizing HNSCC to BET inhibition.
Introduction
Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide and the predominant form of head and neck cancer (Marur et al., 2010; Zhou et al., 2016). In the United States, over 60,000 new HNSCC cases and more than 13,000 HNSCC deaths are reported per year (Zhou et al., 2016). HNSCC covers a wide variety of anatomical sites, including the oral cavity, oropharynx, hypopharynx, and larynx (Forastiere et al., 2001). The prognosis for HNSCC is overall poor with a 5-year survival of approximately 50%, which has remained relatively unchanged for decades (Marur et al., 2010). This is largely attributed to factors such as late stage at initial presentation and high rates of primary tumor recurrence (Bonner et al., 2006; Pickering et al., 2013). Treatment for HNSCC involves combinations of surgery, chemotherapy, and radiotherapy, with exact treatment plans depending on tumor location and TNM stage (Marur et al., 2010; Zhou et al., 2016).
To date, studies on HNSCC have focused largely on genomic characterizations such as exome sequencing and copy number alterations. The most common alterations, such as mutations in TP53 at 17p13 and alterations in p16 at 9p21, have been known for decades (Forastiere et al., 2001; Zhou et al., 2016). More recent comprehensive analyses of HNSCC tumors have supported these previous findings, in addition to identifying common alterations in the Notch1 pathway and cell cycle genes (Agrawal et al., 2011; Pickering et al., 2013; Cancer Genome Atlas, 2015). Unfortunately, very few of these studies have resulted in clinically actionable findings. There are, however, some disputed exceptions, such as the EGFR inhibitor cetuximab, which showed benefit when combined with radiotherapy (Bonner et al., 2006).
One interesting result of these and other studies is the notion of molecular subtypes of disease. Inspired by similar studies in other tumors such as breast, lung, and brain, transcriptomic data from patient samples was used to classify head and neck tumors into 4 subtypes: Atypical (AT), Basal (BA), Classical (CL), and Mesenchymal (MS) (Chung et al., 2004; Walter et al., 2013; Cancer Genome Atlas, 2015). These studies have largely focused on the relationship of these subtypes to genomic alterations, such as mutation patterns, copy number changes, or alterations in key transcription factor expression, and clinical features, such as progression free survival and lymph node metastasis at time of diagnosis. However, there have been very few studies describing the epigenomic features of the subtypes. The importance of chromatin modification states in HNSCC is further evidenced by the finding that global levels of certain histone tail modifications correlate with clinical measurements such as tumor stage, cancer-specific survival, and disease-free survival in oral squamous cell carcinoma (Chen et al., 2013). Because there are currently only a sparse number of HNSCC epigenomics datasets, particularly in the realm of histone modifications and chromatin regulation, there remains an unmet need to investigate these aspects of gene regulation and leverage newly discovered biology to better define the disease and develop new therapeutic approaches (Castilho et al., 2017; Serafini et al., 2020).
Since HNSCC subtypes are defined by their transcriptomic signatures, it stands to reason they would also have unique epigenomic features, such as enhancer landscapes, that may, in part, be driving the defining transcriptomic signatures. Through mapping of H3K27ac-marked active enhancers in 28 HPV-negative HNSCC cell lines, we demonstrate that the AT subtype is characterized by high enhancer activity. Consistently, the AT subtype was associated with resistance to enhancer-blocking bromodomain inhibitors (BETi). BETi resistance pathways specifically showed high enhancer activity as measured by nascent enhancer RNAs (eRNAs) and enhancer-promoter contacts, providing mechanistic insights into the aggressive nature of the AT subtype. Overall, our data suggests high enhancer activity as an epigenetic feature of atypical HNSCCs.
Methods
Cell culture
Human HNSCC cell lines were acquired and characterized as previously described (Zhao et al., 2011). Briefly, cell lines were cultured in DMEM supplemented with 10% FBS, L-glutamine, sodium pyruvate, nonessential amino acids, vitamins, and 1% penicillin-streptomycin. All cell lines were cultured at 37°C in an atmosphere of 5% CO2.
RNA-Sequencing processing and analysis
RNA-seq data for cell line subtype assignments were obtained as raw counts, processed as previously described (data available at GEO accession GSE122512) (Kalu et al., 2017; Gleber-Netto et al., 2019). To assign HNSCC cell lines to their representative subtypes, we used the HNSCC gene list templates generated by Yu et al. (Yu et al., 2019) and utilized the CMScaller workflow and implementation of the NTP algorithm to find the closest matching subtype for each cell line (FDR <0.1) based on their transcriptomic profiles (Hoshida, 2010; Eide et al., 2017). Upregulated genes for each subtype were computed using CMScaller in a one-vs-rest fashion.
For the PLX51107 treatment RNA-seq experiments, representative cell lines were selected for the AT subtype (HN4) and a non-AT subtype (MDA1186, CL subtype) and treated with DMSO, GR50 of the MDA1186, or GR50 of HN4 for 24 h prior to RNA isolation. RNA extraction was performed using an RNeasy Mini Kit per manufacturer’s instructions (Qiagen). Isolation of mRNA was performed using NEBNext Poly(A) mRNA Magnetic Isolation Module and libraries were prepared using the NEBNext Ultra II Directional RNA Library Prep kit (New England BioLabs). Library quality was checked on an Agilent TapeStation 4150 and quantified by Qubit 2000 fluorometer (Invitrogen). Libraries were pooled in equimolar ratios and sequenced on Illumina NovaSeq6000 SP runs with paired-end 100-bp reads at The Advanced Technology Genomics Core (ATGC) at MD Anderson Cancer Center.
PLX51107 treatment RNA-seq raw reads were processed using the provided pipeline: https://github.com/sccallahan/QUACKERS_RNAseq-pipeline. In brief, raw reads were aligned to the hg19 genome using STAR v2.7.2b (Dobin et al., 2013) and quality checked using FastQC v0.11.8 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Counts were generated using featureCounts from subread v1.6.3 80 (Liao et al., 2013). Downstream normalization and differential expression analysis were performed using DESeq2, with size factors being calculated using data-driven housekeeping gene method as implemented in the CustomSelection R package (Love et al., 2014; Dos Santos et al., 2020). Pathway enrichment analyses were performed using GSEA’s pre-ranked list option (Subramanian et al., 2005). Overlaps of HN4 and MDA1186 low dose PLX51107 differentially expressed genes were performed using the VennDiagram package in R, and the HN4 uniquely upregulated gene list was subjected to pathway enrichment analysis using the gsea-msigdb online tool (http://www.gsea-msigdb.org/gsea/msigdb/annotate.jsp).
CCLE RNA-seq data were downloaded as raw counts from the DepMap download portal (https://depmap.org/portal/download/). Subtype assignment and downstream analysis were carried out as above.
Whole exome sequencing processing and analysis
Whole exome sequencing (WES) data was processed as previously described and obtained as a MAF file from the authors (Kalu et al., 2017; Gleber-Netto et al., 2019). To cluster the cell lines based on mutation background, all mutation calls were binarized to 1 or 0 to represent “mutated” or “not mutated,” respectively. The Jaccard distance matrix was then computed, and the resulting matrix was clustered using Ward’s minimum variance method. Total mutational burden was calculated by summing the number of mutations per sample, then grouping the samples based on their assigned molecular subtype. Data for cell line tissue of origin and “source” were obtained as previously described (Zhao et al., 2011).
ChIP-Sequencing processing and analysis
ChIP assays were performed as described previously (Terranova et al., 2018). Briefly, approximately 2 × 107 cells were harvested by scraping. Samples were cross-linked with 1% (wt/vol) formaldehyde for 10 min at 37°C with shaking. After quenching with 150 mM glycine for 5 min at 37°C with shaking, cells were washed twice with ice-cold PBS and frozen at −80°C for further processing. Cross-linked pellets were thawed and lysed on ice for 30 min in ChIP harvest buffer (12 mM Tris-Cl, 1 × PBS, 6 mM EDTA, 0.5% SDS) with protease inhibitors (Sigma). Lysed cells were sonicated with a Bioruptor (Diagenode) to obtain chromatin fragments (∼200–500 bp) and centrifuged at 15,000 × g for 15 min to obtain a soluble chromatin fraction. In parallel with cellular lysis and sonication, antibodies (5 μg/3 × 106 cells) were coupled with 30 μL of magnetic protein G beads in binding/blocking buffer (PBS +0.1% Tween +0.2% BSA) for 2 h at 4°C with rotation. The antibody used for ChIP was anti-H3K27ac (Abcam; ab4729). Soluble chromatin was diluted five times using ChIP dilution buffer (10 mM Tris-Cl, 140 mM NaCl, 0.1% DOC, 1% Triton X, 1 mM EDTA) with protease inhibitors and added to the antibody-coupled beads with rotation at 4°C overnight. After washing, samples were treated with elution buffer (10 mM Tris-Cl, pH 8.0, 5 mM EDTA, 300 mM NaCl, 0.5% SDS), RNase A, and Proteinase K, and cross-links were reversed overnight at 37°C. Immune complexes were then washed five times with cold RIPA buffer (10 mM Tris–HCl, pH 8.0, 1 mM EDTA, pH 8.0, 140 mM NaCl, 1% Triton X-100, 0.1% SDS, 0.1% DOC), twice with cold high-salt RIPA buffer (10 mM Tris–HCl, pH 8.0, 1 mM EDTA, pH 8.0, 500 mM NaCl, 1% Triton X-100, 0.1% SDS, 0.1% DOC), and twice with cold LiCl buffer (10 mM Tris–HCl, pH 8.0, 1 mM EDTA, pH 8.0, 250 mM LiCl, 0.5% NP-40, 0.5% DOC). ChIP DNA was purified using SPRI beads (Beckman Coulter) and quantified using the Qubit 2000 (Invitrogen) and TapeStation 4150 (Agilent). Libraries for Illumina sequencing were generated following the New England BioLabs (NEB) Next Ultra DNA Library Prep Kit protocol. Amplified ChIP DNA was purified using double-sided SPRI to retain fragments (∼200–500 bp) and quantified using the Qubit 2000 and TapeStation 4150 before multiplexing.
Raw fastq reads for all ChIP-seq experiments were processed using a Snakemake based pipeline https://github.com/crazyhottommy/pyflow-ChIPseq. Briefly, raw reads were first processed using FastQC and uniquely mapped reads were aligned to the hg19 reference genome using Bowtie version 1.1.2 (Langmead et al., 2009). Duplicate reads were removed using SAMBLASTER (Faust and Hall, 2014) before compression to bam files. To directly compare ChIP-seq samples, uniquely mapped reads for each mark were downsampled per condition to 15 million, sorted, and indexed using samtools version 1.5 (Li et al., 2009). To visualize ChIP-seq libraries on the IGV genome browser, we used deepTools version 2.4.0 (Ramirez et al., 2016) to generate bigWig files by scaling the bam files to reads per kilobase per million (RPKM) and WiggleTools (Zerbino et al., 2014) to create average profile plots for each molecular subtype.
Peak overlaps were performed by first generating consensus peak files for each subtype, defined as any peak found in at least 2 samples from the subtype. The resulting 4 peaksets (one per subtype) were then used as input for intervene’s (Khan and Mathelier, 2017) upset module to generate upset plots and common/unique peaksets for further analysis. Gene linkage was performed using previously published enhancer-promoter linkage data from Cao et al. (Cao et al., 2017), and the resulting gene list was used as input for pathway enrichment analysis using the gsea-msigdb online tool. Enrichment plots were generated using two definitions of common peaks. The first method uses DiffBind (Stark and Brown, 2020) (https://bioconductor.org/packages/release/bioc/html/DiffBind.html) to define a peakset using any peak found in at least 2 samples, irrespective of subtype. The second uses the peaks found in all subtypes when overlapped using the intervene package mentioned previously. Both resulting peaksets were used as input for ngs.plot (Shen et al., 2014) to generate figures.
Drug response assays
Cell confluence and proliferation were measured using the IncuCyte ZOOM system (Essen Biosciences). For each cell line, seeding density was optimized such that the cells would be in their exponential growth phase for the duration of drug treatment. On day 0, cells were seeded into 96 well plates and left in the incubator overnight. On day 1, media containing either drug (PLX51107 or OTX015) or DMSO was added to the wells. Plates were then left in the IncuCyte ZOOM with treated media for 72 h, at which point cell confluence was measured. Each assay was performed in biological duplicate with technical triplicate wells. Drug response metrics were calculated using GR metrics (Hafner et al., 2016) using cell confluence as a proxy for growth rate, and GRAOC values from each cell line were combined based on their molecular subtype for statistical analysis.
CCLE drug response data were downloaded and processed using the PharmacoGx “auc_recomputed” dataset (Smirnov et al., 2016). Compounds which were missing data for more than 25% of samples were excluded. For the JQ1 analysis, 1 sample was missing data and was imputed using predictive mean matching [as implemented in the mice package (Buuren and Groothuis-Oudshoorn, 2011) (https://www.jstatsoft.org/article/view/v045i03)] on the complete, filtered drug response matrix.
PRO-Seq processing and eRNA analysis
Extraction of nuclei and precision run-on reaction was carried out as described previously (Mahat et al., 2016). Nuclei were isolated from approximately 10 million cells after treating with 12 ml of ice-cold swelling buffer (10 mM Tris-HCl pH 7.5, 2 mM MgCl2, 3 mM CaCl2) for 10 min and scraping out the cells. After spinning at 600 × g for 10 min at 4°C, the supernatant was removed and the cells were lysed in 10 ml of lysis buffer (10 mM Tris-HCl pH 7.5, 2 mM MgCl2, 3 mM CaCl2, 10% glycerol, 0.5% NP-40, 4 U/ml SUPERase inhibitor) on ice for 5 min. The lysate was spun at 600 × g for 8 min and the nuclei were collected. The nuclei were then resuspended in 1 ml of freezing buffer (50 mM Tris-HCl pH 8.0, 5 mM MgCl2, 0.1 mM EDTA, 40% glycerol) and spun at 900 × g for 10 min. For performing precision nuclear run-on reaction, nuclei were resuspended in 100 µL of freezing buffer and added to 100 µL of NRO-reaction mix - NRO-reaction buffer (10 mM Tris-HCl pH 8.0, 5 mM MgCl2, 300 mM KCl), 1 mM DTT, 100 U/ml SUPERase-In, 1% Sarkosyl, 250 M ATP, 250 M GTP, 50 M biotin-11-UTP, 50 M biotin-11-CTP. Reaction was carried out at 29°C for 4 min. RNA was extracted using TRIzol. Base hydrolysis was carried out by heat denaturing briefly at 65°C for 40 s following by cooling on ice and treatment with 1N NaOH for 6 min on ice. The sample with fragmented RNA was neutralized with 1 M Tris-HCl pH 6.8 and isolated by passing through P-30 column (Biorad, #732-6250). The NRO-reaction products containing biotinylated RNA was purified using Streptavidin C1 beads which were washed thrice with wash buffer (0.1 N NaOH, 50 mM NaCl) and twice with 100 mM NaCl. The washed beads were resuspended in binding buffer (10 mM Tris-HCl pH 7.4, 150 mM NaCl, 0.1% Triton X-100) and added to the sample and incubated at room temperature for 30 min on a rotator. After removing the supernatant using a magnetic stand, beads were washed twice with high salt wash buffer (50 mM Tris-HCl pH 7.4, 2 M NaCl, 0.5% Triton X-100), once with low salt wash buffer (10 mM Tris-HCl pH 7.4, 300 mM NaCl, 0.1% Triton X-100) and twice with no salt wash buffer (5 mM Tris-HCl pH 7.4, 0.1% Triton X-100). Beads were then resuspended in TRIzol and RNA was extracted. The bead purification of biotinylated RNA was performed once more, and RNA was extracted using TRIzol to improve the purity of the sample.
Libraries were generated based on previously described protocol (Van Nostrand et al., 2016). Isolated RNA samples were dephosphorylated by FastAP (ThermoFisher) and T4 Polynucleotide kinase (NEB). Samples were cleaned up using MyONE Silane beads and RNA was isolated with RLT buffer (Qiagen). To the eluted RNA, a barcoded RNA adapter (RiL19) was ligated to the 3′ end using T4 RNA ligase (NEB). The 3′ adaptor ligated RNA was again cleaned up as mentioned above. RNA was then reverse transcribed with AR17 primer and AffinityScript reverse transcriptase (Agilent). cDNA was then cleaned up by treating with ExoSAP-IT (Affymetrix) to remove excess oligonucleotides. Excess RNA was removed from cDNA by treating with 1 M NaOH at 70°C for 12 min and neutralizing with 1 M HCl. cDNA was then cleaned up with MyONE Silane beads and RLT buffer and eluted in 5 mM Tris-Cl, pH 7.5. A second 5′adaptor (rand3Tr3) was ligated to cDNA with T4 RNA ligase in an overnight reaction at room temperature. The adaptor ligated cDNA was then cleaned up with MyONE Silane beads and RLT buffer and eluted in 10 mM Tris-Cl, pH 7.5. cDNA samples were then PCR amplified using NEBNext® Ultra™ II Q5® Master Mix multiplexing was done with D50X and D70X primers. Libraries were size selected and purified using SPRI beads. Final libraries were quantified using D1000 tapestation and Qubit™ dsDNA HS Assay Kit (Thermo Fisher Scientific) and sequenced using NovaSeq6000 with 100 nt paired-end format.
Fastq files from precision nuclear run-on sequencing (PRO-seq) experiments were processed using the previously described PEPPRO pipeline (Smith et al., 2021). Briefly, fastq files first undergo pre-processing steps of adapter removal, read deduplication, read trimming, and reverse complementation. The resulting files are then “pre-aligned” to the human rDNA genome to siphon off these unwanted reads. The rDNA-removed files are then aligned to the human hg19 genome using bowtie2 (Langmead and Salzberg, 2012). After quality control assessment, 5 samples from the AT subtype (representing the 3 unique cell lines used in the drug response assays) and 3 samples from the CL subtype (representing the 2 unique cell lines used in the drug response assays) were carried forward for further analysis. The aligned, sorted bam files for these samples were used as input for downstream analysis using the previously described NRSA downstream analysis pipeline (Wang et al., 2018). In brief, NRSA uses bidirectional transcription in intergenic regions to identify and call enhancers/eRNAs. A raw counts table for these eRNAs is then generated and fed into the DESeq2 tool for differential expression analysis. Identified enhancers are assigned to their nearest genes to generate a list of genes with upregulated eRNA expression in the AT subtype, which was then used as input for pathway enrichment analysis.
TCGA RPKM expression levels of numerous eRNAs from typical enhancers were downloaded from publicly available datasets (https://bioinformatics.mdanderson.org/Supplements/Super_Enhancer/TCEA_website/parts/3_eRNA_quantification.html) based on previously published work (Chen et al., 2018; Chen and Liang, 2020). TCGA mRNA-seq for subtype assignment was downloaded from FireBrowse (http://firebrowse.org/). TCGA samples were assigned to molecular subtypes by first generating templates based on previously assigned molecular subtypes from the initial HNSCC TCGA cohort (Cancer Genome Atlas, 2015). These assigned subtypes were then expanded to the current cohort of samples by using the CMSCaller functionality described above. As before, samples were only retained for further analysis if they possessed an assignment FDR <0.1. Samples were then grouped into “Atypical” or “Other” based on their molecular subtype, and significant differential expression of eRNAs was determined by > 1.5 fold-change in expression and FDR <0.05. As with the PRO-seq data, these eRNAs were linked to their nearest gene using bedtools via the bedr R package (Quinlan and Hall, 2010), and the resulting gene list was used as input for pathway enrichment analysis. Intersections of the TCGA eRNA enriched pathways and PRO-seq eRNA enriched pathways were performed using the Venn Diagram package in R and significance was calculated using the hypergeometric overlap method, with a Universe size set to the number of unique pathways in a particular gene set. Only pathways with a p.adjust <0.25 with a maximum of 20 enriched pathways per gene set were included in the analysis.
HiChIP protocol, processing, and analysis
HiChIP was performed as described (Mumbach et al., 2016). Briefly, 1 × 107 cells for each HNSCC cell line (1 unique cell line per HNSCC subtype) were crosslinked. In situ contacts were generated in isolated and pelleted nuclei by DNA digestion with MboI restriction enzyme, followed by biotinylation of digested DNA fragments with biotin–dATP, dCTP, dGTP, and dTTP. DNA was then sheared with Bioruptor (Diagenode); chromatin immunoprecipitation was done for H3K27Ac with use of anti-H3K27ac antibody. After reverse-crosslinking, 150 ng of eluted DNA was taken for biotin capture with Streptavidin C1 beads followed by transposition with Tn5. In addition, transposed DNA was used for library preparation with Nextera Ad1_noMX, Nextera Ad2.X primers, and Phusion HF 2XPCR master mix. The following PCR program was performed: 72°C for 5 min, 98°C for 1 min, then 11 cycles at 98°C for 15 s, 63°C for 30 s, and 72°C for 1 min. Afterward, libraries were two-sided size selected with AMPure XP beads. Libraries were paired-end sequenced with reading lengths of 76 nucleotides.
Using HiC-Pro (Servant et al., 2015), HiChIP paired-end reads were aligned to the hg19 genome with duplicate reads removed, assigned to MboI restriction fragments, and filtered for valid interactions. Interaction matrices were then generated with the same software. To generate anchor points for downstream looping analysis, outputs from HiC-Pro were used as inputs for peak calling in HiChIP-Peaks (Shi et al., 2020). To ensure loops were called from similar background enhancers, peaks from HiChIP-Peaks were concatenated into a single file and used as anchor point inputs for loop calling via hichipper (Lareau and Aryee, 2018). HiChIP loop visualization was performed using DNAlandscapeR (https://molpath.shinyapps.io/dnalandscaper/).
Statistical analysis
Statistical analyses, including generation of graphs and plots, were performed using R versions 3.4.4 and 3.6.0. Significance levels are * = p < 0.05, ** = p < 0.01, and ***p < 0.005 unless otherwise indicated in figure legends. Statistical tests utilized are as indicated in respective text and figure legends.
Results
Head and neck squamous cell carcinoma cell lines assigned to known molecular subtypes have similar mutation profiles and tissue origins
Identifying inter-tumor heterogeneity can help better understand the diversity of biological mechanisms driving the neoplastic phenotypes within pathology-based tumor types (e.g., breast cancer, colon cancer, etc.) and discover targeted therapies for specific patient populations (Guinney et al., 2015; Fragomeni et al., 2018; Collisson et al., 2019; Vasaikar et al., 2019). We sought to define heterogeneity within HNSCC patients, especially at the epigenetic level. To this end, we first leveraged published work by Yu et al., which extensively studied CCLE cancer cell lines and their appropriateness as models of human cancer by comparing them to corresponding TCGA tumors (Yu et al., 2019). As part of this work, the group generated “templates” of gene expression values for numerous subtypes in 9 different tumor types. Using the HNSCC templates from this study (one per molecular subtype) and RNA-seq data from the panel of cell lines available to us (Supplementary Table S1), we utilized the nearest template prediction method, as implemented in the CMScaller R package, to assign our cell lines to their most representative molecular subtype (Figure 1A) (Hoshida, 2010; Eide et al., 2017). After selecting only samples with an assignment FDR <0.1, twenty-eight HPV-negative HNSCC cell lines were successfully matched to a molecular subtype, resulting in 7 AT samples, 9 BA samples, 5 CL samples, and 7 MS samples (Figure 1B). Analysis of RNA-seq data in one-versus-rest comparisons demonstrated varying levels of differential expression based on subtype, with a large number of upregulated genes (FC > 3, n = 756) in the AT subtype (Figure 1C).
FIGURE 1. Cell line subtype assignments and characteristics. (A) Schematic of workflow used to assign HNSCC cell lines to their respective subtypes using RNA-seq data. (B) Table of subtype assignments for each of the 28 cell lines used in this study. (C) Heatmap of gene expression modules in each molecular subtype, defined as FC > 3 in a one-vs-rest comparison. (D) Hierarchical clustering of the 28 cell lines based on Jaccard distance metrics obtained from binarized mutation counts from WES data. (E) Boxplots demonstrating total number of mutations in each sample, grouped by molecular subtype (p = NS for each comparison). (F) Stacked barplot showing distribution of cell line anatomic location for each molecular subtype. (G) Pie chart showing percentage of samples in each molecular subtype that came from primary, recurrent, or metastatic lesions.
To further investigate a potential genomic basis that could be driving the transcriptomic partitioning into molecular subtypes, we investigated WES data on our panel of cell lines. We first clustered our lines based on their mutational background, and, interestingly, we did not observe any clustering of samples from the same molecular subtype. In fact, the only examples of tight clusters in the data came from 3 matched pairs of cell lines in which the samples were either from a primary or metastatic lesion of the same patient (Figure 1D). Similarly, we did not observe any significant differences between subtypes based on total mutational burden (Figure 1E). This finding is largely consistent with HPV-negative TCGA data, with the singular exception of the AT vs. MS comparison showing a significant difference in mutation number in TCGA (p.adj = 0.031) (Supplementary Figure S1) (Cancer Genome Atlas, 2015). Importantly, observed clustering was neither associated with the anatomic site of origin of the primary tumor from which each the cell lines were derived, nor with the type of tumor the sample was from (e.g., primary vs. recurrence). With the possible exception of the MS subtype, all of the HNSCC subtypes had a fairly equal distribution of samples from the oral cavity, oropharynx, and larynx (Figure 1F). We note this varies from the findings in the TCGA data, where subtype was correlated with anatomic location (Cancer Genome Atlas, 2015). With respect to cell line source, the AT subtype was the only one to contain cell lines from all 3 groups of samples (primary, recurrence, and metastasis), while BA contained only primary and recurrence, and CL and MS contained only primary and metastasis (Figure 1G). Taken together, these results demonstrate that HNSCC molecular subtypes can be successfully assigned to cell lines using RNA-seq data, and that despite their transcriptomic differences, the unique HNSCC subtypes do not have significantly different mutational backgrounds, overall mutation burden, or tissues of origin from one another in our cell line models.
Head and neck squamous cell carcinoma molecular subtypes are associated with distinct enhancer landscapes
Transcription of a gene is regulated by concerted action of multiple complexes on specific epigenetic elements located in cis or trans the gene promoter. Enhancers are a major component of the gene regulation circuits and known to be deregulated in cancers (Lee and Young, 2013; Herz et al., 2014). They act as binding platforms for transcription factors that, upon various environmental cues relayed by the cell surface signaling pathways, cooperate with chromatin modifying and remodeling machinery to activate target genes (Lee and Young, 2013). We hypothesized that these transcriptional subtypes could have underlying differences in gene regulatory landscapes that could partly explain the observed transcriptomic differences. To investigate these differences, we generated enhancer profiles for each cell line by performing ChIP-seq for the H3K27ac histone mark, which is widely used as a marker of active enhancers (Creyghton et al., 2010; Rada-Iglesias et al., 2011). We next generated consensus peak sets for each subtype by overlapping the enhancer regions of all cell lines within a subtype and taking the set of enhancers that occurred in 2 or more samples of the subtype. This resulted in 4 total consensus peak sets, each representing a unique subtype (i.e., one consensus peak set per subtype). We observed distinct enhancer peak enrichment among the four molecular subtypes (Figure 2A) (Lex et al., 2014; Khan and Mathelier, 2017). Notably, we discovered the AT subtype has a much larger set of consensus enhancers than any of the other subtypes, and the most common subset of enhancers in our analysis is the set that is unique to the AT subtype (Figure 2A). To investigate the function of the 5,683 enhancers unique to the AT subtype, we utilized data from Cao et al. (Cao et al., 2017), which constructed enhancer-target networks across multiple cancer and sample types, to assign each of these enhancers to their target genes (Figures 2B,C). These genes were then used for pathway enrichment analysis, which revealed enrichment for pathways involved in lipid metabolism, MYC signaling, and MAPK signaling (Figures 2B,C, Supplementary Table S2).
FIGURE 2. The Atypical subtype is associated with unique enhancer peaks regulating genes related to lipid metabolism and MAPK signaling. (A) UpSet plot showing the total number of H3K27ac typical enhancer peaks in each molecular subtype (pink horizontal barplot), as well as the number of peaks in each possible intersection of peaksets (black bars and dot plot). (B,C) Visualization of mean bigWig signal for each subtype at (B) MAP3K8 and (C) IGFBP3 enhancer loci containing H3K27ac peaks unique to the AT subtype (green bar/grey shading). (D,E) H3K27ac ChIP-seq enrichment plots of enhancer loci common to all HNSCC molecular subtypes, defined as (D) any peak contained within 2 or more individual samples or (E) the 3,404 peaks shared among all consensus peaksets in (A), demonstrating the strongest signal in the AT subtype. (F) UpSet plot showing the total number of super enhancer peaks in each molecular subtype (blue horizontal barplot), as well as the number of super enhancer peaks in each possible intersection of peaksets (black bars and dot plot). (G) Visualization of mean bigWig signal for each subtype at a MAP3K12 super enhancer containing peaks unique to the AT subtype (green bar/grey shading).
In addition to identifying unique enhancers, we investigated the total H3K27ac signal enrichment across enhancers shared by all 4 subtypes to determine if, in addition to the largest number of H3K27ac peaks, the AT subtype also had greater signal enrichment overall. Indeed, using two separate methods to arrive at a “shared” H3K27ac peak set, we observed that the AT subtype had more enrichment of H3K27ac signal across enhancers shared among all HNSCC subtypes (Figures 2D,E). In agreement with our typical enhancer analysis, we found that the AT subtype also harbored the largest number of called super-enhancers (Figures 2F,G). Further, linking of these super-enhancers to their gene targets not only displayed enrichment for MAPK signaling and lipid metabolism, but also identified enrichment for PI3K and WNT-β-catenin signaling (Supplementary Table S3). These results demonstrate that the AT subtype of HNSCC is enriched for H3K27ac-marked typical enhancers and super-enhancers compared to other HNSCC, and these enhancer regions may activate important cell signaling pathways that are associated with aggressive HNSCCs.
The Atypical subtype is more resistant to bromodomain inhibition
Recently, the realm of “epigenetic” therapies for cancer has been of major interest for both research and in clinical applications (Castilho et al., 2017; Cheng et al., 2019; Bates, 2020). Targeting epigenetic modifications and the proteins that regulate their placement and/or removal is a particularly attractive approach to cancer therapeutics since these modifications are generally considered to be reversible, particularly when compared to more “permanent” changes such as mutations and copy number alterations. One class of compounds with numerous clinical trials for a variety tumor types is BRD and extraterminal domain (BET) inhibitors, which function by inhibiting the “reader” proteins responsible for recognizing and propagating the signal of acetylated histone residues. These inhibitors have been used in prior studies as enhancer-blocking agents, and the pathways we found to be activated by AT-specific enhancers and super-enhancers (e.g., MAPK and PI3K signaling) are well-characterized mechanisms of BET inhibitor resistance (Supplementary Tables S2, S3) (Rathert et al., 2015; Kurimchak et al., 2016; Iniguez et al., 2018; Cochran et al., 2019; Loganathan et al., 2019; Tonini et al., 2020; Yan et al., 2020). Hence, we hypothesized that the AT subtype may be differentially responsive to BET inhibition.
We first investigated HNSCC cell line response to JQ1 using the publicly available CCLE drug response data (Basu et al., 2013; Seashore-Ludlow et al., 2015; Rees et al., 2016). We used the available CCLE RNA-seq data to assign samples to their respective HNSCC molecular subtype, then compared their response to BET inhibition in this dataset. Interestingly, we found that the AT samples have a lower JQ1 AOC (Area Over the Curve, where lower values indicate resistance) than those in the non-AT group (p = 0.0503, Welch’s t-test) (Figure 3A, Supplementary Figure S2A), indicating the AT subtype is more resistant to JQ1 treatment. To extend this analysis, we selected 2 compounds currently being evaluated in clinical trials, OTX015 (Birabresib) and PLX51107, and performed drug response assays in our HNSCC cell lines (Bates, 2020). For each compound, we selected representative cell lines for each molecular subtype (3 AT, 3 BA, 2 CL, and 3 MS), treated with the respective compound for 72 h, then computed the GRAOC of each cell line using cell confluence as a proxy for cell number. We elected to use the GRAOC metric for drug response since GR metrics have been demonstrated to be more reproducible than traditional metrics, such as Area Under the Curve (AUC) and IC50, when measuring drug sensitivity in cancer cell lines (Hafner et al., 2016). As we anticipated based on our previous analysis, the BET inhibitor PLX51107 demonstrated significantly lower GRAOC values in the AT subtype compared to any other HNSCC subtype, indicating an increased resistance to treatment in that group (Figure 3B). The inhibitor OTX015 also displayed a similar trend towards increased resistance in the AT subtype that was similar to, but more pronounced than, the JQ1 response data in the CCLE database (Supplementary Figures S2A,B). Given that BRD proteins are responsible for mediating gene transcription and are the main targets of BET inhibitors, we investigated the expression of this family of proteins both in our cell line panel and in the HNSCC TCGA data (Rathert et al., 2015; Kurimchak et al., 2016; Iniguez et al., 2018; Cheng et al., 2019; Cochran et al., 2019). This analysis demonstrated a significant elevation in BRD7 in the AT subtype compared to other subtypes in our HNSCC cell lines (Supplementary Figures S2C,D), and a significant elevation in BRD4 expression in AT vs. BA and AT vs. MS comparisons in the HNSCC TCGA data (Supplementary Figures S2E,F).
FIGURE 3. Atypical HNSCC shows increased resistance to BET inhibition and uniquely upregulates genes associated with resistance pathways upon treatment. (A) Atypical samples in the HNSCC CCLE dataset demonstrate lower JQ1 AOC values than non-Atypical samples (p = 0.0503, Welch’s t-test). (B) Drug response assays with the BET inhibitor PLX51107 demonstrate the Atypical subtype is significantly more resistant to BET inhibition than other molecular subtypes (*** adj.p < 0.001, ** adj.p < 0.01). (C) Hierarchical clustering of all genes from HN4 and MDA1186 samples treated with DMSO, PLX51107 at GR50 MDA1186 (low), or GR50 HN4 (high). (D) PCA plot of samples as described in (C), displaying separation on the basis of cell line (PC1) and treatment status (PC2). (E) Overlap of genes upregulated (|log2fold-change| > 1.5 & FDR <0.05) in HN4 and MDA1186 at PLX51107 low concentration; numbers in Venn diagram represent size of set. (F) Horizontal barplots of Hallmark (left) and KEGG (right) pathway enrichment results from the 1,437 genes uniquely upregulated by HN4 in (E); pathways highlighted in blue are associated with known mechanisms of BET inhibitor resistance.
To better understand the mechanisms behind this observed resistance to treatment, we selected one representative cell line from the AT subtype (HN4) and one cell line from the non-AT subtypes (MDA1186, CL subtype), treated with PLX51107 or DMSO, and performed gene expression analysis using mRNA-seq profiling in each condition. MDA1186 was treated with PLX51107 at its own GR50 value (hereafter referred to as “low” concentration), and HN4 was treated at its own GR50 value (hereafter referred to as “high”), as well as the low concentration. To ensure PLX51107 behaved similarly to other published BET inhibitors, we created BET inhibitor response signatures using publicly available data and, using GSEA, confirmed that the response of the AT and non-AT cell lines to PLX51107 was consistent with previously documented BET inhibitor response signatures (Supplementary Figures S3C,D) (Puissant et al., 2013; Picaud et al., 2016). Hierarchical clustering of all genes in the RNA-seq dataset demonstrated 2 major clusters, one per cell line, as well as 2 sub-clusters, either DMSO or PLX51107-treated, per major cluster (Figure 3C). To further examine the differences in response to drug treatment, we performed a principal component analysis (PCA), which revealed a first component driven by the cell line identity, and a second component driven by treatment status (Figure 3D). Examining the results from the hierarchical clustering and PCA analysis together, we noted that the majority of transcriptional response to BET inhibition in the AT cell line occurs at the lower drug concentration, and only a minority of gene expression changes occur between the low and high concentrations (Figures 3C,D). Closer inspection of the gene expression heatmap indicates that the genes that are specifically upregulated in the AT subtype after PLX51107 treatment, but not in the non-AT subtypes, may be responsible for mediating resistance to BET inhibition (Figure 3C).
To investigate these uniquely upregulated genes, we overlapped the set of genes upregulated by the AT subtype and by the non-AT subtype at the low PLX51107 concentration. As we expected, we discovered a set of 1437 genes uniquely upregulated in the AT subtype after BET inhibition (Figures 3C,E). Pathway analysis of these 1437 genes using the Hallmarks and KEGG gene sets reveals enrichment for multiple pathways previously demonstrated to convey resistance to BET inhibition and identified by our previous enhancer-based analysis, including MAPK signaling, WNT-β-catenin signaling, phosphatidylinositol signaling, and lipid metabolism pathways (Figure 3F) (Rathert et al., 2015; Kurimchak et al., 2016; Iniguez et al., 2018; Cochran et al., 2019; Loganathan et al., 2019; Tonini et al., 2020; Yan et al., 2020).
These results support our previous hypothesis that the AT subtype is more resistant to BET inhibition than other HNSCC subtypes and suggest the enrichment of H3K27ac-marked enhancers involved in these pathways is a contributing factor to the observed resistance.
Bromodomain inhibitor resistance is mediated by baseline enhancer activity and chromatin structure
Our previous analyses have indicated that the AT subtype of HNSCC has two intriguing properties with respect to BET inhibition: first, H3K27ac-marked enhancers unique to the AT subtype regulate genes enriched for known BET inhibitor resistance pathways, and, second, the AT subtype is able to uniquely upregulate genes enriched for BET inhibitor resistance pathways after treatment with BET inhibitor (Figures 2B,C,G; Supplementary Tables S2, S3; Figures 3E,F). Because of these observations, we suspected the AT subtype may have a stronger baseline enhancer activity at genes involved in resistance pathways and that these genes may have higher numbers of enhancers-promoter contacts, enabling a more robust response to BET inhibitor treatment.
To investigate the activity of enhancers involved in regulating baseline resistance gene expression, we performed PRO-seq to investigate the eRNA landscape of the AT and non-AT subtypes (Supplementary Table S4) (Mahat et al., 2016). Enhancer RNAs are a recently discovered class of non-coding transcripts found at active enhancers that arise from the transcription of enhancer elements themselves and are involved in functions such as regulating gene transcription and controlling enhancer-promoter looping (Arnold et al., 2019; Zhang et al., 2019; Sartorelli and Lauberth, 2020). We used PRO-seq, with a particular focus on eRNAs, to investigate differential enhancer activity in our AT subtype. For this experiment, we expanded our AT group to include the 3 cell lines from our drug assay, and we expanded the non-AT group to include the 2 cell lines from the CL subtype used in our drug assay. Differential expression analysis of PRO-seq-defined eRNAs revealed 321 differentially expressed eRNAs, with 207 being upregulated and 114 being downregulated (Figure 4A).
FIGURE 4. Enhancers of MAPK signaling, WNT signaling, and Cholesterol Homeostasis genes display increased eRNA transcription and enhancer-promoter looping in Atypical HNSCC. (A) Differential transcription (|log2fold-change| > 1.5 & FDR <0.1) of eRNAs between the Atypical and Classical subtypes as measured by PRO-seq (green dots meet fold-change and FDR thresholds, purple dots meet fold-change threshold only). (B) Hallmark pathway enrichment analysis of genes linked to eRNAs with significantly increased transcription from (A); pathways in blue have been previously associated with BET inhibitor resistance. (C) Overlap of hallmark (left) and KEGG (right) pathway enrichment results between PRO-seq-determined significantly enriched eRNAs from (A) and (B) and TCGA-measured differentially expressed eRNAs between Atypical and non-Atypical samples; p values represent hypergeometric tests of gene set enrichment result overlaps; bolded terms represent shared pathways associated with BET inhibitor resistance. (D) Lollipop plot demonstrating the loop count:anchor count ratio of H3K27ac HiChIP data for each molecular subtype. (E) Volcano plot of differentially transcribed (|log2fold-change| > 1.5 & FDR <0.1) eRNAs between the Atypical and Classical subtypes after filtering transcripts for only those contained within H3K27ac HiChIP anchors (pink dots meet fold-change and FDR thresholds, purple dots meet fold-change threshold only). (F) Joined Hallmark and KEGG pathway enrichment analysis of genes linked to differentially transcribed eRNAs in (E); pathways in blue have been previously associated with BET inhibitor resistance. (G) Visualization of H3K27ac HiChIP loops at the MAP3K8 locus (left) and EGFR locus (right) in all 4 HNSCC molecular subtypes.
To assess the likely functional output of these eRNAs, we assigned each eRNA to its nearest gene and performed pathway enrichment analysis, which demonstrated an enrichment in multiple metabolic pathways, including lipid metabolism and cholesterol homeostasis, and hedgehog signaling (Figure 4B). These findings are largely in agreement with our previous enhancer-based analysis of H3K27ac-linked genes, which displayed enrichment for similar BET inhibitor resistance associated pathways (Supplementary Table S2). To extend this finding to human tumors, we leveraged data from recent publications investigating eRNA expression in TCGA tumors (Cheng et al., 2019; Chen and Liang, 2020). After assigning all the HNSCC TCGA tumor samples to their molecular subtype, we examined eRNA expression in the AT subtype compared to non-AT samples and found the AT subtype upregulated 1,867 eRNAs. After linking these eRNAs to their nearest gene, we found that, in agreement with our PRO-seq data from HNSCC cell lines, these genes were enriched for cholesterol homeostasis, hedgehog signaling, and MAPK signaling function (Figure 4C). Next, we overlapped both our HNSCC cell line PRO-seq data and the TCGA HNSCC eRNA data with the predicted enhancers for Cao et al. (Cao et al., 2017) and found that 52% (3766/7246) of the cell line PRO-seq enhancers and 67% of the TCGA eRNA enhancers (42391/63479) matched a predicted enhancer (Supplementary Figures S4). These results support the hypothesis that the AT subtype has active enhancers, as measured by eRNA expression, enriched for signaling pathways that have been demonstrated to confer resistance to BET inhibition across cancer types.
To assess the enhancer-promoter contacts in our HNSCC cell lines, we performed HiChIP for H3K27ac-marked histones to capture the E-P looping involving this enhancer mark (Supplementary Figure S5) (Mumbach et al., 2016). Consistent with the previous enhancer analyses, we discovered that the AT subtype has the highest ratio of H3K27ac-mediated loops to H3K27ac anchors across all 4 HNSCC subtypes, indicating the AT subtype may have more redundancy in its enhancer architecture than the other subtypes (Figure 4D). To assess if these loops are related to enhancer function, we overlapped our PRO-seq called eRNA enhancer regions with the H3K27ac HiChIP anchor data and performed differential expression analysis of this subset of eRNAs. We identified 48 of 57 differentially transcribed eRNAs as upregulated, and these eRNAs were associated with genes involved in MAPK signaling and lipid metabolism, such as MAP3K8, EGFR, and AGPAT4 (Figure 4E). Inspection of genes identified by this integrative analysis revealed increased contact formation between respective gene promoters and H3K27ac-marked enhancers, supporting the association of eRNA expression with active enhancers and enhancer-promoter loop formation (Figures 4F,G). Further, we compared our HiChIP data for MAP3K8, EGFR, and AGPAT4 to the Cao et al. (Cao et al., 2017) predicted enhancers, which, while demonstrating variable numbers of overlaps depending on the gene queried, maintained the enrichment of enhancer looping in the AT subtype (Supplementary Tables S5, S6). Increasing the loop call stringency by increasing the number of required paired-end tags (PETs) further exaggerated the enrichment of E-P looping in the AT subtype (Supplementary Table S6).
Overall, insights from the eRNA expression and HiChIP data support a model in which the AT subtype has more active enhancers regulating genes associated with lipid metabolism and MAPK signaling, and AT enhancers have, on average, a higher level of redundancy in their control of gene expression than non-AT enhancers by forming larger numbers of enhancer-promoter contacts.
Discussion
Here, we have demonstrated that HNSCC cell line molecular subtypes have largely similar mutational backgrounds, mutational burden, and anatomic sites of origin. In contrast, the enhancer landscapes, marked by histone H3K27 acetylation, are distinct among subtypes. In particular, we discovered the AT subtype has the highest number of enhancers and super-enhancers, as well as the most enhancer signal at common enhancer regions and a global increase in enhancer-promoter loop formation. We also demonstrate that the AT subtype is more resistant to BET inhibition and that, upon treatment with BET inhibitors, the AT subtype is able to uniquely upregulate genes associated with cell growth and BET inhibitor resistance pathways (MAPK signaling, WNT signaling, and lipid metabolism) (Rathert et al., 2015; Kurimchak et al., 2016; Iniguez et al., 2018; Cochran et al., 2019; Loganathan et al., 2019; Tonini et al., 2020; Yan et al., 2020). Further, we demonstrate a significant baseline upregulation of eRNA transcription from the enhancers of genes involved in BET inhibitor resistance pathways such as lipid metabolism and hedgehog signaling in the AT subtype (Cochran et al., 2019). Interestingly, many of these genes with increased eRNA expression in their enhancers were also found to have baseline increased enhancer-promoter looping. Together, our findings suggest that the AT subtype of HNSCC is characterized by high enhancer activity, which likely drives the expression of pathways known to confer resistance to BET inhibition.
Delineation of HNSCC into 4 subtypes was originally proposed by Walter et al. and the TCGA HNSCCC study (Walter et al., 2013; Perez Sayans et al., 2019). These two manuscripts largely focus on genomic alterations, such as copy number alterations and somatic mutations, and only one epigenetic element, in the form of DNA methylation, was assessed in the TCGA paper. As such, limited epigenomic data for HNSCC is available (Serafini et al., 2020). Despite these limitations, interest in therapies that target the epigenome continues to grow, indicating a need for more studies that focus on the epigenome of HNSCC (Alsahafi et al., 2019; Bates, 2020). The work presented here is, to our knowledge, the first to characterize the enhancer landscape of HNSCC based on the Walter/TCGA molecular subtypes. Interestingly, we identified HPV-negative samples belonging to AT subtype, which has traditionally been associated with HPV-positive or “HPV-like” samples, have increased enhancer activity compared to the non-AT subtypes. This activity is measured by increased H3K27ac peak counts, increased H3K27ac signal at common enhancer peaks, global increases in enhancer-promoter looping, and significant upregulation of eRNA expression compared to non-AT samples. This finding suggests that defining features of the AT subtype are enhancer architecture and activity - two key epigenomic aspects of HNSCC subtypes that were not previously explored. Clinical and translational significance of enhancer-based classification was shown by our recent studies in other tumor types like colorectal cancer (Orouji et al., 2021) and MPNST (Kochat et al., 2021).
Unfortunately, BET inhibitors have shown limited promise in clinical studies in solid tumors (Shorstova et al., 2021). However, in specific solid tumor contexts, such as BRD4-NUT midline carcinoma, BET inhibitors have had very encouraging results in clinical trials (Stathis et al., 2016; Piha-Paul et al., 2020; Shorstova et al., 2021). Considering these findings, identifying subsets of patients with tumor biology favorable or unfavorable to BET inhibitor response could improve the clinical utility of these compounds. Since BET inhibitors inherently rely on modulating the reader protein of H3K27ac-marked enhancers in target cells, understanding enhancer landscapes and their role in BET inhibitor response becomes an important first step in sorting patients into “favorable” or “unfavorable” groups (Stathis and Bertoni, 2018; Cheng et al., 2019). In our work, we discovered that the AT subtype is significantly more resistant to BET inhibition than other HNSCC subtypes, and that this resistance seems, at least in part, mediated by increased enhancer activity on pathways associated with lipid and cholesterol metabolism, MAPK signaling, and WNT-β-catenin signaling. Accordingly, we expect that including compounds that target these pathways in combination with BET inhibitors may sensitize otherwise resistant tumors to BET inhibition and expand the current chemotherapeutic repertoire for HNSCC treatment. As such, other enhancer/transcription blocking inhibitors, such as those against CDK9 (Zhang et al., 2018), could be tested in such enhancer-based subtypes.
We acknowledge that a limitation of our work is the focus on cell line models of HNSCC, which has certain limitations compared to studying human tumors directly. In particular, we noticed differences in BRD expression patterns between our cell line RNA-seq data and the HNSCC TCGA dataset, and the possibility of this being at least partially the result of the sample sources cannot be excluded and warrants further investigation in subsequent studies. However, given the relative sparsity of data in the HNSCC enhancer regulation space, our data can serve as a valuable resource as this field continues to grow. The work presented in this manuscript also serves as an early investigation into the enhancer regulatory landscape of HNSCC using multiple methods that can be technically challenging to perform in human tissue because of the amount of sample required and the associated difficulty of acquiring sufficient numbers of human samples. Moving forward, it will be important to perform similar studies in human tumor samples and animal models to corroborate the findings from our work in an in vivo setting.
While our work focused on HPV-negative HNSCC, our findings suggest increased enhancer activity on genes involved in lipid and cholesterol metabolism, MAPK signaling, and WNT-β-catenin signaling may serve as a general mechanism of baseline resistance to BET inhibition. Since enhancer architecture is a critical component of cell identity, it is possible that, moving forward, an assessment of a tumor’s baseline enhancer activity could serve as a potential epigenomic biomarker of response to BET inhibition and aid in tailoring treatment in a patient-specific manner (Hnisz et al., 2013; Kron et al., 2014). This could be especially useful in the case of HNSCC, where subtype-specific and tumor-specific treatments are generally lacking (Alsahafi et al., 2019).
Data availability statement
Processed data is available from the Gene Expression Omnibus under repositories GSE185531 (ChIP-seq), GSE185532 (HiChIP), GSE185533 (PRO-seq), and GSE185534 (mRNA-seq). Please note that the raw sequence data cannot be provided due to these being old cell lines without proper consents. Please contact KR for the raw sequence data. Relevant code for the manuscript can be found at: https://gitlab.com/railab/hnscc_subtypes.
Author contributions
SCC conceived, conceptualized and designed the study, planned and carried out experiments, performed data analysis, prepared figures, and wrote the manuscript. VK planned and carried out PRO-seq experiments. ZL planned and carried out ChIP-seq experiments. ATR aided in study design and analysis. MD assisted with manuscript preparation and formatting. JS, CJT, AKG, and MT provided technical help. FMJ, JW, HDS, and CRP contributed to study design and provided reagents and datasets. JNM and KR conceived, conceptualized and designed the study and wrote the manuscript. All authors edited the manuscript.
Funding
KR acknowledges support from the NIH/NCI (1RO1CA222214) and start-up from MD Anderson Cancer Center. JNM acknowledges support from the NIH/NCI (5P30CA016672-44) and a grant from the Kadoorie Charitable Foundation. CRP acknowledges support from NIH/NIDCR (5R01DE028061). This work was additionally supported by The University of Texas MD Anderson Cancer Center’s HPV-Related Cancers Moon Shot Program.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
Figures were created with Biorender.com.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.936168/full#supplementary-material
Supplementary Figure S1 | Boxplots demonstrating total number of mutations in each sample, grouped by molecular subtype (p.adj = .031 for AT vs. MS, p.adj = NS for all other comparisons).
Supplementary Figure S2 | (A) CCLE-derived JQ1 AOC values grouped by individual HNSC molecular subtype. (B) Response of HNSC cell lines to the BET inhibitor OTX015, reported as GRAOC to adjust for cell line growth rates and grouped by molecular subtype (** adj.p < 0.01). (C) Box and whisker plot of BRD4 expression in HNSCC cell lines (* = p.adj < 0.05, NS = not significant). (D) Box and whisker plot of BRD7 expression in HNSCC cell lines (* = p.adj < 0.05, ** = p.adj < 0.01, and *** p.adj < 0.005). (E) Box and whisker plots of BRD gene expression in TCGA HNSCC samples. (F) Box and whisker plot of BRD4 expression in HNSCC cell lines (* = p.adj < 0.05, NS = not significant).
Supplementary Figure S3 | GSEA analysis of (A) MDA1186 and (B) HN4 response signatures to 24 h of PLX51107 treatment at their respective GR50 values; gene sets for enrichment calculations were generated from previously published BET inhibitor response signatures.
Supplementary Figure S4 | (A) Venn diagram of overlaps between HNSCC cell line PRO-seq eRNA-based enhancers (orange) and Cao et al. (Cao et al., 2017) predicted enhancers (blue). (B) Venn diagram TCGA HNSCC eRNA-based enhancers (magenta), and Cao et al. (Cao et al., 2017) predicted enhancers (blue).
Supplementary Figure S5 | (A–D) HiC-Pro valid interaction and contact range metrics for HN4, CAL27, MDA1186, and UMSCC25 HNSCC cell lines, respectively. (E) Number of intrachromosomal paired-end tags (PETs) and distribution in distance bins per sample as reported by hichipper.
Supplementary Table S1 | List of cell lines used for H3K27ac ChIP-seq, their total H3K27ac peak counts, and their total uniquely mapped read counts.
Supplementary Table S2 | List of top 20 enriched Hallmark and KEGG pathways and their respective FDRs from genes lists generated by linking AT-unique typical enhancers to their target genes.
Supplementary Table S3 | List of top 20 enriched Hallmark and KEGG pathways and their respective FDRs from genes lists generated by linking AT-unique super enhancers (sheet 2) to their target genes.
Supplementary Table S4 | PRO-seq quality metrics for each HNSCC cell line sample as reported by the PEPPRO pipeline.
Supplementary Table S5 | Table showing relationship between Cao et al. (Cao et al., 2017) predicted enhancers and HiChIP anchors for AGPAT4, MAP3K8, and EGFR.
Supplementary Table S6 | Table showing relationship between Cao et al. (Cao et al., 2017) predicted enhancers and HiChIP H3K27ac-based loops. The “threshold” value is the number of required PET interactions used for generating Figures 4F,G (10 for AGPAT4, 5 for MAP3K8 and EGFR).
References
Agrawal, N., Frederick, M. J., Pickering, C. R., Bettegowda, C., Chang, K., Li, R. J., et al. (2011). Exome sequencing of head and neck squamous cell carcinoma reveals inactivating mutations in NOTCH1. Science 333, 1154–1157. doi:10.1126/science.1206923
Alsahafi, E., Begg, K., Amelio, I., Raulf, N., Lucarelli, P., Sauter, T., et al. (2019). Clinical update on head and neck cancer: molecular biology and ongoing challenges. Cell Death Dis. 10, 540. doi:10.1038/s41419-019-1769-9
Arnold, P. R., Wells, A. D., and Li, X. C. (2019). Diversity and emerging Roles of Enhancer RNA in Regulation of Gene Expression and Cell Fate. Front. Cell Dev. Biol. 7, 377. doi:10.3389/fcell.2019.00377
Basu, A., Bodycombe, N. E., Cheah, J. H., Price, E. V., Liu, K., Schaefer, G. I., et al. (2013). An interactive resource to identify cancer genetic and lineage dependencies targeted by small molecules. Cell 154, 1151–1161. doi:10.1016/j.cell.2013.08.003
Bates, S. E. (2020). Epigenetic Therapies for Cancer. N. Engl. J. Med. 383, 650–663. doi:10.1056/nejmra1805035
Bonner, J. A., Harari, P. M., Giralt, J., Azarnia, N., Shin, D. M., Cohen, R. B., et al. (2006). Radiotherapy plus cetuximab for squamous-cell carcinoma of the head and neck. N. Engl. J. Med. 354, 567–578. doi:10.1056/nejmoa053422
Buuren, Sv, and Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in R. J. Stat. Softw. 45, 1–67. doi:10.18637/jss.v045.i03
Cancer Genome Atlas, N. (2015). Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature 517, 576–582. doi:10.1038/nature14129
Cao, Q., Anyansi, C., Hu, X., Xu, L., Xiong, L., Tang, W., et al. (2017). Reconstruction of enhancer-target networks in 935 samples of human primary cells, tissues and cell lines. Nat. Genet. 49, 1428–1436. doi:10.1038/ng.3950
Castilho, R. M., Squarize, C. H., and Almeida, L. O. (2017). Epigenetic Modifications and head and neck Cancer: Implications for Tumor progression and Resistance to Therapy. Int. J. Mol. Sci. 18. doi:10.3390/ijms18071506
Chen, H., Li, C., Peng, X., Zhou, Z., Weinstein, J. N., Cancer Genome Atlas Research, N., et al. (2018). A pan-Cancer analysis of Enhancer Expression in nearly 9000 patient samples. Cell 173, 386–399 e312. doi:10.1016/j.cell.2018.03.027
Chen, H., and Liang, H. (2020). A high-Resolution map of human Enhancer RNA loci Characterizes super-enhancer activities in Cancer. Cancer Cell 38, 701–715 e705. doi:10.1016/j.ccell.2020.08.020
Chen, Y. W., Kao, S. Y., Wang, H. J., and Yang, M. H. (2013). Histone modification patterns correlate with patient outcome in oral squamous cell carcinoma. Cancer 119, 4259–4267. doi:10.1002/cncr.28356
Cheng, Y., He, C., Wang, M., Ma, X., Mo, F., Yang, S., et al. (2019). Targeting epigenetic regulators for cancer therapy: mechanisms and advances in clinical trials. Signal Transduct. Target Ther. 4, 62. doi:10.1038/s41392-019-0095-0
Chung, C. H., Parker, J. S., Karaca, G., Wu, J., Funkhouser, W. K., Moore, D., et al. (2004). Molecular classification of head and neck squamous cell carcinomas using patterns of gene expression. Cancer Cell 5, 489–500. doi:10.1016/s1535-6108(04)00112-6
Cochran, A. G., Conery, A. R., and Sims, R. J. (2019). Bromodomains: a new target class for drug development. Nat. Rev. Drug Discov. 18, 609–628. doi:10.1038/s41573-019-0030-7
Collisson, E. A., Bailey, P., Chang, D. K., and Biankin, A. V. (2019). Molecular subtypes of pancreatic cancer. Nat. Rev. Gastroenterol. Hepatol. 16, 207–220. doi:10.1038/s41575-019-0109-y
Creyghton, M. P., Cheng, A. W., Welstead, G. G., Kooistra, T., Carey, B. W., Steine, E. J., et al. (2010). Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl. Acad. Sci. U. S. A. 107, 21931–21936. doi:10.1073/pnas.1016071107
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). Star: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi:10.1093/bioinformatics/bts635
Dos Santos, K. C. G., Desgagne-Penix, I., and Germain, H. (2020). Custom selected reference genes outperform pre-defined reference genes in transcriptomic analysis. BMC Genomics 21, 35. doi:10.1186/s12864-019-6426-2
Eide, P. W., Bruun, J., Lothe, R. A., and Sveen, A. (2017). CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models. Sci. Rep. 7, 16618. doi:10.1038/s41598-017-16747-x
Faust, G. G., and Hall, I. M. (2014). Samblaster: Fast duplicate marking and structural variant read extraction. Bioinformatics 30, 2503–2505. doi:10.1093/bioinformatics/btu314
Forastiere, A., Koch, W., Trotti, A., and Sidransky, D. (2001). Head and neck cancer. N. Engl. J. Med. 345, 1890–1900. doi:10.1056/nejmra001375
Fragomeni, S. M., Sciallis, A., and Jeruss, J. S. (2018). Molecular subtypes and local-regional control of breast cancer. Surg. Oncol. Clin. N. Am. 27, 95–120. doi:10.1016/j.soc.2017.08.005
Gleber-Netto, F. O., Rao, X., Guo, T., Xi, Y., Gao, M., Shen, L., et al. (2019). Variations in HPV function are associated with survival in squamous cell carcinoma. JCI Insight 4. doi:10.1172/jci.insight.124762
Guinney, J., Dienstmann, R., Wang, X., de Reynies, A., Schlicker, A., Soneson, C., et al. (2015). The consensus molecular subtypes of colorectal cancer. Nat. Med. 21, 1350–1356. doi:10.1038/nm.3967
Hafner, M., Niepel, M., Chung, M., and Sorger, P. K. (2016). Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat. Methods 13, 521–527. doi:10.1038/nmeth.3853
Herz, H. M., Hu, D., and Shilatifard, A. (2014). Enhancer malfunction in cancer. Mol. Cell 53, 859–866. doi:10.1016/j.molcel.2014.02.033
Hnisz, D., Abraham, B. J., Lee, T. I., Lau, A., Saint-Andre, V., Sigova, A. A., et al. (2013). Super-enhancers in the control of cell identity and disease. Cell 155, 934–947. doi:10.1016/j.cell.2013.09.053
Hoshida, Y. (2010). Nearest template prediction: a single-sample-based flexible class prediction with confidence assessment. PLoS One 5, e15543. doi:10.1371/journal.pone.0015543
Iniguez, A. B., Alexe, G., Wang, E. J., Roti, G., Patel, S., Chen, L., et al. (2018). Resistance to epigenetic-targeted therapy engenders tumor cell vulnerabilities associated with enhancer remodeling. Cancer Cell 34, 922–938 e927. doi:10.1016/j.ccell.2018.11.005
Kalu, N. N., Mazumdar, T., Peng, S., Shen, L., Sambandam, V., Rao, X., et al. (2017). Genomic characterization of human papillomavirus-positive and -negative human squamous cell cancer cell lines. Oncotarget 8, 86369–86383. doi:10.18632/oncotarget.21174
Khan, A., and Mathelier, A. (2017). Intervene: a tool for intersection and visualization of multiple gene or genomic region sets. BMC Bioinforma. 18, 287. doi:10.1186/s12859-017-1708-7
Kochat, V., Raman, A. T., Landers, S. M., Tang, M., Schulz, J., Terranova, C., et al. (2021). Enhancer reprogramming in PRC2-deficient malignant peripheral nerve sheath tumors induces a targetable de-differentiated state. Acta Neuropathol. 142, 565–590. doi:10.1007/s00401-021-02341-z
Kron, K. J., Bailey, S. D., and Lupien, M. (2014). Enhancer alterations in cancer: a source for a cell identity crisis. Genome Med. 6, 77. doi:10.1186/s13073-014-0077-3
Kurimchak, A. M., Shelton, C., Duncan, K. E., Johnson, K. J., Brown, J., O'Brien, S., et al. (2016). Resistance to BET bromodomain inhibitors is Mediated by kinome Reprogramming in ovarian Cancer. Cell Rep. 16, 1273–1286. doi:10.1016/j.celrep.2016.06.091
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi:10.1038/nmeth.1923
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10, R25. doi:10.1186/gb-2009-10-3-r25
Lareau, C. A., and Aryee, M. J. (2018). hichipper: a preprocessing pipeline for calling DNA loops from HiChIP data. Nat. Methods 15, 155–156. doi:10.1038/nmeth.4583
Lee, T. I., and Young, R. A. (2013). Transcriptional regulation and its misregulation in disease. Cell 152, 1237–1251. doi:10.1016/j.cell.2013.02.014
Lex, A., Gehlenborg, N., Strobelt, H., Vuillemot, R., and Pfister, H. (2014). UpSet: Visualization of intersecting sets. IEEE Trans. Vis. Comput. Graph 20, 1983–1992. doi:10.1109/tvcg.2014.2346248
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). Genome project data processing S: The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi:10.1093/bioinformatics/btp352
Liao, Y., Smyth, G. K., and Shi, W. (2013). The subread aligner: Fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 41, e108. doi:10.1093/nar/gkt214
Loganathan, S. N., Tang, N., Holler, A. E., Wang, N., and Wang, J. (2019). Targeting the IGF1R/PI3K/AKT pathway sensitizes ewing sarcoma to BET bromodomain inhibitors. Mol. Cancer Ther. 18, 929–936. doi:10.1158/1535-7163.mct-18-1151
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi:10.1186/s13059-014-0550-8
Mahat, D. B., Kwak, H., Booth, G. T., Jonkers, I. H., Danko, C. G., Patel, R. K., et al. (2016). Base-pair-resolution genome-wide mapping of active RNA polymerases using precision nuclear run-on (PRO-seq). Nat. Protoc. 11, 1455–1476. doi:10.1038/nprot.2016.086
Marur, S., D'Souza, G., Westra, W. H., and Forastiere, A. A. (2010). HPV-associated head and neck cancer: a virus-related cancer epidemic. Lancet Oncol. 11, 781–789. doi:10.1016/s1470-2045(10)70017-6
Mumbach, M. R., Rubin, A. J., Flynn, R. A., Dai, C., Khavari, P. A., Greenleaf, W. J., et al. (2016). HiChIP: Efficient and sensitive analysis of protein-directed genome architecture. Nat. Methods 13, 919–922. doi:10.1038/nmeth.3999
Orouji, E., Raman, A. T., Singh, A. K., Sorokin, A., Arslan, E., Ghosh, A. K., et al. (2021). Chromatin state dynamics confers specific therapeutic strategies in enhancer subtypes of colorectal cancer. Gut 71 (5), 938–949. doi:10.1136/gutjnl-2020-322835
Perez Sayans, M., Chamorro Petronacci, C. M., Lorenzo Pouso, A. I., Padin Iruegas, E., Blanco Carrion, A., Suarez Penaranda, J. M., et al. (2019). Comprehensive Genomic review of TCGA head and neck squamous Cell Carcinomas (HNSCC). J. Clin. Med. 8. doi:10.3390/jcm8111896
Picaud, S., Leonards, K., Lambert, J. P., Dovey, O., Wells, C., Fedorov, O., et al. (2016). Promiscuous targeting of bromodomains by bromosporine identifies BET proteins as master regulators of primary transcription response in leukemia. Sci. Adv. 2, e1600760. doi:10.1126/sciadv.1600760
Pickering, C. R., Zhang, J., Yoo, S. Y., Bengtsson, L., Moorthy, S., Neskey, D. M., et al. (2013). Integrative genomic characterization of oral squamous cell carcinoma identifies frequent somatic drivers. Cancer Discov. 3, 770–781. doi:10.1158/2159-8290.cd-12-0537
Piha-Paul, S. A., Hann, C. L., French, C. A., Cousin, S., Brana, I., Cassier, P. A., et al. (2020). Phase 1 study of molibresib (GSK525762), a bromodomain and extra-terminal domain protein inhibitor, in NUT Carcinoma and other solid Tumors. JNCI Cancer Spectr. 4, pkz093. doi:10.1093/jncics/pkz093
Puissant, A., Frumm, S. M., Alexe, G., Bassil, C. F., Qi, J., Chanthery, Y. H., et al. (2013). Targeting MYCN in neuroblastoma by BET bromodomain inhibition. Cancer Discov. 3, 308–323. doi:10.1158/2159-8290.cd-12-0418
Quinlan, A. R., and Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842. doi:10.1093/bioinformatics/btq033
Rada-Iglesias, A., Bajpai, R., Swigut, T., Brugmann, S. A., Flynn, R. A., and Wysocka, J. (2011). A unique chromatin signature uncovers early developmental enhancers in humans. Nature 470, 279–283. doi:10.1038/nature09692
Ramirez, F., Ryan, D. P., Gruning, B., Bhardwaj, V., Kilpert, F., Richter, A. S., et al. (2016). deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160–W165. doi:10.1093/nar/gkw257
Rathert, P., Roth, M., Neumann, T., Muerdter, F., Roe, J. S., Muhar, M., et al. (2015). Transcriptional plasticity promotes primary and acquired resistance to BET inhibition. Nature 525, 543–547. doi:10.1038/nature14898
Rees, M. G., Seashore-Ludlow, B., Cheah, J. H., Adams, D. J., Price, E. V., Gill, S., et al. (2016). Correlating chemical sensitivity and basal gene expression reveals mechanism of action. Nat. Chem. Biol. 12, 109–116. doi:10.1038/nchembio.1986
Sartorelli, V., and Lauberth, S. M. (2020). Enhancer RNAs are an important regulatory layer of the epigenome. Nat. Struct. Mol. Biol. 27, 521–528. doi:10.1038/s41594-020-0446-0
Seashore-Ludlow, B., Rees, M. G., Cheah, J. H., Cokol, M., Price, E. V., Coletti, M. E., et al. (2015). Harnessing connectivity in a large-scale small-molecule sensitivity dataset. Cancer Discov. 5, 1210–1223. doi:10.1158/2159-8290.cd-15-0235
Serafini, M. S., Lopez-Perez, L., Fico, G., Licitra, L., De Cecco, L., and Resteghini, C. (2020). Transcriptomics and Epigenomics in head and neck cancer: available repositories and molecular signatures. Cancers Head. Neck 5, 2. doi:10.1186/s41199-020-0047-y
Servant, N., Varoquaux, N., Lajoie, B. R., Viara, E., Chen, C. J., Vert, J. P., et al. (2015). HiC-pro: an optimized and flexible pipeline for hi-C data processing. Genome Biol. 16, 259. doi:10.1186/s13059-015-0831-x
Shen, L., Shao, N., Liu, X., and Nestler, E. (2014). ngs.plot: Quick mining and visualization of next-generation sequencing data by integrating genomic databases. BMC Genomics 15, 284. doi:10.1186/1471-2164-15-284
Shi, C., Rattray, M., and Orozco, G. (2020). HiChIP-peaks: a HiChIP peak calling algorithm. Bioinformatics 36, 3625–3631. doi:10.1093/bioinformatics/btaa202
Shorstova, T., Foulkes, W. D., and Witcher, M. (2021). Achieving clinical success with BET inhibitors as anti-cancer agents. Br. J. Cancer 124, 1478–1490. doi:10.1038/s41416-021-01321-0
Smirnov, P., Safikhani, Z., El-Hachem, N., Wang, D., She, A., Olsen, C., et al. (2016). PharmacoGx: an R package for analysis of large pharmacogenomic datasets. Bioinformatics 32, 1244–1246. doi:10.1093/bioinformatics/btv723
Smith, J. P., Dutta, A. B., Sathyan, K. M., Guertin, M. J., and Sheffield, N. C. (2021). Peppro: Quality control and processing of nascent RNA profiling data. Genome Biol. 22, 155. doi:10.1186/s13059-021-02349-4
Stark, R., and Brown, G. (2020). DiffBind: Differential binding analysis of ChIP-seq peak data, 312. Bioconductor version: Release.
Stathis, A., and Bertoni, F. (2018). BET proteins as Targets for anticancer Treatment. Cancer Discov. 8, 24–36. doi:10.1158/2159-8290.cd-17-0605
Stathis, A., Zucca, E., Bekradda, M., Gomez-Roca, C., Delord, J. P., de La Motte Rouge, T., et al. (2016). Clinical Response of Carcinomas harboring the BRD4-NUT oncoprotein to the Targeted bromodomain inhibitor otx015/MK-8628. Cancer Discov. 6, 492–500. doi:10.1158/2159-8290.cd-15-1335
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550. doi:10.1073/pnas.0506580102
Terranova, C., Tang, M., Orouji, E., Maitituoheti, M., Raman, A., Amin, S., et al. (2018). An integrated platform for genome-wide mapping of chromatin states using high-throughput ChIP-sequencing in tumor tissues. J. Vis. Exp. 134 56972. doi:10.3791/56972
Tonini, C., Colardo, M., Colella, B., Di Bartolomeo, S., Berardinelli, F., Caretti, G., et al. (2020). Inhibition of bromodomain and extraterminal domain (BET) proteins by JQ1 unravels a novel epigenetic modulation to control lipid homeostasis. Int. J. Mol. Sci. 21. doi:10.3390/ijms21041297
Van Nostrand, E. L., Pratt, G. A., Shishkin, A. A., Gelboin-Burkhart, C., Fang, M. Y., Sundararaman, B., et al. (2016). Robust transcriptome-wide discovery of RNA-binding protein binding sites with enhanced CLIP (eCLIP). Nat. Methods 13, 508–514. doi:10.1038/nmeth.3810
Vasaikar, S., Huang, C., Wang, X., Petyuk, V. A., Savage, S. R., Wen, B., et al. (2019). Proteogenomic analysis of human colon cancer reveals new therapeutic opportunities. Cell 177, 1035–1049 e1019. doi:10.1016/j.cell.2019.03.030
Walter, V., Yin, X., Wilkerson, M. D., Cabanski, C. R., Zhao, N., Du, Y., et al. (2013). Molecular subtypes in head and neck cancer exhibit distinct patterns of chromosomal gain and loss of canonical cancer genes. PLoS One 8, e56823. doi:10.1371/journal.pone.0056823
Wang, J., Zhao, Y., Zhou, X., Hiebert, S. W., Liu, Q., and Shyr, Y. (2018). Nascent RNA sequencing analysis provides insights into enhancer-mediated gene regulation. BMC Genomics 19, 633. doi:10.1186/s12864-018-5016-z
Yan, G., Wang, H., Lunab, A., Bozorgui, B., Xubin Li, S M., Dereli, Z., et al. (2020). Adaptive response to BET inhibition induces therapeutic vulnerability to MCL1 inhibitors in breast cancer. biorXiv, 711895.
Yu, K., Chen, B., Aran, D., Charalel, J., Yau, C., Wolf, D. M., et al. (2019). Comprehensive transcriptomic analysis of cell lines as models of primary tumors across 22 tumor types. Nat. Commun. 10, 3574. doi:10.1038/s41467-019-11415-2
Zerbino, D. R., Johnson, N., Juettemann, T., Wilder, S. P., and Flicek, P. (2014). WiggleTools: parallel processing of large collections of genome-wide datasets for visualization and statistical analysis. Bioinformatics 30, 1008–1009. doi:10.1093/bioinformatics/btt737
Zhang, H., Pandey, S., Travers, M., Sun, H., Morton, G., Madzo, J., et al. (2018). Targeting CDK9 reactivates epigenetically silenced genes in cancer. Cell 175, 1244–1258 e1226. doi:10.1016/j.cell.2018.09.051
Zhang, Z., Lee, J. H., Ruan, H., Ye, Y., Krakowiak, J., Hu, Q., et al. (2019). Transcriptional landscape and clinical utility of enhancer RNAs for eRNA-targeted therapy in cancer. Nat. Commun. 10, 4562. doi:10.1038/s41467-019-12543-5
Zhao, M., Sano, D., Pickering, C. R., Jasser, S. A., Henderson, Y. C., Clayman, G. L., et al. (2011). Assembly and initial characterization of a panel of 85 genomically validated cell lines from diverse head and neck tumor sites. Clin. Cancer Res. 17, 7248–7264. doi:10.1158/1078-0432.ccr-11-0690
Keywords: epigenome analyses, head and neck cancer, enhancer regulation, BET inhibitors, drug resistance
Citation: Callahan SC, Kochat V, Liu Z, Raman AT, Divenko M, Schulz J, Terranova CJ, Ghosh AK, Tang M, Johnson FM, Wang J, Skinner HD, Pickering CR, Myers JN and Rai K (2022) High enhancer activity is an epigenetic feature of HPV negative atypical head and neck squamous cell carcinoma. Front. Cell Dev. Biol. 10:936168. doi: 10.3389/fcell.2022.936168
Received: 04 May 2022; Accepted: 27 June 2022;
Published: 19 July 2022.
Edited by:
Chandrima Das, Saha Institute of Nuclear Physics (SINP), IndiaReviewed by:
Borbala Mifsud, Hamad bin Khalifa University, QatarJorge Melendez-Zajgla, Instituto Nacional de Medicina Genómica (INMEGEN), Mexico
Copyright © 2022 Callahan, Kochat, Liu, Raman, Divenko, Schulz, Terranova, Ghosh, Tang, Johnson, Wang, Skinner, Pickering, Myers and Rai. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jeffrey N. Myers, am15ZXJzQG1kYW5kZXJzb24ub3Jn; Kunal Rai, a3JhaUBtZGFuZGVyc29uLm9yZw==