IL-10 and integrin signaling pathways are associated with head and neck cancer progression
© Bornstein et al. 2016
Received: 30 April 2015
Accepted: 28 December 2015
Published: 8 January 2016
Head and neck cancer is morbid with a poor prognosis that has not significantly improved in the past several decades. The purpose of this study was to identify biological pathways underlying progressive head and neck cancer to inform prognostic and adjuvant strategies. We identified 235 head and neck cancer patients in The Cancer Genome Atlas (TCGA) with sufficient clinical annotation regarding therapeutic treatment and disease progression to identify progressors and non-progressors. We compared primary tumor gene expression and mutational status between these two groups.
105 genes were differentially expressed between progressors and nonprogressors (FDR < 0.05). Pathway analyses revealed deregulation (FDR < 0.05) of multiple pathways related to integrin signaling as well as IL-10 signaling. A number of genes were uniquely mutated in the progressor cohort including increased frequency of truncating mutations in CTCF (P = 0.007). An 11-gene signature derived from a combination of unique mutations and differential expression was identified (PAGE4, SMTNL1, VTN, CA5A, C1orf43, KRTAP19-1, LEP, HRH4, PAGE5, SEZ6L, CREB3). This signature was associated with decreased overall survival (Logrank Test; P = 0.03443). Cox modeling of both key clinical features and the signature was significant (P = 0.032) with the greatest prognostic improvement seen in the model based on nodal extracapsular spread and alcohol use alone (P = 0.004).
Molecular analyses of head and neck cancer tumors that progressed despite treatment, identified IL-10 and integrin pathways to be strongly associated with cancer progression. In addition, we identified an 11-gene signature with implications for patient prognostication. Mutational analysis highlighted a potential role for CTCF, a crucial regulator of long-range chromatin interactions, in head and neck cancer progression.
The 5-year survival rate for primary locally advanced head and neck squamous cell carcinoma (HNSCC) is approximately 50 % , however recurrent disease carries a dismal prognosis of 10.1 months with first line chemotherapy . HNSCC recurs ~30 % of the time, most often within the first 1–2 years of definitive treatment. Pathways associated with progression have been identified using array-based gene expression analysis; however these studies are limited by the lack of rigor using older analysis techniques and normalization techniques, and heterogeneously treated patients. Identification of specific pathways linked to progression after radiation has the promise of informing targeted strategies to improve the prognosis of head and neck cancer.
The Cancer Genome Atlas (TCGA) is a joint effort of the National Cancer Institute and the National Human Genome Research Institute and has revolutionized the ability of investigators to ask prognostic questions about tumor biology previously limited by suboptimal sample sizes, non-standard sample normalization, and outdated techniques. HNSCC lends itself to this study, as there are diverse epidemiologic risk factors (e.g. HPV vs. tobacco), anatomic subsites (e.g. larynx vs. oral cavity), and issues with heterogeneity even within subpopulations of tumor cells within a given patient. To this end, TCGA has generated whole exome sequencing, SNP array, DNA methylation, RNA-Seq, and miRNA-Seq data for a large collaborative cohort of HNSCC patients. Their results have recently been published highlighting distinct subgroups within newly diagnosed HNSCC patients e.g. different mutational profiles between HPV-driven and tobacco-related tumors . TCGA annotated data on these patients provides an unprecedented opportunity to determine which molecular pathways are most associated with disease progression and survival, in order to gain insight into potentially targetable biology. We sought to determine what molecular alterations were unique among HNSCC progressors in TCGA to help inform future patient stratification and adjuvant treatment.
Results and discussion
Demographics for TCGA HNSCC patients analyzed in this study (progressors and nonprogressors)
Patients (n = 235)
173(74 %)/62(26 %)
212(90 %)/12(4.5 %)/4(2 %)/2(1 %)/5(2.5 %)
181(77 %)/52(22 %)/2(1 %)
169(71.5 %)/61(26 %)/5(2.5 %)
HPV p16 or ISH (+/−/NA)
19(8 %)/49(21 %)/167(71 %)
135(57 %)/41(18 %)/59(25 %)
T Stage (T1-T2/T3/T4/TX/NA)
80(34 %)/64(27 %)/87(37 %)/3(1.5 %)/1(0.5 %)
N Stage (N0/N+/NA)
107(45.5 %)/127(54 %)/1(0.5 %)
Tumor Stage (I-III/IV/NA)
96(41 %)/135(57 %)/4(2 %)
Margin Status (+/−/Close/NA)
17(7 %)/158(67 %)/22(10 %)/38(16 %)
Nodal Extracapsular Spread (GE/ME/NE/NA)c
13(6 %)/34(14 %)/114(49 %)/74(31 %)
Curated Therapy: Therapy (C/R/CR/CRTM/CRTMV/NA)d
5(2 %)/60(26 %)/95(40 %)/1(0.5 %)/1(0.5 %)/73(31 %)
Radiation data: Radiation Dose cGy (mean)
Follow-up data: Follow-up Days (median)
Follow-up data: Mortality (Living/Deceased)
164(69.5 %)/71(30.5 %)
Follow-up data: Days to Death (median)
Follow-up data: Days to New Tumor (median)
Radiation Treatment (Y/N)
165(70 %)/70(30 %)
68(29 %)/167(71 %)
Differentially expressed genes
Radiation treatment assignment
We examined differential expression in progressors and nonprogressors assigned to radiation treatment to determine the extent of overlap in enriched pathways between this subset and the entire group. There were 60 progressors assigned radiation treatment and 105 nonprogressors. From the 492 differentially expressed genes (based on FDR < 0.05, data not shown), the overlap of differentially expressed genes from the 105 list with the 492 list was 15.2 %., With respect to shared pathways, we identified that the KEGG Complement and Coagulation Cascades pathway was also enriched in the progressors assigned to radiation treatment (Additional file 3: Table S4 and S5). Other pathways enriched in progressors assigned to radiation included MAPK signaling (KEGG), Cell Adhesion (KEGG), FGFR Ligand Binding and Activation, GPCR ligand binding, PI-3 K Cascades, and Cell-Cell Junction Organization (all Reactome) (Additional file 3: Table S5). Of note, in this cohort as well, interactions with the extracellular environment appeared to be important for progression. While this was an interesting finding, given our modest sample size after stratification by treatment assignment, we carried out the rest of our analyses with the entire cohort. We next moved to look at somatic mutations in the entire cohort, comparing progressors and nonprogressors.
Pathways enriched for genes with increased frequency of mutations in the Progressor cohort (Differential > 5 %) compared to NonProgressors
FDR adjusted P-value
Pre-NOTCH Transcription and Translation (Reactome)
CREBBP, NOTCH2, TP53*
Ion transport by P-type ATPases (Reactome)
ATP10B, ATP2C1, ATP8B4**
ECM Receptor Interaction (KEGG)
RELN, LAMA2, ITGA7
Glycosaminoglycan metabolism (Reactome)
CHSY3, CSGALNACT1, B3GNT7
Pathways enriched for DE genes (Progressor vs NonProgressor) with CTCF Binding Sites
FDR adjusted P-value
Beta3 integrin Cell Surface Interactions (PID)
FGG, FGB, VTN a
P130Cas linkage to MAPK signaling for Integrins (Reactome)
FGG, FGB, VTN a
Response to Elevated Platelet Cytosolic CA2+ (Reactome)
FGG, FGB, ALB
Fibrinolysis Pathway (Biocarta)
Integrin signaling deregulation in HNSCC progressors
As mentioned above, the Reactome Integrin Cell Surface Interaction pathway was significantly enriched for the putative differentially expressed candidate genes (FDR adjusted P-value = 0.00424, Fig. 1). In addition, the frequency of uniquely mutated genes in progressors ranged from 1.4 to 5.8 %, with 26.5 % of the progressors having at least one gene mutated in the pathway (of those, the median was 4.8 % and the range was 4.8–9.5 %). For the progressors, 45.6 % had overexpression of at least one gene in this pathway representation (of those, the median was 9 % of the pathway overexpressed with a range of 4.5 %–27.3 %). When evaluating combined mutation and overexpression, 61.8 % of the progressors had at least one gene aberrant in this pathway. In addition to this specific Reactome pathway, four other Reactome pathways relating to Integrin signaling were deregulated as well as the PID Beta1, Beta2, and Beta3 Integrin Cell Surface Interactions pathways (Additional file 3: Table S4). As previously noted, the radiation treatment assignment cohort also exhibited gene expression aberrancies in extracellular matrix interactions indicating similar biology in this subgroup (Additional file 3: Table S5). This supports the concept that microenvironmental interactions involving integrins are essential for HNSCC progression.
IL-10 signaling alterations in HNSCC progressors
Notably, the IL-10 Anti-inflammatory Signaling pathway was significantly enriched for putative differentially expressed genes among HNSCC progressors (FDR adjusted P-value 0.042, Fig. 2, Additional file 3: Table S4). In addition, HNSCC progressors also harbored unique mutations in several members of the pathway including IL-6, STATs including STAT5A, and BLVRB (Fig. 2). Specifically, 41.2 % of the progressors had at least one gene overexpressed in this pathway (with a maximum of 38.5 % of the pathway overexpressed observed in any progressor) and 11.8 % of the progressors had at least one gene mutated. This was intriguing given the early promise of immune-based therapies in HNSCC. In addition, multiple clotting pathways were found to be differentially expressed (FDR < 0.05) including the Fibrinolysis Pathway (Biocarta), the Intrinsic Prothrombin Activation Pathway (Biocarta), Genes involved in Platelet Aggregation (Reactome), and Complement and Coagulation Cascades (KEGG). These pathways are linked to inflammation  as well, and could also point to the importance of this microenvironmental alteration in HNSCC progression. Again, similar gene expression deregulation was seen in these pathways for the radiation treatment assignment cohort.
Gene signature predicting survival
Cox modeling of molecular and clinical data
Given the prognostic ability of the gene signature above, we were interested in whether we could model survival based on these molecular changes and historically important clinical factors. We examined key clinical features (nodal extracapsular spread, alcohol use, tobacco smoking history, gender, margin status) as well as the combined mutation and expression gene signature using a multivariate Cox proportional hazards modeling approach. While the overall model was significant (P = 0.032), the key factors were interestingly nodal extracapsular spread and alcohol use (P = 0.004 for model with those factors alone).
Mining the TCGA database provides unprecedented opportunities to unravel unique feature of tumorigenesis in HNSCC. Recent published analysis of primary HNSCC molecular alterations in patient data from TCGA reiterated known HNSCC drivers and uncovered distinct molecular alterations between HPV and tobacco-driven tumors . Tumor heterogeneity regardless of the pathways involved was linked to reduced overall survival in a companion study of this population, highlighting the genomically unstable nature of this cancer . The heterogeneity of HNSCC, in part, has limited improved therapeutic targeting of this disease. In fact, cetuximab, targeting EGFR overexpression, is the only targeted agent used in the treatment of HNSCC. However single agent response rates are low, and in combination with standard chemotherapy for progressive disease, overall survival remains less than one year . Many groups are harnessing the power of the TCGA data to characterize molecular changes that might predict survival, exemplified by a recent study suggesting an 11-gene signature was able to predict nodal extracapsular spread and also overall survival . Our study was designed to identify the genomic differences between progressors and nonprogressors at both the DNA and RNA level in order to highlight important pathways associated with progression. Interestingly, we uncovered a significant increase in deleterious mutations of CTCF, which is a master chromatin regulator associated with genomic instability and cancer progression [7, 8]. Deregulation of this gene could be a contributor to the genomic instability and heterogeneity in HNSCC although further mechanistic studies would be required for evaluation. Interestingly, progressors displayed differentially expressed genes harboring CTCF binding sites that participated Integrin-related pathways. This indicates that at least one potential downstream effect of CTCF deregulation could be aberrant microenvironmental interactions involving Integrins facilitating HNSCC progression.
We identifed Integrin and IL-10 signaling as unique prognostic pathways for HNSCC progression. Microenvironmental interaction aberrancies were confirmed by both mutational and expression analysis, and was revealed in progressor cohorts irrespective of their radiation treatment assignment. This implicates tumor microenvironment interactions in the driving biology of tumor progression for all HNSCC patients including those that require radiation as part of their treatment. Intriguingly, both of these pathways have potential promise for guiding targeting therapies. Recently, targeting both EGFR overexpression and Integrin B1 signaling was shown to radiosensitize HNSCC cells, building on previous literature demonstrating Integrin aberrations in HNSCC . Further, cilengitide, an αvβ3 and αvβ5 Integrin inhibitor, has been tested in the recurrent/metastatic HNSCC setting in combination with cytotoxic chemotherapy, however there was no improvement in progression free survival with addition of cilengitide . This should be reevaluated with the improved biomarkers identified in this study or in the definitive rather than metastatic setting. Alterations in IL-10 signaling uncovered in this study suggests to an interesting therapeutic angle. Inflammation is a hallmark of HNSCC progression based on both animal and human studies [11, 12]. IL-10 signaling plays a key role in regulation of cancer-associated inflammation including regulation of CD8 T cells . As targeting of CD8 cells with PD1 (programmed death 1) pathway inhibitors has shown significant promise in multiple similar tumor types, it has emerged as a attractive targetable pathway in HNSCC . Potentially, deregulation of the IL-10 pathway could be used as a biomarker to stratify patients more likely to respond to this therapy. Finally, we identified an 11-gene signature to predict for progression. In addition to alcohol use and nodal extracapsular spread, a known poor prognostic pathologic factor utilized by other groups , this pathway was very powerful in stratifying patients with poor prognosis. Our novel gene signature could be used to identify patients that could benefit from intensified therapy (either concurrent with definitive therapy or adjuvant). The limitations of our study include that we were unable to stratify by HPV status given incomplete clinical annotation within the TCGA dataset, we did not have access to recurrent tumor tissue, and we did not stratify by stage. Nevertheless we were able to uncover significant aberrant pathways that, after further mechanistic validation, have potential to open new avenues for therapeutic treatment of recurrent HNSCC.
Selection of patients and study design
TCGA HNSCC data used for this analysis were time stamped August 13th 2014 and downloaded from the TCGA Research Network: http://cancergenome.nih.gov/. Data types utilized were the clinical data (patient demographics, drug therapy, radiation therapy, and follow-up), RNA-Seq V2 (Level 3; Illumina HiSeq 2000), and somatic mutations (Level 2; Illumina Genome Analyzer DNA Sequencing). All data was mapped to genome build hg 19.
Patients were first classified as progressor or nonprogressor based on follow-up annotation, specificaly the presence or absence of a new tumor event. We required annotation to confirm the tumor event (days to new tumor and/or new tumor anatomical location). All patients were required to have treatment annotation in addition to the follow-up data. All samples used in this study were collected from initial pre-treatment diagnosis (the samples had not been exposed to chemotherapy or radiation).
In-house workflows in the R Statistical Programming environment were used for all QA/QC and analysis . The clinical data was checked for duplicated rows, blank fields and other quality checks. All of the clinical data sets were merged together by the common BCR Patient Barcodes. Differential expression (DE) analysis between progressors and nonprogressors was conducted by fitting linear models using the edgeR framework . As a secondary analysis, we also examined differential expression among progressors and nonprogressors with radiation treatment assignment. For all DE, P-values were False Discovery Rate Adjusted. Genes with low expression in all samples (<1 (count per million) were flagged and filtered out. Somatic mutation data was also filtered out if there were any tuples with no known gene symbols in RNA-Seq V2. All gene symbols were verified to have approved symbols or synonyms. Cytoscape was used for stylized pathway visualization  specifically the Reactome FI Cytoscape Plugin 4 .
Somatic mutations for progressors and nonprogressors were evaluated by first assessing gene symbol, chromosome, and start and stop. The distribution of truncating mutations (nonsense, nonstop, frameshift deletion, frameshift insertion, and splice site) as well as missense mutation rates was compared between progressors and nonprogressors. For the entire TCGA HNSCC cohort, MutSig2CV annotation (ranking and significance) was examined to assess unique mutations. MutSig analyzes mutations to identify genes that were mutated more often than expected by chance given background mutation processes . Fisher’s exact tests were performed to examine differences in mutational frequency by mutational class in those genes with overall mutational frequency difference of 5 % or more between progressors and nonprogressors. Lollipop figures of mutational type by gene were generated by the MutationMapper visualization tool (courtesy of Memorial Sloan-Kettering Cancer Center). Somatic mutations were counted both as the number of total of variants, as well as summarized at the gene level as the total number of genes mutated. For mutations unique to progressors and unique to nonprogressors (noting the caveat this can in some cases be due to sampling), we computed the ratio of variants to gene level mutation (in both cases only for those mutations unique to each group) = # of somatic mutations (individual variants) / # of genes with somatic mutations. This allowed us to assess the number of mutations relative to the number of genes mutated to understand the impact on the genome. A large ratio could be indicative a concentration of highly mutated genes.
Both differentially expressed genes as well as somatic mutation data were annotated to pathway models from Reactome, KEGG, Pathway Interaction Database (PID), and Biocarta from MSIGDB . We then examined if there was significant enrichment of these candidate genes in the pathways. As with the differential expression analysis, all enrichment P-values were False Discovery Rate Adjusted.
Mutation and gene expression data was overlaid to identify an aggregate gene signature (based on both differentially expression and unique mutation in progressors only). Overall Survival Kaplan-Maier estimates based on alternations in this signature were examined. Both clinical features (nodal extracapsular spread, alcohol use, tobacco smoking history, gender, margin status) as well as the aggregated gene signature were examined by a multivariate Cox proportional hazards modeling approach.
The results published here are in whole or part based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/. Support for this work was provided by NIH/NCI (5P30CA06533), the Oregon Health & Science University Knight Cancer Institute Cathy and Jim Rudd Career Development Award for Cancer Research, and the Oregon Health & Science University Medical Research Foundation Early Clinical Investigator Grant. We also wish to thank the reviewers and editors for their insightful comments and suggestions.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Bonner JA, Harari PM, Giralt J, Cohen RB, Jones CU, Sur RK, et al. Radiotherapy plus cetuximab for locoregionally advanced head and neck cancer: 5-year survival data from a phase 3 randomised trial, and relation between cetuximab-induced rash and survival. Lancet Oncol 2010, 11(1):21-28.Google Scholar
- Vermorken JB, Mesia R, Rivera F, Remenar E, Kawecki A, Rottey S, et al. Platinum-based chemotherapy plus cetuximab in head and neck cancer. N Engl J Med 2008, 359(11):1116-1127.Google Scholar
- Cancer Genome Atlas N: Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature 2015, 517(7536):576-582.Google Scholar
- Oikonomopoulou K, Ricklin D, Ward PA, Lambris JD. Interactions between coagulation and complement--their role in inflammation. Seminars in immunopathology 2012, 34(1):151-165.Google Scholar
- Mroz EA, Tward AM, Hammon RJ, Ren Y, Rocco JW. Intra-tumor Genetic Heterogeneity and Mortality in Head and Neck Cancer: Analysis of Data from The Cancer Genome Atlas. PLoS medicine 2015, 12(2):e1001786.Google Scholar
- Wang W, Lim WK, Leong HS, Chong FT, Lim TK, Tan DS, et al. An eleven gene molecular signature for extra-capsular spread in oral squamous cell carcinoma serves as a prognosticator of outcome in patients without nodal metastases. Oral Oncol 2015.Google Scholar
- Libby RT, Hagerman KA, Pineda VV, Lau R, Cho DH, Baccam SL, et al. CTCF cis-regulates trinucleotide repeat instability in an epigenetic manner: a novel basis for mutational hot spot determination. PLoS genetics 2008, 4(11):e1000257.Google Scholar
- Marshall AD, Bailey CG, Rasko JE. CTCF and BORIS in genome regulation and cancer. Current opinion in genetics & development 2014, 24:8-15.Google Scholar
- Eke I, Zscheppang K, Dickreuter E, Hickmann L, Mazzeo E, Unger K, et al. Simultaneous beta1 integrin-EGFR Targeting and Radiosensitization of Human Head and Neck Cancer. J Natl Cancer Inst 2015, 107(2).Google Scholar
- Vermorken JB, Peyrade F, Krauss J, Mesia R, Remenar E, Gauler TC, et al. Cisplatin, 5-fluorouracil, and cetuximab (PFE) with or without cilengitide in recurrent/metastatic squamous cell carcinoma of the head and neck: results of the randomized phase I/II ADVANTAGE trial (phase II part). Ann Oncol 2014, 25(3):682-688.Google Scholar
- Bonomi M, Patsias A, Posner M, Sikora A. The role of inflammation in head and neck cancer. Adv Exp Med Biol 2014, 816:107-127.Google Scholar
- Bornstein S, White R, Malkoski S, Oka M, Han G, Cleaver T, et al. Smad4 loss in mice causes spontaneous head and neck cancer with increased genomic instability and inflammation. The Journal of clinical investigation 2009, 119(11):3408-3419.Google Scholar
- Oft M: IL-10: master switch from tumor-promoting inflammation to antitumor immunity. Cancer immunology research 2014, 2(3):194-199.Google Scholar
- Zandberg DP, Strome SE. The role of the PD-L1:PD-1 pathway in squamous cell carcinoma of the head and neck. Oral Oncol 2014, 50(7):627-632.Google Scholar
- R Core Team (2014). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.
- Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.PubMedPubMed CentralView ArticleGoogle Scholar
- Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome research 2003, 13(11):2498-2504.Google Scholar
- Wu G, Feng X, Stein L. A human functional protein interaction network and its application to cancer data analysis. Genome Biol. 2010;11(5):R53.PubMedPubMed CentralView ArticleGoogle Scholar
- Lawrence MS, Stojanov P, Polak P, Kryukov GV, Cibulskis K, Sivachenko A, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499(7457):214–8.PubMedPubMed CentralView ArticleGoogle Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.PubMedPubMed CentralView ArticleGoogle Scholar