Transcriptome analysis of EpEX- and EGF-induced EGFR activation in HNSCC cell lines

Kyse30 squamous cell esophageal carcinoma cells and FaDu squamous cell pharyngeal carcinoma were chosen to delineate effects of EpEX and EGF on EGFR signaling based on their responsiveness to EGFR-mediated EMT. Early effects of treatment have been observed at the level of ERK1/2 phosphorylation up to six hours, whereas morphologic and functional changes occurred after 48-72 h [21]. Bulk 3´RNA-seq transcriptome analysis of 6 and 72 h treatments established an EGFR-mediated EMT gene signature, which formed the basis of a signature transfer to clinical cohorts to assess prognostic, therapeutic, and predictive markers (Fig. 1A).

Fig. 1

RNA-seq analysis of EGFR-mediated EMT. A Workflow: Kyse30 and FaDu cells were treated with EGF and EpEX inducing proliferation or EMT, respectively. Early transcriptional changes were assessed after 6 h and transcriptional differences associated with cellular effects after 72 h. EGFR-mediated EMT-associated genes serve to define prognostic gene signatures in clinical cohorts. Potential therapeutic and predictive markers in the EGFR-mediated EMT gene signature are explored via functional characterization and bioinformatic approaches. B Copy number variation and gene expression values of cell lines of the upper aerodigestive tract (esophageal carcinoma and HNSCC) were extracted from CCLE and are depicted as dot plot. C Kyse30 and FaDu cells were stained with EGFR-specific antibodies in combination with FITC-labeled secondary antibody. Shown are representative histograms (left panels) and mean expression values with SD of n = 5 independent experiments; * p-value < 0.05 (t-test). D-E Representative cytological pictures of Kyse30 and FaDu cells treated with the indicated components at 6 h and 72 h are shown (n = 4 independent experiments). Ctrl.: control treatment under serum-free conditions; Fc: recombinant immunoglobulin Fc region; EGFlow: 1.8 nM EGF; EGFhigh: 9 nM EGF; EGFhigh + EpEX: 9 nM EGF in combination with 50 nM EpEX-Fc. Scalebars represent 100 µm. F Principal component analysis of 3´-RNA-seq samples in Kyse30 and FaDu cells with the indicated treatments are shown (n = 4)

FaDu cells have neither amplifications nor mutations in EGFR [26], while Kyse30 show amplification of EGFR without mutations [27]. Copy number variations of Kyse30 (CNV + 1) and FaDu cells (CNV + 0) were assessed in datasets of the cancer cell line encyclopedia (CCLE) (Supplementary Fig. 1). CNV correlated with increased EGFR gene expression in Kyse30 compared to FaDu (Fig. 1B) and a 3.2-fold higher expression of EGFR in Kyse30 cells compared to FaDu cells was observed (Fig. 1C).

Serum-starved Kyse30 and FaDu cells were treated with EpEX-Fc (50 nM), EGFlow (1.8 nM), EGFhigh (9 nM), or a combination of EpEX-Fc and EGFhigh. As controls, cells remained either untreated under serum-free conditions (controls for EGF treatment) or were treated with recombinant immunoglobulin Fc-region that served as control for EpEX-Fc. Fc allowed to express EpEX-Fc primarily in a dimeric form, which represents its natural state [21, 28]. Treatment with EGFhigh resulted in the induction of a mesenchymal morphology after 72 h and co-treatment with EpEX inhibited EGFhigh-induced EMT (Fig. 1D-E). Principal component analysis showed distinct clustering of 0 and 6 h control-treated cells (0–6 h Ctrls.), 6 h samples independently of the treatment (6 h Treatment), 72 h treatments (72 h non-EMT), except EGFhigh-treated cells, which clustered separately (72 h EMT) (Fig. 1F).

EGFhigh induces sustained transcriptional activation

Gene set enrichment analysis (GSEA) was conducted for each treatment condition compared to controls, i.e. serum withdrawal in absence or presence of Fc. Significantly enriched gene ontology (GO)-terms were conserved for all four treatment modalities at 6 h and they primarily addressed ribosome biogenesis in Kyse30 cells and cell migration and motility in FaDu cells (Supplementary Fig. 2A). At 72 h, significantly enriched GO-terms were only identified for EGFhigh treatment in Kyse30 cells, and for EGFhigh and EGFlow treatments in FaDu cells. Significantly suppressed GO-terms in Kyse30 cells were related to DNA replication and nuclear division, whereas activated terms were related to keratinization, vesicle-mediated transport and cytokine signaling. Activated GO-terms in FaDu cells covered cell adhesion, cytoskeleton organization, and wound healing, while cell cycle-associated GO-terms were suppressed (Supplementary Fig. 2A). A similar activation of ribosome biogenesis at 6 h and a suppression at 72 h was observed in Kyoto encyclopedia of genes and genomes (KEGG) terms for Kyse30 cells. KEGG term “DNA replication” was suppressed in FaDu cells at 72 h, while KEGG terms significantly regulated at 6 h were heterogeneous and primarily associated with cytokine signaling (Supplementary Fig. 2B).

Molecular signatures database (MSigDB) hallmark gene sets of biological processes or states from the BROAD institute were further chosen for GSEA. MSigDB hallmarks are composed of genes identified in well-characterized, ground-truth-based cellular and animal systems for 50 major biological processes in which the selected genes display an interrelated expression. “Epithelial Mesenchymal Transition” and “Kras Signaling Up” were activated in both cell lines at 6 h for all treatments (Fig. 2A). At 72 h, Kyse30 cells were characterized by activated “Glycolysis” (EGFhigh and EpEX), “Kras Signaling Up”, “Heme metabolism”, and “P53 Pathway” (EGFhigh), and by suppressed “G2M Checkpoint” and “E2F Targets” (EGFhigh and EGFlow). Suppression of “E2F Targets” and “G2M Checkpoint” was confirmed in FaDu cells treated with high-dose EGF. Additionally, FaDu cells activated MSigDB hallmarks “Inflammatory Response”, “Kras Signaling Up”, and “Epithelial Mesenchymal Transition” following treatment with EGFhigh, EGFlow, and a combination of EGFhigh with EpEX (Fig. 2A).

Fig. 2

EGF- and EpEX-treatment associated GSEA and DEGs. A GSEA were performed for Kyse30 and FaDu cells at 6 h and 72 h after the indicated treatment (EGFlow: EGFL, EGFhigh: EGFH, EpEX, EGFhigh with EpEX: EGF_w_EpEX) using MSigDB hallmark gene sets. Significantly activated or suppressed hallmarks are depicted with gene ratios and adjusted p-values. B Venn diagram of DEGs (|log2FC|> 0.5; adjusted p-value ≤ 0.05) at 6 h and 72 h for Kyse30 and FaDu cells under indicated treatments. C Upset plot displays exclusive DEGs in combinations (connected bullet points). Numbers of exclusive DEGs are indicated

Regulated MSigDB hallmarks were comparable across all four cell treatments at 6 h and identified an induction of EMT by all treatments. However, independent, published results from our group identified morphologic, molecular, and functional traits of EMT only following EGFhigh treatment using morphology features, EMT-TF gene expression patterns, and migration behavior [21]. Therefore, differentially expressed genes (DEGs) of treated cells versus controls were identified to compare the strength of gene regulation across treatments (|log2FC|> 0.5; adjusted p-value ≤ 0.05). EGFhigh treatment resulted in highest numbers of DEGs and EpEX treatment in lowest numbers of DEGs after 6 h (Fig. 2B, left panels). These differences in gene regulation were exacerbated at 72 h, when EGFhigh induced n = 1208 and n = 1536 DEGs in Kyse30 and FaDu cells, respectively (Fig. 2B). All other treatments resulted in no DEGs (Kyse30, EGFlow and EGFhigh + EpEX; FaDu, EpEX), few DEGs (Kyse30, EpEX n = 2; FaDu, EGFlown = 4), or moderate regulation (FaDu, EGFhigh + EpEX n = 103) (Fig. 2B).

To filter genes that are potentially associated with the differential outcome of EpEX, EGFlow, and EGFhigh treatments, DEGs exclusively regulated by each treatment in Kyse30 and FaDu cells were visualized in an upset plot (Fig. 2C). Numbers of exclusive DEGs for each separate treatment at 6 h ranged from n = 1 (FaDu cells treated with EpEX) to a maximum of n = 207 (FaDu cells treated with EGFhigh + EpEX). Highest numbers of exclusive DEGs were observed following treatment with EGFhigh after 72 h in FaDu (n = 919) and Kyse30 cells (n = 709), respectively. Thus, induction of EMT by EGFhigh is associated with sustained and strong gene regulation compared with treatments that failed to induce morphologic and functional EMT.

EpEX induced gene transcription

Effects of EpEX on gene transcription were compared to EGF-dependent regulation at 6 h, since no or only two DEGs were retrieved after 72 h in Kyse30 and FaDu cells, respectively (Fig. 2B). EpEX-Fc induced n = 137 and n = 36 DEGs in Kyse30 and FaDu cells, respectively (Fig. 3A). In Kyse30 cells, n = 73/137 and n = 77/137 EpEX DEGs were shared with EGFlow and EGFhigh, respectively. Shared DEGs were similarly regulated with Spearman ρ-values of 0.925 and 0.942 and high statistical significance (Fig. 3B). In FaDu cells, n = 33/36 and n = 34/36 EpEX DEGs were shared with EGFlow and EGFhigh, respectively, and ρ-values were 0.884 and 0.854 (Fig. 3C). Hence, none of the shared DEGs were counter-regulated or differed substantially regarding the magnitude of regulation. In Kyse30 cells, EpEX induced n = 64 and n = 60 exclusive DEGs compared to EGFlow and EGFhigh, respectively. Exclusive EpEX DEGs were very restricted in FaDu cells (n = 3 and n = 2). Furthermore, EGFlow and EGFhigh treatments induced a substantially higher number of exclusive DEGs compared to EpEX-Fc in both cell lines (Kyse30: EGFlown = 329, EGFhighn = 540; FaDu: EGFlown = 763, EGFhighn = 995) (Fig. 3B-C).

Fig. 3

EpEX regulates gene expression and competes with EGF for binding to EGFR. A Volcano plots of DEGs in Kyse30 and FaDu cells 6 h after EpEX treatment. Numbers of DEGs for each cell line are indicated (|log2FC|> 0.5; adjusted p-value ≤ 0.05). B-C Shared and exclusive DEGs between EpEX, EGFlow, and EGFhigh are depicted in scatter/bar plots with Spearman correlation and p-values. D Schematic representation of the EGF/EpEX competition assay. Cells were incubated with fluorescence-labeled EGF alone or in combination with non-labeled EGF or EpEX in equimolar conditions. E–F Fluorescence of Kyse30 and FaDu cells incubated with fluorescence-labeled EGF and combinations with non-labeled EGF and EpEX was quantified by flow cytometry. E Representative histograms of EGFR staining in a competition assay in Kyse30 and FaDu cells. F Mean values with SD of n = 3 independent experiments are shown. P-values from one-way ANOVA with Tukey´s multiple comparison test. * ≤ 0.05; ** ≤ 0.01; *** ≤ 0.001; **** ≤ 0.0001. Control = unstained cells; EGF-Fluo: fluorescence-labeled EGF; EGF-Fluo + EGF: fluorescence-labeled EGF plus unlabeled EGF; EGF-Fluo + EpEX: fluorescence-labeled EGF plus EpEX

We concluded that EpEX-Fc is a weak ligand of EGFR that regulates a subset of EGF-dependent genes and few unique genes. Thus, EpEX-Fc-dependent repression of EGFR-mediated EMT is most likely not the result of a transcriptional repression but of a competition with EGF for binding to EGFR. To test this hypothesis, Kyse30 and FaDu cells were incubated with fluorescently labeled EGF in absence or presence of equimolar amounts of unlabeled EGF or EpEX-Fc, and fluorescence intensities were measured by flow cytometry (Fig. 3D). Incubation with labeled EGF increased mean fluorescence intensities (MFI) in both cell lines (38.8-fold and 7.0-fold in Kyse30 and FaDu, respectively). Differences in fluorescence induction were congruent with EGFR expression levels on the respective cell line. Co-treatment with unlabeled EGF reduced the fluorescence by 66.7% and 47.2% in Kyse30 and FaDu cells, respectively. Co-treatment with EpEX-Fc reduced the fluorescence by 47.1% and 28.8% in Kyse30 and FaDu cells, thus confirming a competitive binding of EGF and EpEX to EGFR (Fig. 3E-F). Hence, EpEX is a ligand of EGFR that induces lesser transcriptional changes than EGF and competes with EGF for binding to EGFR.

EGFR-mediated EMT differs from EMT hallmark and pEMT signatures

Intersected exclusive DEGs induced by EGFhigh in Kyse30 and FaDu cells were extracted (n = 181, Fig. 2C) and 171 genes were similarly regulated in both cell lines. This EGFR-mediated EMT signature (n = 171) was subjected to an over-representation analysis (ORA). GO-terms activated in the EGFR-mediated EMT signature were related to cell migration and invasion (“focal adhesion”, “cell-substrate junction”, “cell leading edge”, and “cadherin binding”) (Fig. 4A, left panel). Suppressed GO-terms were related to DNA replication and cell division, pinpointing a reduced cell proliferation following EGFR-mediated EMT (Fig. 4A, middle panel). This finding was corroborated by the suppressed KEGG terms “cell cycle” and “DNA replication”. Additionally, EGFR-mediated EMT was associated with the KEGG terms of “cellular senescence” and “p53 signaling pathway” (Fig. 4A, right panel). Genes included in enriched GO-terms were extracted and linkages with GO-terms are depicted in a gene-concept network. HNSCC cancer stem cell (CSC) marker CD44, integrin beta 4 (ITGB4), and small GTPase Rac2 were linked to the terms “Cell leading edge”, “Cell-substrate junction”, and “Focal adhesion”, while further genes showed linkages to one or two GO-terms (Fig. 4B). Gene-concept networks of suppressed GO-terms and KEGG-terms are depicted in chord plots in Supplementary Fig. 3. Regulation of Polo-like kinase 1 (PLK1) and Kinesin family members was linked to numerous GO-terms related to cell cycle and nuclear division.

Fig. 4

EGFR-mediated EMT-dependent 5-gene prognostic signature for HNSCC. A Genes of the EGFR-mediated EMT signature (n = 171) were subjected to an over-representation analysis. Significantly activated and suppressed pathways in GO and KEGG are depicted with gene counts and adjusted p-value. B Genes from the enriched GO-terms are depicted in a gene-concept network. C Forest plot of the multivariable Cox PH regression model incorporating the 5-gene signature in the training cohort (TCGA) of n = 240 HPV-negative HNSCC with event numbers, log-rank p-value, AIC, and concordance index. D Gene expression of DDIT4, FADD, ITGB4, NCEH1, and TIMP1 in normal mucosa (n = 34) and HNSCC (n = 238) of the training cohort (TCGA); p-values *** ≤ 0.001; **** ≤ 0.0001. E Principal component analysis and hierarchical clustering of the top 25% most strongly regulated DEGs for matched pairs of normal mucosa and HNSCC of the TCGA cohort are shown (n = 34). F Volcano plot of DEGs from matched pairs of normal mucosa and HNSCC of the TCGA cohort (n = 34; |log2FC|> 0.5, p-value ≤ 0.05). Genes of the 5-gene signature are significantly up-regulated DEGs. G Stratification of HPV-negative HNSCC of the training cohort (TCGA; n = 240) with a 5-gene signature-based risk score (median cut-off) for overall survival (time in months). Numbers at risk, HR, 95% CI, and p-value are indicated in the Kaplan–Meier survival curve. H Validation of the risk score in the MD Anderson Cancer Center (MDACC; n = 74) and the Fred Hutchinson Cancer Research Center (FHCRC; n = 97) HNSCC cohorts. Median cut-off threshold of the training cohort served to dichotomize the MDACC and FHCRC cohorts (risk − , risk +) for overall survival. Numbers at risk, HR, 95% CI, and p-value are indicated in the Kaplan–Meier survival curve

Analysis of overlapping and unique DEGs of the EGFR-mediated EMT signature with the MSigDB EMT hallmark and the pEMT [10] gene sets revealed only few overlapping genes (n = 7/171 (4%) and n = 4/171 (2.3%), respectively), whereas pEMT and EMT hallmark showed more similarity (n = 36/100 (36%)) (Supplementary Fig. 4A and Supplementary Table 1). These results suggested molecular heterogeneity across all three signatures potentially defining differing states of EMT. Two publicly available datasets (Gene Expression Omnibus accession numbers GSE23952 and GSE32254), which were used by Broad MSigDB for refining and validating their ‘EMT hallmark’ gene set, were used here as ground-truth to define EMT. Gene set variation analysis (GSVA) was conducted with the MSigDB hallmark gene set (n = 200 genes), the pEMT signature (n = 100 genes) [10], and the EGFR-mediated EMT signature (n = 171 genes). All three EMT signatures were significantly enriched in cells that had undergone EMT following TGFβ (GSE23952) and TNFα/TGFβ treatment (GSE32254) (Supplementary Fig. 4B).

Malignant cells (n = 2,176) from the scRNA-seq HNSCC dataset GSE103322 were subjected to a GSVA of the pEMT signature, the MSigDB hallmark gene set, and the EGFR-mediated EMT signature. Enrichment of the pEMT signature across all ten HNSCC patients was regarded as ground-truth in this analysis and GSVA results were highly comparable to original findings, which were re-computed in analogy to Puram et al. using Seurat R AddModuleScore (Supplementary Fig. 4C-D, Spearman ρ = 0.95, p-value 2.2e-16). Comparison of pEMT, EMT, and EGFR-mediated EMT signature enrichments at the single cell level showed overlap as well as differences across all three signatures. Cells enriched for the EGFR-EMT were more frequent across patients and were also detected in HNSCC that have not been deemed in pEMT (P6, P20). Furthermore, single cells within patients differed in the enrichment of EGFR-mediated EMT, whereas pEMT and EMT enrichment were more concurrent (see for example P5, P16, P22, and P28) (Supplementary Fig. 4C). Correlation analysis of all three EMT scores confirmed a higher concordance between pEMT and the MSigDB EMT hallmark and lower concordance with EGFR-mediated EMT (Supplementary Fig. 4D). Lastly, comparing pEMT, EMT, and EGFR-mediated EMT scores at the individual patient level confirmed a more homogeneous distribution of EGFR-mediated EMT across patients. However, a more heterogeneous degree of EGFR-mediated EMT was observed across single cells of a given patient, with two or three separate major sub-populations of single cells defined in violin plots (Supplementary Fig. 4E). Hence, EGFR-mediated EMT defines a meta-program in single malignant cells of HNSCC that is partly redundant yet distinct to pEMT and EMT.

EGFR-mediated EMT-dependent prognostic risk score in HNSCC

The EGFR-mediated EMT signature was transferred to HPV-negative patients of the TCGA-HNSCC cohort [5] (n = 240; Supplementary Table 2) as training cohort to develop a prognostic risk score. Expression values for n = 170 genes from the EGFR-mediated EMT signature were available and subjected to univariate Cox regression analysis of the correlation to OS. Assuming that EGFR-mediated EMT is detrimental to the patients´ survival, up-regulated genes with a HR > 1 and down-regulated genes with HR < 1 were filtered (n = 53 genes). Following forward feature selection using rbSurv, a multivariate Cox regression model defined a 5-gene signature composed of NCEH1 (Neutral Cholesterol Ester Hydrolase 1), DDIT4 (DNA Damage Inducible Transcript 4), ITGB4, FADD (Fas-Associated protein with Death Domain), and TIMP1 (Tissue Inhibitor of Metaloproteinases 1) as prognostic risk score (Fig. 4C).

Gene expression of the 5-gene signature was assessed with HTSeq-counts for n = 238 HNSCC and n = 34 normal samples from TCGA. The expression of all five genes was induced following EGFhigh-treatment in Kyse30 and FaDu cells and were significantly up-regulated in HNSCC versus normal samples (Fig. 4D). Thirty-four matched pairs of normal tissue and HNSCC were identified within TCGA, which separated in two major clusters in a PCA and a hierarchical clustering heatmap based on the top 25% most highly expressed genes (Fig. 4E). All five genes of our prognostic signature (NCEH1, DDIT4, ITGB4, FADD, and TIMP1) were significantly up-regulated DEGs from tumors in matched pairs of HNSCC and normal tissue from TCGA data (Fig. 4F; |log2FC|> 0.5, p-value ≤ 0.05).

HPV-negative TCGA-HNSCC patients were stratified according to the 5-gene signature-based risk score (median threshold) and survival probabilities depicted in Kaplan–Meier (KM) curves. High risk scores correlated with significantly reduced 5-year OS (Fig. 4G; HR 2.41, 95% CI: 1.64–3.55, p-value = 4.38e-06). mRNA coefficients and the median threshold of the prognostic model trained in the TCGA cohort were transferred to validation cohorts (MDACC-HNSCC, n = 74; FHCRC, n = 97; Supplementary Table 2) and confirmed the correlation of the EGFR-mediated EMT-based risk score with poor OS (Fig. 4H; MDACC: HR 3.92, 95% CI: 1.9–8.11, p-value = 6.13e-05; FHCRC: HR 3.32, 95% CI: 1.73–6.38, p-value = 8.89e-05). Thus, the EGFR-mediated EMT-based risk score prognosticates OS in HNSCC.

The prognostic value of the EGFR-mediated EMT signature to predict OS was compared to published EMT signatures for HNSCC. The MSigDB EMT hallmark signature, the HNSCC pEMT signature by Puram et al. and the HNSCC EMT signatures by Jung et al. and Vallina et al. served as comparisons [10, 29, 30]. Except for the Vallina et al. signature that is composed of six genes, which were extracted in a meta-analysis of eight independent cohorts, all remaining larger signatures were subjected to feature selection using rbsurv and multivariate Cox regression model within the TCGA HNSCC cohort. The MSigDB EMT hallmark signature (n = 200 genes) retrieved a 4-gene risk score (NT5E, PVR, COL8A2, and APLP1), the pEMT signature (n = 100 genes) a 5-gene risk score (IGFBP7, EMP3, SLC39A14, CALU, and SLC38A5), and the Jung et al. signature (n = 82 genes) a 2-gene risk score (GLT8D2 and COL6A1). Only the pEMT- and the EGFR-mediated EMT-based risk scores showed a significant Log-Rank p-value in the multivariate Cox model, with a concordance index of the pEMT-based risk score that was slightly inferior to the EGFR-mediated EMT-based risk score (0.59 vs. 0.63) (Supplementary Fig. 5A-D and Fig. 4C). Stratification of the patients according to the median risk score for each model showed a significant prognosis of OS for the EGFR-mediated EMT-, the MSigDB EMT hallmark- and the pEMT-based risk scores in KM curves (EMT hallmark: HR 1.46; 95% CI 1.01–2.11; p = 0.044; pEMT: HR 1.69; 95% CI 1.16–2.46; p = 0.0054) (Supplementary Fig. 5A-D and Fig. 4C). Analysis of the area under the curve (AUC) for false and true positives for risk scores generated from all five models confirmed superior performances of the EGFR-mediated EMT- and the pEMT-based risk scores to prognosticate 3- and 5-year OS (Supplementary Fig. 5E).

ITGB4-associated gene expression in HNSCC

NCEH1, DDIT4, ITGB4, FADD, and TIMP1 expressions in the TCGA cohort were correlated with major signaling pathway activities inferred using PROGENy (Pathway RespOnsive GENes for activity inference). PROGENy compiles ground-truth validated core responsive genes from 14 major cellular pathways, allowing to compute pathway activity scores from bulk and scRNA-seq datasets [31]. ITGB4 was significantly correlated with EGFR and mitogen-activated protein kinase (MAPK) pathway activities (Fig. 5A-B), which are both required for EGFR-mediated EMT [21]. ITGB4 co-regulated genes were identified by batch correlation analysis in the HPV-negative TCGA cohort. The highest Spearman correlation was observed for Integrin alpha 3 (ITGA3), which is a binding partner of integrin beta 1 (ITGB1) that forms a heterodimeric laminin receptor. LAMA3, LAMB3, and LAMC2, together encoding the ITGB4 ligand laminin 5, were in the top ten co-regulated genes and are part of the common pEMT gene signature [10] that prognosticates OS in HNSCC [14]. Focal adhesion-associated adapter protein paxillin was likewise strongly correlated to ITGB4 expression. Top ten genes that showed a counter-regulated expression to ITGB4 were more heterogeneous in function, including transcriptional activators/repressors and replication factors (CREG1, REPIN1, AFF3, WDSUB1, EID3), membrane-associated enzymes and functional proteins (CA5B, GULP1), methyltransferase 9, complement inhibitor SUSD4, and one gene of unclear function (C12ORF57) (Supplementary Table 3).

Fig. 5

Integrin beta 4 (ITGB4) expression in HNSCC. A Pathway activities were inferred in HPV-negative TCGA-HNSCC samples using PROGENy (n = 240). Spearman correlations and p-values of the 5-gene signature with pathway activities are depicted. * ≤ 0.05; ** ≤ 0.01; *** ≤ 0.001. B Scatter plots of ITGB4 correlations with EGFR and MAPK pathways in HPV-negative HNSCC of TCGA patients (n = 240) are shown with Spearman correlation and p-value. C ITGB4-correlated genes were subjected to a GSEA using MSigDB hallmark gene sets. Top 15 regulated hallmarks are depicted with gene counts and adjusted p-values. D ITGB4 expression was compared in matched pairs of normal mucosa and HPV-negative HNSCCs of TCGA (n = 34). Matched expression values are shown in a box plot whiskers graph (t-test **** ≤ 0.0001). E Representative examples of ITGB4 expression in normal mucosa (n = 64), primary tumor (n = 80) belonging to n = 84 patients, and in triplets including lymph node metastases (n = 26) from HPV-neg. and -pos. HNSCC of the LMU cohort. F-G IHC scoring of ITGB4 protein expression is shown in scatter dot plots with means and SD for all samples and stratified according to HPV-status. Ns: not significant, **** ≤ 0.0001 (t-test and One-way ANOVA). H Representative examples of ITGB4 and laminin 5 expression in consecutive sections of normal mucosa and HNSCC are shown. I Double immunofluorescence staining of ITGB4 (red), laminin 5 (green), and DAPI (blue) in HNSCC with edge-localized (left) or homogeneous ITGB4 expression (right)

ITGB4-associated hallmarks were inferred via GSEA of all genes from TCGA-HNSCC ranked by their Spearman correlation with ITGB4 and using MSigDB hallmarks gene sets. Activated hallmarks were related to cell proliferation (“Hallmark of E2F targets”, “Hallmark of G2M checkpoint”, Hallmark mitotic spindle”), glycolysis, apical junction, TNF alpha signaling via NF kappa B, and EMT. Suppressed hallmarks referred to bile acid metabolism, oxidative phosphorylation, and fatty acid metabolism (Fig. 5C).

ITGB4 and its ligand laminin 5 are up-regulated in HNSCC

Matched pairs of normal mucosa and HNSCC from the TCGA-HNSCC cohort confirmed a significant over-expression of ITGB4 mRNA in tumor samples (Fig. 5D). ITGB4 protein expression was analyzed in an in-house cohort of HPV-negative and -positive HNSCC patients (Supplementary Table 2). Primary HNSCCs (n = 80) in comparison to normal mucosa (n = 64) proved a significant over-expression of ITGB4 in tumors (Fig. 5E and F left panel). Representative IHC staining in Fig. 5E showed a strong expression of ITGB4 at the basal cell membrane adjacent to the basal lamina and in suprabasal cell layers of normal mucosa and the over-expression in tumor areas. Matched triplets (n = 26) of normal mucosa, primary tumor, and lymph node metastases demonstrated a significant over-expression of ITGB4 in primary tumors and lymph node metastases (Fig. 5E and F right panel). Stratification of patients according to the HPV-status showed a significant up-regulation of ITGB4 in HPV-negative (n = 43) and -positive (n = 41) patients (Fig. 5G).

ITGB4 and laminin 5 were assessed in serial sections of HNSCC with weak and strong ITGB4 expression. In normal mucosa, ITGB4 and laminin 5 were co-localized at the basal lamina and ITGB4 was additionally expressed in suprabasal cell layers (Fig. 5H). In ITGB4-negative or ITGB4low HNSCC, laminin 5 was either absent or only expressed in non-malignant cells representing endothelial cells and leukocytes. Differing localization of ITGB4 was observed within tumor areas including an exclusive expression at the interface between tumor and non-malignant stromal tissue, a marginal expression at the edges of tumor areas, and a more homogeneous expression throughout tumor cells. Laminin 5 co-localized with ITGB4 at the tumor-stroma interface (Fig. 5H). Co-localization patterns of ITGB4 and laminin 5 were confirmed in double immunofluorescence staining (Fig. 5I).

Normal mucosa (n = 34), and HNSCC (n = 238) and matched pairs of normal mucosa and HNSCC (n = 34) from the TCGA cohort were assessed for the expression of integrin alpha 6 (ITGA6) and laminin 5 genes LAMA3, LAMB3, and LAMC2. Both sub-cohorts demonstrated a significant up-regulation of ITGA6, LAMA3, LAMB3, and LAMC2 in tumors (Supplementary Fig. 6A). ITGB4 and ITGA6, LAMA3, LAMB3, and LAMC2 showed a positive correlation in HNSCC of the TCGA cohort and in single malignant HNSCC cells (Supplementary Fig. 6B-C).

ITGB4 correlates with EGFR activity in malignant single HNSCC cells

The 5-gene signature expression was analyzed in malignant single HNSCC cell data from Puram et al. [10]. Ten oral cavity carcinomas with a total of n = 2176 single cell transcriptomes from the scRNA-seq dataset GSE103322 were included in the analysis (tSNE plot in Fig. 6A). Expression analysis of the 5-gene signature revealed highest expressions of ITGB4 and DDIT4 at the individual patient level and highest percentages of positive single cells (Fig. 6A and B).

Fig. 6

5-gene signature expression in malignant single HNSCC cells. A T-distributed stochastic neighbor embedding (t-SNE) plots of malignant single HNSCC cells (n = 2176; GSE103322; left plot). Expression of NCEH1, DDIT4, ITGB4, FADD, and TIMP1 are displayed in t-SNE plots in n = 2176 malignant single cells (right plots). B 5-gene signature expression is represented for individual patients with percent of malignant single cell expression and average expression values. C Pathway activities of malignant single cells from GSE103322 were inferred using PROGENy. Spearman correlations and p-values of the 5-gene signature with pathway activities are depicted. * ≤ 0.05; ** ≤ 0.01; *** ≤ 0.001. D Scatter plots of ITGB4 correlations with the EGFR pathway in n = 2176 malignant single cells from GSE103322 are shown with Spearman correlation and p-value. E ITGB4 cell surface protein expression is shown in Kyse30 and FaDu cells treated with 50 ng/mL EGF for 72 h. Mean values with SD are shown in scatter dot plots of n = 3 independent experiments. ** ≤ 0.01 (t-test). F ITGB4 mRNA expression in Kyse30 and FaDu cells treated with 50 ng/mL EGF or combinations of EGF with Cetuximab, MEK inhibitor, or AK inhibitor. Mean values with SD of qPCR measurements are shown in scatter dot plots after 72 h of the indicated treatment of n = 3 independent experiments. * ≤ 0.05; ** ≤ 0.01; *** ≤ 0.001 (one-way ANOVA)

NCEH1, DDIT4, ITGB4, FADD, and TIMP1 expression was correlated to pathway activities inferred in single cells using PROGENy. Highest significant correlation was observed for ITGB4 and the EGFR pathway activity (Fig. 6C-D). Additionally, ITGB4 expression was analyzed in scRNA-seq datasets of different cancer entities using TISCH (Tumor Immune Single-Cell Hub) [32]. TISCH is an online scRNA-seq database composed of n = 79 datasets and n = 2,045,746 cells, allowing for the expression analysis of genes of interest in malignant and non-malignant cells of n = 18 individual cancer types including HNSCC (http://tisch.comp-genomics.org/). Cancer entities with highest ITGB4 expression in single malignant cells were compared (log(TPM/10 + 1) > 0.5). ITGB4 was substantially expressed in single cells of basal cell carcinoma, cholangiocarcinoma (CHOL), colorectal carcinoma CCRC), glioma, HNSCC, non-small cell lung cancer (NSCLC), ovarian serous cystadenocarcinoma (OV), and pancreatic adenocarcinoma (PAAD). HNSCC were characterized by the highest ITGB4 expression in malignant single cells and the strongest differential expression compared to immune and stromal cells (Supplementary Fig. 7A, upper panel). Expression analyses in subtypes of immune and stromal cells showed additional expression of ITGB4 in endothelial cells, only (Supplementary Fig. 7B). ITGB4 expression was represented in violin plots for the three tumor entities with highest expression, i.e. HNSCC, PAAD, and CRC, confirming the strong and most selective expression in malignant HNSCC cells (Supplementary Fig. 8C). ITGA6, LAMA3, LAMB3, and LAMC2 followed the expression pattern of ITGB4 and was highest in malignant single HNSCC cells und most differential to immune, stromal, and other cells (Supplementary Figs. 8– 11).

Induction of ITGB4 cell surface expression by EGFR activation was validated using flow cytometry. ITGB4 protein expression at the plasma membrane was enhanced 2.53-fold and 4.62-fold in Kyse30 and FaDu cells upon treatment with EGFhigh, respectively (Fig. 6E). ITGB4 mRNA expression was induced by 5.12-fold and 1.99-fold in