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).
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).
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, EGFlow n = 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: EGFlow n = 329, EGFhigh n = 540; FaDu: EGFlow n = 763, EGFhigh n = 995) (Fig. 3B-C).
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.
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).
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).
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 Kyse30 and FaDu cells upon EGFhigh treatment. Cetuximab and MEK inhibitor co-treatment blocked this induction, whereas AKT inhibition had no impact on ITGB4 mRNA induction (Fig. 6F). Hence, ITGB4 is over-expressed and associated with enhanced EGFR and MAPK activity in malignant single cells of HNSCC.
ITGB4 promotes migration and invasion
ITGB4 expression was knocked down (KD) in Kyse30 and FaDu cells using specific and scrambled shRNA in lentiviral vectors. KD cells generated with different MOI for each cell line were further analyzed (n = 2). Wildtype and scrambled shRNA controls (Neg) expressed similar levels of ITGB4, whereas ITGB4-KD resulted in 84% and 87% reduction in Kyse30 cells and in 73% and 68% reduction in FaDu cells (Supplementary Fig. 12). Migration and invasion of ITGB4-KD cells was tested in Boyden chambers in the absence or presence of Matrigel. Treatment of Kyse30 and FaDu cells with EGFhigh promoted migration and invasion. ITGB4-KD induced a substantial reduction or abrogation of EGFhigh-mediated migration and invasion. Cetuximab inhibited migration and invasion, thus confirming the specificity for EGFR activity (Fig. 7A-B). EGFhigh treatment also increased wound closure in Kyse30 and FaDu control cells, whereas ITGB4-KD cells migrated only weakly or not at all. Co-treatment of all cell lines with Cetuximab entirely abrogated migration as compared to negative controls (Fig. 7C and Supplementary Fig. 13).
EGFR-mediated local invasion was addressed in a 3D model. Spheroids of FaDu control and ITGB4-KD cells were transferred into Matrigel before addition of EMT-inducing concentrations of EGF (EGFhigh). Treatment of control cells induced a strong invasion into the ECM, which was efficiently inhibited by Cetuximab (Fig. 7D). ITGB4-KD cells showed a reduction of invasion, which was further inhibited by Cetuximab (Fig. 7D). Quantification of invasive area and migratory distance confirmed a strong induction of EGFR-mediated local invasion that was significantly reduced in ITGB4-KD cells. It must further be noted that cell numbers in the invasive area were also considerably diminished as visualized in Fig. 7D. Invasive area and distance were efficiently inhibited by Cetuximab (Fig. 7E). Hence, ITGB4 expression supports EGFR-mediated migration and local invasion.
ITGB4 and the proliferation marker Ki67 were visualized by immunofluorescence staining and confocal imaging in spheroids maintained in Matrigel. Owing to the staining and imaging method, we concentrated on cells of the outer layers to address ITGB4 availability. Control-treated spheroids showed a weak expression of ITGB4, where only the most outer cells expressed higher levels of the protein in areas of contact to Matrigel. Cells at the leading edges of invasion following EGF treatment expressed high levels of ITGB4 and reduced levels of proliferation marker Ki67 (Fig. 7F). Hence, ITGB4 is strongly expressed and available as a target in locally invasive cells.
In primary HNSCC, initial steps of local invasion were assessed as tumor budding, which is defined as single cells or clusters composed of less than 5–10 cells. Budding was observed in primary HNSCC with a homogeneous and an edge localization of ITGB4-positive tumor cells (Fig. 7G). Budding was quantified as negative, weak, intermediate, and strong according to the frequency of tumor buds (Supplementary Fig. 14). Edge localization of ITGB4 was observed in 31.1% of cases (33/106) (Fig. 7H) and was associated with significantly more tumor budding (Fig. 7I) and higher budding intensity than the homogeneous ITGB4 localization (Fig. 7J).
Antagonizing ITGB4 inhibits local invasion.
ITGB4-antagonizing antibody ASC8 was tested in 2D and 3D models. EGFhigh treatment induced migration and invasion in a Boyden chamber in the absence or presence of Matrigel, which were both blocked upon co-treatment with Cetuximab. Co-treatment with ASC8 antibody had no significant effect on cell migration but blocked invasion comparably to Cetuximab (Fig. 8A-B).
Effects of ASC8 antibody on local invasion were further addressed in spheroids. Wildtype FaDu spheroids embedded in Matrigel were treated with EGFhigh, EGFhigh plus Cetuximab, EGFhigh plus ASC8, and EGFhigh plus an anti-GFP antibody as a control. EGFhigh treatment induced local invasion of FaDu cells into Matrigel, which was inhibited by Cetuximab treatment. Antagonizing ITGB4 antibody resulted in strong inhibition of invasion, while the control anti-GFP antibody had no effect (Fig. 8C). Quantification of the invasive area (IA) and invasive distance (ID) confirmed a significant and strong increase following treatment with EGFhigh, which was efficiently inhibited by Cetuximab (IA: 76.19% reduction; ID: 37.9%). ASC8 antibody treatment decreased the IA by 40.5% and the ID by 21%, anti-GFP antibody had no significant effect (IA: 6.5%; ID: 2.17%) (Fig. 8D).
EGFR-mediated EMT and response to Cetuximab
EGFR-mediated EMT may have various repercussions on tumor progression. EMT is generally associated with resistance to Cetuximab [33] and EGFR-mediated EMT promotes tumor invasion. To address potential roles of EGFR-mediated EMT in the regulation of treatment response, unique and overlapping DEGs (|log2FC|> 0.5, p-value < 0.05) following Cetuximab treatment in the sensitive HNSCC cell lines SCC1, SCC6 and SCC25 (GSE137524) were determined in the EGFR-mediated EMT signature. Overlapping DEGs were defined as genes with identical directionality of regulation (up- or down-regulated) in SCC cell lines following Cetuximab treatment and the EGFR-mediated EMT signature. A total of 48, 20, and 19 overlapping DEGs with EGFR-mediated EMT (n = 171 genes) were determined for SCC1, SCC6, and SCC25, respectively (Supplementary Table 4).
Overlapping DEGs from SCC1, SCC6, and SCC25 were separately assigned to hallmark gene sets of MSigDB using the functions “Investigate Gene Sets” and “Compute Overlaps with Hallmarks”. Overlapping DEGs were significantly enriched in “HALLMARK_E2F-TARGETS”, “HALLMARK_G2M-CHECKPOINT”, and “HALLMARK_MITOTIC_SPINDLE”. All overlapping DEGs contributing to significantly enriched hallmarks involved in cell cycle regulation were down-regulated genes. Intersection of overlapping DEGs from SCC1, SCC6, and SCC25 (MND1, MNS1, CCNB1, HMGB2, PLK1, and CDC20) represents common DEGs between Cetuximab-treated SCC cell lines and the EGFR-mediated EMT. These DEGS confirmed cell cycle down-regulation within “HALLMARK_E2F_TARGETS” and “HALLMARK_G2M_CHECKPOINT”. Hence, genes significantly regulated early upon Cetuximab treatment and in EGFR-mediated EMT suggested a common feature of inhibition of cell cycle progression, which may improve resistance to therapeutic treatment in general and, specifically, upon Cetuximab treatment.
Furthermore, overlapping DEGs between the EGFR-mediated EMT and Cetuximab resistant versus sensitive cells were extracted from the GSE21483 dataset (Cetuximab-sensitive SCC1 cells and their resistant derivative 1Cc8). A total of 22 overlapping DEGs were identified (Supplementary Table 4), however the assignment of these DEGs to MSiGDB gene sets failed to identify significantly enriched hallmarks.
Thus, we conclude that DEGs common to the EGFR-mediated EMT signature and Cetuximab-treated cells are rather related to an early response to treatment and the inhibition of cell cycle progression. Within an established resistance to Cetuximab, DEGs overlapping between EGFR-mediated EMT and Cetuximab resistance are more heterogeneous and did not contribute to any significant enrichment of MSigDB hallmarks.
Next, we concentrated on the potential involvement of genes of the EGFR-mediated EMT signature in functional aspects of tumor progression. As part of a prognostic 5-gene signature of EGFR-mediated EMT, ITGB4 is functional in invasion and its up-regulation is blocked by Cetuximab. Treatment with Cetuximab remains a central tool in multimodal palliative management of progressed metastatic HNSCC, but markers indicative of response to therapy are lacking. Bossi et al. (2016) have reported a whole-genome cDNA analysis of a clinical cohort of recurrent-metastatic HNSCC (RM-HNSCC; n = 40) that were all treated with chemotherapy in combination with Cetuximab. N = 14 patients were characterized by long PFS (median = 19 months) and n = 26 patients showed a short PFS (median = 3 months [34]). Assuming ITGB4 might represent a surrogate marker for EGFR-mediated EMT, its association with PFS was assessed. The odd ratio of ITGB4 expression for PFS estimation was adjusted by available demographic and clinical features in a multivariate logistic regression model. Low ITGB4 expression (median cut-off) was significantly associated with higher relative risk of shortened PFS (HR: 6.95, 95% CI 1.25–55.68, p-value 0.04) (Fig. 8E). These results suggest that patients with higher ITGB4 expression profited more from Cetuximab treatment.