Molecular dissection of the oncogenic role of ETS1 in the mesenchymal subtypes of head and neck squamous cell carcinoma

Abstract

Head and Neck Squamous Cell Carcinoma (HNSCC) is a heterogeneous disease of significant mortality and with limited treatment options. Recent genomic analysis of HNSCC tumors has identified several distinct molecular classes, of which the mesenchymal subtype is associated with Epithelial to Mesenchymal Transition (EMT) and shown to correlate with poor survival and drug resistance. Here, we utilize an integrated approach to characterize the molecular function of ETS1, an oncogenic transcription factor specifically enriched in Mesenchymal tumors. To identify the global ETS1 cistrome, we have performed integrated analysis of RNA-Seq, ChIP-Seq and epigenomic datasets in SCC25, a representative ETS1high mesenchymal HNSCC cell line. Our studies implicate ETS1 as a crucial regulator of broader oncogenic processes and specifically Mesenchymal phenotypes, such as EMT and cellular invasion. We found that ETS1 preferentially binds cancer specific regulator elements, in particular Super-Enhancers of key EMT genes, highlighting its role as a master regulator. Finally, we show evidence that ETS1 plays a crucial role in regulating the TGF-β pathway in Mesenchymal cell lines and in leading-edge cells in primary HNSCC tumors that are endowed with partial-EMT features. Collectively our study highlights ETS1 as a key regulator of TGF-β associated EMT and reveals new avenues for sub-type specific therapeutic intervention.

Fig 1. ETS1 is highly expressed in HNSCC tumors and is associated with poor patient survival.

Boxplots displaying the distribution of ETS1 probe intensity, across Normal, Dysplastic and Cancerous Tissue, in the (A) Peng et al. (2011), (p = 5.16e-10) and the (B) Chen et al. (2008) microarray datasets, (p = 3.6e-6, p = 2.46e-4). (C) Bar-graph displaying the expression of ETS1 in the HNSC TCGA dataset. Transcriptomic data is from tumor samples with matched normal mucosa (n = 43) The bar-graph is displayed in ascending order (left-right) based on the expression levels of ETS1. (D) Kaplan-Meier plot of the overall survival of patients in the Chen et al. (2008) dataset (n = 97), based on the expression of ETS1 (ETS1 High, n = 49, ETS1 Low, n = 48). (E) Bar-graph displaying the distribution of HNSCC clinical stages based on the expression of ETS1 in the Chen et al. (2008) dataset, (p = 0.038, Chi-Square). X-axis denotes the number of cases.

Fig 2. ETS1 is enriched in the EMT^high Mesenchymal subtype of HNSCC across various datasets.

(A) Heatmap displaying the hierarchical clustering of the expression values of the top 50 genes specifically enriched in each subtype present in the TCGA HNSCC Dataset. (B) Bargraph highlighting enriched Biological Processes (GO Consortium) in Mesenchymal HNSCC as determined by the subtype specific DESeq2 analysis (Hypergeometric Test, FDR < 0.05). X-axis: -log(Adjusted p-value). (C) Boxplots displaying the HNSCC subtype-specific distribution of ETS1 expression and (D) EMT score in the TCGA HNSCC dataset. (E) Boxplots displaying the HNSCC subtype-specific distribution of ETS1 expression and (F) EMT score in the HNSCC PDX dataset.

Fig 3. SCC25, a representative HNSCC cell line expresses high levels of ETS1.

(A) Heatmap showing the cross-correlation value of the top 1500 most variably expressed genes in 8 HNSCC cell lines of different subtypes. The correlation matrix was reorganized via hierarchical clustering (Pearson Correlation, complete linkage). The green box highlights the two Mesenchymal classified cell lines, SCC25 and SCC4. (B) Western blot analysis of ETS1 across HNSCC cell lines. SE, short exposure; LE, long exposure; GAPDH, loading control.

Fig 4. ChIP-Seq analysis reveals direct ETS1 targets and their epigenomic features in SCC25 cells.

(A) Heatmap of the ChIP-Seq signal at the IDR-generated consensus ETS1 binding sites across three ChIP replicates. The individual ChIP-Signal matrices were subjected to hierarchical clustering and presented as a dendrogram (B) The top de novo motif, derived from the ETS1 ChIP-Seq peak-set as generated by MEME. The most similar JASPAR motif (2018 Database) as determined by TOMTOM analysis is shown above the MEME motif. (C) Distribution pattern of ETS1 binding sites across the genome. (D) Heatmap of histone modification signal density using k-means clustering on ETS1 binding sites showing three different groups. (E) Bargraphs displaying selected top enriched gene-sets (GO Biological Processes and MSigDB pathways) associated with genes annotated to 3 clusters of ETS1 ChIP-Seq peaks.

Fig 5. Super-Enhancers in SCC25 and enriched binding of ETS1.

(A) Line plot displaying the ranked H3K27Ac ChIP-Seq signal in SCC25. Representative genes marked by Super-Enhancer (SE) are shown. (B) Boxplot of RNA-Seq expression values of genes associated with typical enhancers and Super-Enhancers. (C) Bargraph of SE marked top enriched MSigDB pathway. (D) Top enriched DNA-Binding motifs of transcription factors found within the Nucleosome Free Regions (NFR) of SCC25 SEs. (E) The distribution of binding events of ETS1 with SEs displayed as a histogram. (F) Linegraph showing the differences in the average ETS1 ChIP-Seq signal measured at Nucleosome Free Regions present in either SEs or typical enhancers.

Fig 6. Loss of ETS1 affects SCC25 cell proliferation and migration.

(A) Western Blot analysis of expression of ETS1 and ETS2 in SCC25 cells either expressing either ETS1-targeting shRNAs or a non-targeting shRNA (shCON). GAPDH, loading control (B) Line plot displaying the differences in cellular proliferation between SCC25-shETS1 and SCC25-shCON cells, as determined by the MTT assay, (p = 0.0213, ANOVA, Tukey Post-Hoc Test). Wound scratch-healing assays for assessing cell migration after mitomycin treatment. (C) Percent area closure of the initial wound area of the SCC25-shETS1 and SCC25-shCON cells. (D) Representative images of the wound area after 0, 24 and 48 hours are displayed in the right (p = 0.0318, ANOVA, Tukey Post-Hoc). White hash marks denote the boundary of the wound.

Fig 7. Identification and analysis of global ETS1 target genes.

Sankey plot showing target genes of ETS1 based on the integration of RNA-Seq and ChIP-Seq data-sets. The tails of the plot display selected ETS1-dependent pathways revealed by gene-set enrichment analysis.

Fig 8. ETS1 is a key modulator of TGF-β signaling activated EMT.

(A) Overlap between selected ETS1 target genes with the TGFβ-induced EMT signature. (B) Western blot showing the effect of loss of ETS1 in SCC25 cells on TGFβ signaling and activation.

Fig 9. ETS1 is a key modulator of EMT.

(A) Heat map of direct transcriptional ETS1 target genes that overlapped with published EMT Signature (Tan et al., 2014). (B) Heatmap showing the EMT score and the relative transcript levels of EMT and MET markers within SCC25 cells expressing shETS1 compared to shCON cells. (C) Western blot showing the expression of selected markers in SCC25 cells either expressing shETS1-2 or shETS1-3 as compared to cells expressing shCON. GAPDH, loading control.