The functional GRHL3-filaggrin axis maintains a tumor differentiation potential and influences drug sensitivity.

Current therapies for treating heterogeneous cancers such as head and neck squamous cell carcinoma (HNSCC) are non-selective and are administered independent of response biomarkers. Therapy resistance subsequently emerges, resulting in increased cellular proliferation that is associated with loss of differentiation. Whether a cancer cell differentiation potential can dictate therapy responsiveness is still currently unknown. A multi-omic approach integrating whole-genome and whole-transcriptome sequencing with drug sensitivity was employed in a HNSCC mouse model, primary patients' data, and human cell lines to assess the potential of functional differentiation in predicting therapy response. Interestingly, a subset of HNSCC with effective GRHL3-dependent differentiation was the most sensitive to inhibitors of PI3K/mTOR, c-Myc, and STAT3 signaling. Furthermore, we identified the GRHL3-differentiation target gene Filaggrin (FLG) as a response biomarker and more importantly, stratified HNSCC subsets as treatment resistant based on their FLG mutational profile. The loss of FLG in sensitive HNSCC resulted in a dramatic resistance to targeted therapies while the GRHL3-FLG signature predicted a favorable patient prognosis. This study provides evidence for a functional GRHL3-FLG tumor-specific differentiation axis that regulates targeted therapy response in HNSCC and establishes a rationale for clinical investigation of differentiation-paired targeted therapy in heterogeneous cancers.

Introduction

Head and neck squamous cell carcinoma (HNSCC) is the sixth most prevalent cancer worldwide, with 600,000 new cases reported annually. HNSCC is recognized as a highly heterogeneous cancer due to patients being exposed to a variety of risk factors that include smoking tobacco, frequent alcohol use, chewing betel quid, radiation exposure, genetic predisposition, and infection with high-risk human papilloma viruses (HPVs)., Despite this heterogeneity, current conventional HNSCC treatment regimens such as radiation and chemotherapy are non-selective and are administered independent of etiology or molecular drivers. While immunotherapy has recently shown effectiveness in subgroups of HNSCC with high Programmed Death-Ligand 1 (PD-L1) expression, additional biomarkers beyond PD-L1 are needed to inform the choice of therapy. The clinical introduction of the epidermal growth factor receptor (EGFR)-inhibitor Cetuximab for treating HNSCC ushered in a new era of targeted therapy. Unfortunately, only a small subset of patients responded to Cetuximab and showed increased overall survival, while STAT3 activation played a significant role in the resistance to this therapy. Similar data showed a STAT3-driven resistance mechanism for PI3K inhibitors that is independent of the presence or absence of PIK3CA activating mutations and reinforces the urgent need for additional stratification methods for the treatment of HNSCCs.

Large-scale genomic and transcriptomic sequencing of HNSCC tumors has shown very high (∼90%) inactivating mutations in tumor suppressor genes. While TP53 is the most frequently mutated gene (TP53) in HNSCC (up to 83%), mutations of squamous differentiation genes (i.e., TP63, RIPK4, IRF6) are dominant and co-exist in the same cancer.,, These mutations are likely to drive more proliferative basal-like HNSCC phenotypes and correlate with poor patient survival., Surprisingly, the incidence of oncogene-activating mutations is low (∼20%), suggesting that dysregulation of differentiation acts as a primary driver in TP53 HNSCC. A tremendous challenge remains in translating genomic information into functional evidence for the separation of driver from passenger mutations, the identification of key impaired differentiation pathways, and the testing of biomarker-driven therapy responsiveness.

Balanced proliferation and differentiation of squamous epithelial cells maintains tissue homeostasis and prevents transformation. Activation of differentiation programs results in terminal differentiation that effectively suppresses the proliferation of squamous cells., Interestingly, the majority of the head and neck differentiation factors such as p63, RIPK4, IRF6,, and the tumor suppressor TP53 converge to regulate expression of the terminal differentiation factor Grainyhead-like 3 (GRHL3)., GRHL3 is a highly conserved epidermal-specific developmental transcription factor that functions as a major tumor suppressor in mouse and human HNSCC., Germline deletion of Grhl3 in mouse embryos results in a markedly hyperproliferative oral and skin epithelium, dysregulated epidermal differentiation, and defective wound healing and skin barrier defects, with newborn pups dying of dehydration. Interestingly, epithelial-specific conditional deletion of Grhl3 (Grhl3cKO) in mice using a keratin (K)-14 driven Cre recombinase induced spontaneous head and neck tumor development. Furthermore, ∼90% of human tumors and HNSCC cell lines show reduced levels of GRHL3 downstream of the oncogenic microRNA-21 (miR-21), resulting in hyperactivation of phosphatidylinositol 3-kinase (PI3K)/mammalian target of rapamycin (mTOR) and c-Myc signaling and conferring a poor prognosis., Since inhibitors of PI3K/mTOR and c-Myc signaling are not universally effective in TP53 HNSCC, response biomarkers are required to stratify patients likely to derive benefit and exclude those unlikely to respond. These observations prompted the analysis of a functional GRHL3-dependent differentiation pathway to identify a clinically relevant molecular vulnerability in heterogeneous HNSCC.

Results

A GRHL3-linked differentiation signature balances oral epithelial proliferation and differentiation

To explore key regulators of oral epithelial homeostasis, we performed transcriptome sequencing of normal, dysplastic, and differentiated cultured oral keratinocytes and assigned differentially expressed genes (DEGs) to four clusters (Figure 1A). Cluster 1 (C1) genes were high in differentiated cells and low in dysplastic cells, with the inverse DEG signature observed for C4. Gene Ontology (GO) functional enrichment analysis illustrated that C1 (high in differentiated cells) was enriched for keratinocyte differentiation and keratinization, while C4 (high in dysplastic cells) was enriched for G1/S transition, RNA pol II transcription, and DNA replication. These transcriptomic transitions (Figure 1B) further support suppressed differentiation and concomitant activation of proliferation. Intriguingly, within DEGs significantly enriched in the top 5 GO biological processes during differentiation (Figure 1C), GRHL3 and GRHL1 were the only transcription factors significantly upregulated in the differentiated stage (Figure 1D, right). Meanwhile, GRHL3 was significantly downregulated in dysplastic cells but no change was observed for GRHL1 expression (Figure 1D, left). Analysis of gene-gene interactions between GRHL3 and DEGs in C1 and C4 revealed a direct association with differentiation-related genes while loss of GRHL3 was linked to proliferation-related factors, respectively (Figure S1A). Knockdown of GRHL3 using short hairpin RNA (shRNA) in differentiated normal human epidermal keratinocytes (Figure S1B) and heatmaps of gene expression (Figure S1C) confirm the in-silico analysis and indicate a shift to proliferation-related DEGs (Figure S1D). These data imply that GRHL3-dependent differentiation plays an important role in the regulation of oral epithelial homeostasis and its loss promotes oral hyperproliferation and transformation.

Figure 1

A GRHL3-linked differentiation signature balances oral epithelial proliferation and differentiation

(A) Heatmap showing unsupervised hierarchical clustering of gene-expression clusters (C1 to C4) specifically expressed in dysplastic, normal, and differentiated oral epithelial cells, identified from the analysis of differentially expressed genes (DEGs) between the three stages. Top significantly enriched functional terms in each gene cluster are listed on the right. GO, Gene Ontology. (B) Principal-component analysis (PCA) plot showing gene expression transitions between dysplastic, normal, and differentiated oral epithelial cells. (C) Top 5 significantly enriched GO terms for the DEGs between normal and differentiated oral epithelial cells. (D) Volcano plots of DEGs between normal and differentiated cells (right) and between dysplastic and normal cells (left). Values with an absolute logFC > 2 and p value < 0.05 are highlighted in pink. DEGs related to epithelial differentiation are highlighted in purple and GRHL3 in red.

Loss of GRHL3 potentiates STAT3 activation in mouse and human HNSCC

The hyperactivation of STAT3 signaling is a key resistance mechanism against therapies in HNSCC. Active STAT3 is known to upregulate miR-21 in HNSCC, and miR-21 directly targets and reduces the tumor suppressor GRHL3. Expression of miR-21 was assessed in oral epithelial tissue of gp130 knockin mice, which express a mutant gp130 signal transducer resulting in constitutive Stat3 signaling. We found ∼4-fold upregulation of miR-21 in 6-month-old gp130 mice compared to wild-type (WT) control animals (Figure S2A), demonstrating a conserved correlation between STAT3 activation and elevated miR-21 in humans and mice. However, Grhl3 levels were not reduced (Figure S2B) and the gp130 mice did not develop spontaneous oral tumors. This indicates that an inverse correlation between Grhl3 levels and STAT3 activation is HNSCC-specific. We then examined normal whole-tongue tissue and spontaneous tumors from Grhl3-deficient mice using whole transcriptome sequencing. Interestingly, gene set enrichment analysis (GSEA) identified a significant enrichment for genes involved in the inflammatory response (Figure 2A) and cytokines and chemokines associated with STAT3 signaling in Grhl3cKO tumors (Figure 2B). Analysis by immunohistochemistry (IHC) revealed strong nuclear staining for p-Stat3 in Grhl3cKO tumors, whereas the normal adjacent cells were negative (Figure 2C). Western blot (WB) analysis demonstrated upregulation of active Stat3 (pY705) in tumor samples (Figure 2D) along with increased mRNA expression of the STAT3 target genes, Cish, Icam1, and Irf1 in the Grhl3cKO tumors (Figure 2E). The correlation between miR-21/GRHL3 expression and STAT3 pathway activation was further assessed in human samples using The Cancer Genome Atlas (TCGA) Research Network-HNSCC dataset (n = 279). Compared to normal adjacent oral tissues, the matched primary tumor samples displayed a significant ∼2-fold upregulation of miR-21 and ∼4-fold reduction in GRHL3 gene expression (Figure S2C). Pearson analysis revealed a significant (R = –0.26, p < 0.0001) negative correlation in human HNSCC (Figure S2D). Furthermore, a similar enrichment for genes involved in the inflammatory response and STAT3 signaling was identified in patients with high miR-21 and low GRHL3 expression (Figure 2F) along with significant overexpression of the STAT3 target genes CISH, ICAM1, and IRF1 compared to patients with low miR-21 and high GRHL3 expression (Figure 2G). Expression of these genes was not significantly altered between the group of patients with low miR-21 and low GRHL3 and those with high miR-21 and high GRHL3 (Figure S2E). These data indicate that loss of GRHL3 downstream of miR-21 correlates with the activation of STAT3 signaling in mouse and human HNSCCs.

Figure 2

Loss of GRHL3 potentiates STAT3 activation in mouse and human HNSCC

(A) A significant enrichment for inflammatory genes is observed in Grhl3cKO tumors compared to normal adjacent tissues. The significance was empirically determined by 1,000 gene-set permutations (normal adjacent tissues, red; tumors, blue). (B) Expression of cytokines/chemokines is shown as a heatmap. (C) p-Stat3 IHC demonstrates hyperactive Stat3 in tumors but not in adjacent tissues from Grhl3cKO mice. Scale bars correspond to 50 μm. (D) WB for p-Stat3 on spontaneous tongue tumor tissues compared to adjacent control tissues from 6-month-old Grhl3cKO mice. (E) Upregulation of the STAT3 target genes Cish, Icam1, and Irf1 was quantitated by qPCR for tumors from Grhl3cKO mice normalized to Gapdh and relative to adjacent (normal) tissue. The results are presented as mean ± SEM. The statistical values were considered significant at p < 0.05. ∗p ≤ 0.05, and ∗∗p ≤ 0.01. (F) GSEA of the TCGA-HNSCC RNA-seq data. Patients were divided into two groups: high miR-21, low GRHL3 (n = 79, in blue), and low miR-21, high GRHL3 (n = 62, in red). The significance was empirically determined by 1,000 gene-set permutations. (G) Analyses of STAT3 target gene expression in the same patient groups. The comparison between the two groups was performed using an unpaired Student’s t test (∗p < 0.05; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001). Data are shown as mean ± SEM. Patients HNSCC with high miR-21, low GRHL3 expression show STAT3 pathway activation compared to patients with low miR-21, high GRHL3 expression. (H and I) WB analyses of GRHL3 expression and STAT3 activation in cell lines treated for 2 h with (H) 100 ng/mL IL-6 or (I) 100 ng/mL EGF. While STAT3 activity is undetectable in SCC25 (boxed), these cells do not respond to IL-6 or EGF treatment. (J and K) OKF-6 cells with shRNA-mediated GRHL3 knockdown respond to (J) IL-6 stimulation, or (K) EGF stimulation, and show a significant increase in STAT3 activation compared to cells transduced with a scrambled control. β-actin was used as a loading control.

To evaluate the effects of GRHL3 loss on STAT3 activation, we characterized the molecular profiles of human HNSCC cell lines cultured under normal conditions. Expression of GRHL3 was measured by qPCR for a series of cell lines (SCC4, CAL27, SCC25, SCC9, A253, and FaDu). The gene expression data mirrors the degree of heterogeneity observed in primary human HNSCC with the highest GRHL3 expression identified for SCC25 (Figure S2F, top panel). We next examined the commonly activated signaling pathways by WB (Figure S2F, lower panel). Consistent with previous findings, PI3K/mTOR signaling was hyperactive (phospho-S6 positive) and c-Myc was overexpressed in all HNSCC lines. Furthermore, SCC25 showed increased expression and activation of the translation initiation factor eIF4E, a downstream effector of PI3K/mTOR and c-Myc signaling. ERK and STAT3 were also examined under basal culture conditions. Compared to normal OKF6 cells, p-ERK was consistently low in all SCC lines examined. While p-STAT3 levels were variable, they were negative in SCC25 and SCC9 (Figure S2F) and these lines showed differential expression of pro-inflammatory genes related to STAT3 signaling (Figure S2G). We next stimulated the HNSCC lines with interleukin-6 (IL-6) and EGF to activate the STAT3 pathway, as this can induce adaptive resistance to targeted therapy. Interestingly, SCC25 was the only line non-responsive to 100 ng/mL of IL-6 (Figure 2H) or EGF (Figure 2I) (i.e., complete lack of p-STAT3), suggesting a potential inhibitory mechanism unique to SCC25, which has high GRHL3 expression and lacks STAT3 activation. To validate this, OKF-6 cells transduced with either a scrambled control or GRHL3 shRNA and treated with IL-6 (Figure 2J) or EGF (Figure 2K) showed a significant increase in p-STAT3 in shGRHL3 compared to control cells. Furthermore, we show that the SCC1 cancer cell line has comparable GRHL3 expression to the normal oral cell line OKF-6 (Figure S3A). Interestingly, shRNA-mediated knockdown of GRHL3 in SCC1 (Figure S3B) induced significant upregulation of p-STAT3 that was further increased following treatment with IL-6 but not EGF (Figures S3C and S3D). Collectively, our data support a role for GRHL3 in the negative regulation of IL-6/STAT3 signaling and indicate that loss of GRHL3 promotes hyperactive STAT3 in both normal and HNSCC cell lines.

Mutational analysis and drug sensitivity of HNSCC cell lines

We investigated the role of GRHL3 in acquisition of drug sensitivity by integrating whole-genome mutational analysis, whole-transcriptome and drug sensitivity data for HNSCC patient samples, and cell lines retrieved from the TCGA Research Network, Cancer Cell Line Encyclopedia (CCLE; Broad Institute), COSMIC, and Genomics of Drug Sensitivity in Cancer (GDSC) databases. The inhibitors of PI3K/AKT/mTOR and c-Myc downstream signaling were prioritized because these pathways are hyperactive in all the cell lines (Figure S2F). While a generally broad range of the half maximum inhibitory concentration (IC50) values was observed (Figure 3A), SCC25 showed the highest sensitivity to drugs in all of the classes. Data retrieved from the CCLE database showed that SCC25 harbors the lowest number of somatic mutations, while A253 and FaDu have the highest burden (Figure 3B). Interestingly, SCC9 is the only cell line with WT TP53 (Figure 3B, lower panel; Table S1) demonstrating sensitivity to most inhibitors while the other resistant TP53 mutant lines mirror patients with TP53 mutant variants who have poor treatment responses and reduced overall survival.

Figure 3

Mutational analysis and drug sensitivity of HNSCC cell lines

(A) Heatmap of drug sensitivity and resistance to inhibitors of STAT3, PI3K/AKT/mTOR, c-Myc, and their downstream kinases in the HNSCC cell lines. The drug log-IC50 was adapted from the Genomics of Drug Sensitivity in Cancer (GDSC) database. (B) The top panel shows the type and rate of mutations from the Cancer Cell Line Encyclopedia (CCLE) Cell Line mutation and COSMIC databases. The bottom panel depicts the most common somatic mutations in TCGA-HNSCC patients that were also present in the HNSCC cell lines. (C) qPCR analysis of the highly mutated genes (TP53, CDKN2A, and LRP1B) in the HNSCC cell lines treated with the selected inhibitors for 24 h. mRNA expression levels were normalized to untreated cells. The results are presented as mean ± SEM. The statistical values were considered significant at p < 0.05. ∗p ≤ 0.05, ∗∗p ≤ 0.01, ∗∗∗p ≤ 0.001, and ∗∗∗∗p ≤ 0.0001.

To uncover predictive therapeutic vulnerabilities, we assessed 13 genes (ASTN1, BRD1, CARD10, CEP78, EPHB3, FANCM, LRRK2, MUC20, RYR3, SAP130, SRRM2, TMX4, and TTN) showing common hotspot mutations in the HNSCC lines (Figure S4A) for their basal expression. Gene-expression profiling did not show a pattern that correlates with their mutational status in the mutant cell lines (Figure S4B). Furthermore, treatment with the small molecule inhibitors did not consistently induce (or repress) their expression in treated cells (Figure S5).

Interestingly, within the TP53 mutant HNSCC cell lines, SCC25 was the most sensitive to all tested inhibitors (Figure 3A). This drug-sensitive cell line is the only one with a WT filaggrin (FLG) gene (Table S1), and FLG is within the top commonly mutated genes (∼16%) in TCGA Research Network-HNSCC patient samples (Figure 3B). To assess whether FLG is associated with SCC25’s sensitivity to small molecule inhibitors, we quantitated FLG expression after treatment with drugs targeting the PI3K/AKT/mTOR (Dactolisib and SEL201), c-Myc (JQ1 and CX-5461), and STAT3 (Niclosamide) signaling pathways. Strikingly, FLG expression was robustly induced in response to these inhibitors only in SCC25 compared to the other FLG mutant HNSCC cells (Figure 3C), and more importantly, FLG levels correlated with the drug sensitivity of SCC25 to each of the PI3K/AKT/mTOR and c-Myc inhibitors (Figure 3A). These findings suggest that WT FLG is a predictive biomarker of responsiveness to PI3K/AKT/mTOR and c-Myc inhibitors.

A functional GRHL3-FLG differentiation pathway promotes sensitivity to targeted therapy

FLG is an important terminal differentiation gene in squamous tissues and its deregulation results in defective epithelial barrier formation and atopic dermatitis. Under normal (non-stimulated) cell culture conditions, FLG expression was the highest in SCC25 (Figure S6A). To assess whether FLG-dependent differentiation plays a role in HNSCC, we induced cellular differentiation by treating the cell lines with 2 mM CaCl2 (24 h). Both GRHL3 and FLG were significantly upregulated following Ca2+-induced differentiation in SCC25 only, reflecting a functional differentiation response in these cells (Figures 4A and 4B). This was further confirmed by growing SCC25 to >100% confluency (14 days), which resulted in differentiation-dependent GRHL3-FLG induction (Figure S6B). However, the expression of other terminal differentiation markers, such as GRHL1 and TGM1, did not increase in response to calcium treatment, indicating that functional differentiation in SCC25 might be entirely restricted to (or only dependent on) GRHL3 and FLG expression (Figure S6C).

Figure 4

A functional GRHL3-FLG pathway promotes sensitivity to targeted therapy

(A and B) qPCR expression of (A) GRHL3, and (B) FLG in HNSCC cell lines treated with 2 mM CaCl2 (24 h) and normalized to GAPDH, relative to control untreated cells show significant induction only in SCC25. (C) Illustration of the CRISPR-dCas9 GRHL3-activation SAM system with the potential GRHL3 binding site in FLG promoter. (D) Significant induction of GRHL3 and FLG mRNA, but not of TGM1 (not significant, n.s.) in 293T cells transduced with the CRISPR-dCas9 GRHL3-activation SAM system. mRNA expression levels were normalized to GAPDH relative to control empty sgRNA vector. (E) FLG expression following Dox (1 μM, 72 h) treatment of SCC25 cells transduced with an inducible lentivirus CRISPR-Cas9 system to mediate FLG knockdown. Two single guide RNAs (sgRNAs) were employed to target the FLG gene at different sites. Compared to control (Dox negative) cells, downregulation of FLG (~50%) was achieved with both sgFLG.1 and sgFLG.2 following Dox induction. (F) SCC25 cells grown for 96 h in the presence of Dox (1 μM) did not show any proliferative difference to untreated cells. (G) Downregulation of FLG renders the treatment-sensitive SCC25 resistant to all small molecule inhibitors. (H) WB analyses of p-STAT3 in SCC25 treated with the inhibitors at IC50 in the absence of FLG for 3 days. p-STAT3 was increased in cells treated with c-Myc and PI3K/mTOR inhibitors and absent in those treated with Niclosamide. β-actin was used as the loading control. The results are presented as mean ± SEM. The statistical values were considered significant at p < 0.05. ∗p ≤ 0.05, ∗∗p ≤ 0.01, ∗∗∗p ≤ 0.001, and ∗∗∗∗p ≤ 0.0001.

We next sought to identify whether GRHL3 transcriptionally activates the FLG gene. The GRHL3 DNA consensus-binding site (5′-AACCGGTT-3′) has been conserved across 700 million years of evolution. Using the Jaspar database, we identified a potential GRHL3-binding site in the FLG promoter (5′-TACAGGTT-3′) located 173 bp upstream of the TSS, which is conserved in human and mouse (Figure 4C). To test whether GRHL3 activates FLG in vitro, we induced GRHL3 expression from its endogenous promoter using a two-component CRISPR-activation SAM system. A nuclease dead Cas9 (dCas9) fused to a VP64 transcriptional activator domain promoted GRHL3 overexpression in human HEK293T cells, leading to specific induction of its downstream FLG target, but not of its TGM1 target (Figure 4D).

Direct binding of GRHL3 to the FLG promoter was confirmed in vitro by electrophoretic mobility shift assay (EMSA; Figure S6D). Furthermore, RNA sequencing (RNA-seq) analysis of Grhl3-deficient mouse tongues showed reduced FLG levels (∼4 fold) only in tumors when compared to adjacent tissues and WT mice (Figure S6E). These data indicate that FLG is a critical downstream target of GRHL3 in HNSCC differentiation.

In parallel, we measured the enzymatic activity of the aldehyde dehydrogenase (ALDH) stem cell marker in the HNSCC cells using the AldeRed detection assay (Figures S7A and S7B). SCC25 and FaDu, which have the highest ALDH activity among the HNSCC lines, were grown on ultra-low attachment plates for 7 days. Compared to FaDu, SCC25 formed relatively small spheroids (Figure S7C), suggesting that functional differentiation may limit the 3D-growth of SCC25. We then examined whether loss of functional GRHL3-FLG-dependent differentiation affects SCC25 drug sensitivity. Dox-induced CRISPR-Cas9 with 2 separate single guide RNAs against FLG (sgFLG), confirmed by qPCR (Figure 4E), did not affect the proliferation of SCC25 (Figure 4F). Interestingly, the loss of FLG in SCC25 resulted in a dramatic increase in resistance to the PI3K/AKT/mTOR and c-Myc inhibitors compared to sensitive control cells (Figure 4G). WB analysis showed increased p-STAT3 in sgFLG cells treated with these inhibitors (Figure 4H). We next validated these findings in SCC1 cells that are FLG WT (CCLE and COSMIC databases). This cell line was able to differentiate in response to calcium treatment and showed significant upregulation of both GRHL3 and FLG (Figure S3E). When FLG was knocked down in SCC1 (Figure S3F), similar to SCC25, there was no effect on cell proliferation (Figure S3G), but importantly, a significant reduction in the sensitivity to PI3K/AKT/mTOR and c-Myc inhibitors was observed (Figure S3H). Moreover, increased resistance to these inhibitors correlated with increased p-STAT3 (Figure S3I) compared to control untreated cells. These results support the data in SCC25 and indicate that FLG is a biomarker for treatment response and its loss induces resistance to targeted therapy that correlates with STAT3 activation.

A GRHL3-FLG signature predicts favorable prognosis for HNSCC patients

Expression levels of GRHL3 and FLG were significantly higher in both moderate and well differentiated HNSCC compared to poorly differentiated patient HNSCC (Affymetrix array; Figure 5A), who have worse survival outcomes., Interestingly, expression levels of the terminal differentiation markers GRHL1 and TGM1 were not consistently correlated with moderate/well differentiated HNSCC compared to the poorly differentiated ones (Figure S6F), suggesting that increased expression of GRHL3 and FLG is favorably associated with tumor differentiation. To determine whether the functional GRHL3-FLG axis is clinically relevant, RNA-seq expression data of human HNSCC was extracted from the TCGA Research Network-HNSCC dataset and revealed that while a significant positive correlation exists for GRHL3 and FLG, their expression is negatively correlated to CISH, ICAM1, and IRF1 (Figure 5B). Supervised clustering showed that clusters with high or low GRHL3-FLG levels are inversely correlated with STAT3 target gene expression (CISH, ICAM1, and IRF1) in patient HNSCC (Figure 5C). Furthermore, the survival analyses of patients (n = 271) with primary HNSCC displayed no significant difference between GRHL3 (1/3 of patients with the lowest GRHL3 expression) and GRHL3 (compared to adjacent tissue) or between FLG and FLG; but importantly, a survival advantage of the GRHL3-FLG cohort was evident in comparison to the GRHL3-FLG patients (Figure 5D).

Figure 5

A GRHL3-FLG signature predicts favorable prognosis for HNSCC patients

(A) The expression of GRHL3 and FLG is significantly higher in moderate (n = 38) and well-differentiated (n = 19) compared to poorly differentiated (n = 21) HNSCC samples (Affymetrix HG133 arrays) and normalized to GAPDH. The results are presented as mean ± SEM. The statistical values were considered significant at p < 0.05. The results are presented as mean ± SEM. The statistical values were considered significant at p < 0.05. ∗p ≤ 0.05, and ∗∗p ≤ 0.01. (B)Significant correlation between GRHL3 and FLG and STAT3-target (CISH, ICAM1, and IRF1) gene expression from the TCGA-HNSCC patient database (n = 279; Pearson’s correlation). (C) Supervised clustering of high GRHL3 and FLG levels shows a clear subtype (boxed) with an inverse correlation to STAT3-target gene (CISH, ICAM1, and IRF1) expression. Rows are centered and unit variance scaling is applied to rows. (D) Kaplan-Meier survival analyses of HNSCC patients stratified based on their GRHL3 expression level and FLG mutation status. The difference in survival rates between patients with GRHL3 (n = 24) and GRHL3 (n = 247) or between patients with FLG (n = 227) and FLG (n = 44) was n.s. The survival of the cohort with GRHL3 and FLG (n = 15) was significantly superior to that of patients with GRHL3 and FLG (n = 35; p = 0.04771, log-rank test).

Collectively, these data establish a functional GRHL3-FLG differentiation axis as a favorable prognostic predictor and a response biomarker in HNSCC and indicate that loss of GRHL3-FLG promotes an acquired STAT3-dependent resistance to targeted therapy.

Discussion

Many studies have focused on genomic biomarker-driven targeted therapies, which have not proven totally effective for predicting responsiveness. Key examples include EGFR-T790M mutations and resistance to EGFR inhibitors in EGFR mutant lung cancer, ESR1 mutations in estrogen-receptor-positive breast cancer treated with endocrine therapy and reversions of pathogenic mutations in BRCA1- and BRCA2-deficient cancers treated with poly-ADP ribose polymerase (PARP) inhibitors. Differentiation-paired targeted therapy has not been investigated, and whether it promotes therapy response remains to date largely unknown. Mutations in TP53 synergizes with the loss of differentiation pathways to facilitate tumor progression. In agreement with this, upregulation of TP53 was shown to induce differentiation and suppress cutaneous SCC in mice, and inhibitors of PI3K/mTOR signaling provided anti-tumor activity through the induction of WT TP53 in human xenograft and murine models of HNSCC. This is recapitulated in the TP53 SCC9 cells that showed high sensitivity to multiple inhibitors compared to TP53 cells (harboring FLG; Figure 3A). Differentiation mechanisms are therefore promising response biomarkers for the stratification of therapies in both loss-of-function (LOF) and gain-of-function (GOF) TP53 HNSCC.

Substantial efforts have been made to stratify HNSCC patients into subtypes based on their HPV status, EGFR expression and activation, and CD8+ lymphocyte tumor infiltration., However, the survival rate of these patients is still significantly low. While HPV-positive patients generally show better responses to therapies and favorable prognosis, an inflammatory response may contribute to the resistance of HPV-negative HNSCC to therapies. These HNSCC are commonly induced by heavy smoking and alcohol consumption and have been associated with chronic inflammation, cytokine-activated receptors, and STAT3 pro-inflammatory signaling. Recent studies have shown that cytokines are also secreted by epithelial cells., Importantly, SCC25 cells are non-responsive to IL-6 (or EGF) and display minimal STAT3 phosphorylation. We propose that the developmental factor GRHL3, which is highly expressed in SCC25 cells, directly regulates STAT3-related phosphatase gene transcription. A GRHL3 consensus binding motif was identified in the type-1 protein phosphatase (PP1) regulatory subunit (PPP1R3F) genes’ proximal first intron, and this is conserved across the placental mammals. We showed direct binding of GRHL3 to the PPP1R3F putative regulatory region by a chromatin immunoprecipitation (ChIP) assay (Figures S8A and S8B). Moreover, reduced GRHL3 expression, both in vivo and in vitro, correlated with reduced PPP1R3F mRNA expression (Figures S8C and S8D). Therefore, a GRHL3-PPP1R3F-STAT3 regulatory mechanism may operate (Figure S9) since PP1 inhibition induces phosphorylation of STAT3.

Our multi-omic approach integrating whole-genome, whole-transcriptome, and drug sensitivity data strongly implicates the downregulation of GRHL3 and FLG with STAT3 pathway activation. While STAT3 is active in most HNSCC samples, STAT3 GOF mutations have not been observed nor have activating mutations in upstream regulators., Conversely, FLG is mutated in ∼16% of HNSCC. Our findings demonstrate that loss of GRHL3-FLG-mediated differentiation in HNSCC potentiates a pre-existing STAT3 pro-inflammatory pathway, resulting in resistance to small molecule inhibitors. This is shown in cells initially sensitive to targeted therapies (SCC25) that became resistant along with an increase in p-STAT3 (Figures 4G and 4H). On the other hand, expression of FLG was significantly induced in sensitive cells treated with the inhibitors. Collectively, these data indicate that the terminal differentiation GRHL3-FLG axis predicts response to PI3K/AKT/mTOR and c-Myc inhibition, and its induction by targeted therapy dramatically reduces cell growth. Since Cetuximab-resistant SCC25 cells also express low GRHL3-FLG, this mechanism of action may be valid for the sensitization of HNSCC cells to a broader range of inhibitors.

Recently, immunotherapy has resulted in improved HNSCC patient survival by exploiting inhibitory checkpoint pathways that suppress anti-tumor T cell responses. However, only a minority of patients derive benefit for non-selective immunotherapies. Because Stat3 signaling was shown to regulate PD-1/PD-L1 expression and the anti-tumor immune response, future investigations will determine whether combining immunotherapy with differentiation-paired targeted therapy may provide greater therapeutic benefit.

In summary, this study identifies a functional GRHL3-FLG differentiation axis in a subset of HNSCC. GRHL3-FLG upregulation in response to targeted therapy efficiently reverses cancer cell proliferation, while disruption of this axis confers resistance to PI3K/AKT/mTOR and c-Myc inhibitors (Figure S9). These findings have exciting therapeutic implications with the potential to evaluate this axis as a predictive prognostic and response biomarker in the clinic for the therapeutic stratification of heterogeneous HNSCC, notably those with TP53 mutations for which there are currently modest clinical response rates to most available treatment options.

Materials and methods

RNA-seq analysis

3′ RNA-seq was conducted by the Molecular Genomics Core (PMCC, Melbourne). The total RNA quantity was measured using Qubit RNA HS (Thermo Fisher Scientific). 500 ng total RNA was used for library preparation according to standard protocols (QuantSeq 3′ mRNA-seq FWD, Lexogen). Indexed libraries were pooled and sequenced on a NextSeq500 (Illumina). 5–15 million single-end 75 bp reads were generated per sample. Sequenced reads were trimmed and aligned to hg38 genome via Cutadapt and HISAT2. Gene counts were obtained from featureCounts. Differential expression was performed using Limma. Genes were considered as DEGs if the absolute log-fold change was >2 and p value <0.05. All analysis packages were operated within the Galaxy suite environment. The functional enrichment analysis of annotated terms from GO was performed with the online tool DAVID using the human genome as the background. The gene-gene interaction analysis was completed using the GeneMania app in Cytoscape 3.

Drug-sensitivity assays

HNSCC cells grown at log phase were harvested, counted, and 8,000 cells seeded in 96-well plates. Cells were allowed to attach and grow overnight and were then treated for 3 days with drugs at IC50 values. To calculate the IC50 of each drug, we calculated the percentage of cell growth normalized to a vehicle control using the following equation: 100 × (T–T0) ⁄ (C–T0), where C is the optical density (OD) of the control. For SCC25, IC50 values were as follow for Niclosamide, 0.16 μM; JQ1, 0.8 μM; Dactolisib, 0.1 μM; SEL201, 4 μM; and CX-5461, 1.2 μM. Cell viability was measured using Alarma Blue according to the manufacturer’s instructions. To calculate the cell-killing rate, we normalized absorbance readings to controls (without cells or without treatment).

For gene-expression analysis, HNSCC cells were seeded in 12-well plates at a density of 100,000 cells/well and treated with drugs at IC25 or IC50 for 24 h followed by RNA extraction. For western blotting, SCC25 with/without reduced FLG were seeded in 6-well plates at a density of 200,000 cells/well and treated with drugs at IC50 for 72 h followed by protein extraction.

Drug sensitivity and proliferation of the HNSCC cells were measured using the Sulforhodamine B (SRB) assay. Cells were seeded into 96-well plates (3,000–8,000 cells/well) and 24 h later serial dilutions of the drugs dissolved in DMSO were added independently. One plate per cell line was used to estimate the OD at the start (T0) in the absence of any treatment (DMSO vehicle only). After 3 days, absorbance readings of treated plates were measured. The plates were fixed and then stained with SRB, and the OD values of cells exposed to each drug concentration were measured at a wavelength of 564 nm (T).

Animal studies

Grhl3cKO and gp130(Y757F) mice were maintained on a C57BL/6J background. All animal studies followed the National Health and Medical Research Council (NHMRC) Australian Code of Practice for the Care and Use of Animals for Scientific Purposes and were approved by the Animal Experimentation Ethics Committee at the PMCC (E587). Mice were euthanized by CO2 asphyxiation, the tongues dissected, and tumors and adjacent tissues removed for pathological and gene-expression analysis. gp130 knockin mice were a kind gift from Dr. Michael Buchert (Olivia Newton-John Cancer Research Institute, Heidelberg, VIC, Australia).

Bioinformatic datasets and analyses

Publicly available datasets were retrieved from TCGA Research Network using the TCGA-biolinks R package. Gene expression (RNA-seq and miR-seq), mutation status, and patient survival information for 279 HNSCC samples were analyzed. Gene mutation profiling of the cell lines were retrieved from the CCLE (release date 02-Jan-2019). Publicly available HNSCC microarray data from ArrayExpress: E-TABM-302 based on Affymetrix HG133 arrays were retrieved from ArrayExpress. Publicly available RNA-seq data (GEO: GSE37049) of normal human keratinocytes transduced with shRNA GRHL3 were retrieved from Gene Expression Omnibus.

Statistical analyses

Statistical analysis was performed using Prism 8 (GraphPad). Statistical significance was assessed using the unpaired Student’s t test or ANOVA and multiple comparison test or log-rank (Mantel-Cox) test for survival analysis. The results are presented as mean ± SEM. The statistical values were considered significant at p < 0.05. ∗p ≤ 0.05, ∗∗p ≤ 0.01, ∗∗∗p ≤ 0.001, and ∗∗∗∗p ≤ 0.0001.