Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Identification of immune‐related genes contributing to head and neck squamous cell carcinoma development using weighted gene co‐expression network analysis

Identification of immune‐related genes contributing to head and neck squamous cell carcinoma... INTRODUCTIONHead and neck squamous cell carcinoma (HNSCC) is a common malignancy, with an annual incidence of 500 000 cases worldwide.1,2 It is well‐established that most patients with HNSCC are diagnosed with advanced or metastatic disease at diagnosis, with a 5‐year overall survival (OS) rate of only 50% and a recurrence or metastasis rate of 30%.3,4 Indeed, once recurrence and metastasis occur, the prognosis is poor, and the median survival time is only 6 months. These findings emphasize the need to identify novel approaches to reduce recurrence and metastasis rates and improve the survival of patients with recurrent and metastatic HNSCC.In recent years, tumor immunity research has gathered momentum given the efficacy of immune checkpoint inhibitor (ICI) treatment in various tumors.5 An increasing body of evidence suggests that immune cell infiltration is essential for the efficacy of ICIs6–8 and the cause and predictor of cancer recurrence and metastasis.6,7 It is well‐established that ICIs are efficient in HNSCC8; however, predicting the efficacy of ICIs and the mechanism underlying ICI resistance remain unclear. Current evidence suggests that the tumor immunosuppressive microenvironment can inhibit immune cell function and infiltration levels.9 Accordingly, it is necessary to explore the mechanisms underlying the inhibition of immune cell infiltration in HNSCC.The rapid advent of bioinformatics has enhanced our understanding of the biological functions of tumors, providing the foothold for the discovery and verification of tumor biomarkers. There is ample evidence showing that immune‐related genes in HNSCC can predict tumor prognosis and are potential biomarkers.10,11 However, these genes are known to be inversely related to immunity. Importantly, bioinformatic analysis provides a novel approach to uncover genes related to immune cell infiltration and explore their biological functions and their potential clinical translational value as diagnostic and prognostic markers.Herein, weighted correlation network analysis (WGCNA) was conducted to identify highly correlated gene clusters using module clustering.12 We sought to identify specific immune‐related genes associated with the development or prognosis of HNSCC from The Cancer Genome Atlas (TCGA). Furthermore, we explored key genes related to the occurrence of HNSCC, the biological process involved, and their prognostic value in this patient population.MATERIALS AND METHODSData collection and preprocessingThe TCGA (https://tcgadata.nci.nih.gov/tcga/) HNSCC transcriptomic and clinical data, and two Gene Expression Omnibus (GEO) datasets (GSE65858 and GSE41613) were downloaded. A total of 528 HNSCC patients were identified in the TCGA database, consisting predominantly of males (73.11%) with no previous history of cancer (94.89%). We used DESeq213 to analyze the differential expression of genes in the transcriptomic data of TCGA HNSCC database. The screening criteria for differentially expressed genes (DEGs) included: |log2 (fold change)|>2 and p < .05.In addition, the expression levels of key genes of two HNSCC datasets from the Oncomine database (Ye Head‐Neck,14 and Peng Head‐Neck),15 were included in our analyses.Weighted correlation network analysisWGCNA16 was used to identify co‐expressed gene modules in DEGs and convert the adjacency matrix into a topological overlap matrix (TOM). According to the TOM‐based dissimilarity, DEGs were finally divided into 13 different modules. Using the InnateDB17 database, a database of genes involved in innate immunity, we identified 824 immune‐related genes, of which 72 exhibited significant changes in HNSCC. Furthermore, genes from the brown module were selected as the gene set for subsequent analysis since it contained the most significant number of immune‐related DEGs (n = 223).Gene ontology function enrichment analysisMetascape18 (https://metascape.org/), an online analysis tool, was used to perform Gene ontology (GO) function enrichment analysis (Biological Process, Molecular Function, Cellular Component) on DEGs, with the parameters set to min overlap = 3, p‐value cutoff = .01, and min enrichment = 1.5.19 Specifically, the immunological signature (Immunologic Signature) was selected for functional enrichment analysis, and the parameters remain unchanged.Gene set enrichment analysisThe Gene set enrichment analysis (GSEA)20 method was used to screen the biological state or process related to the occurrence of HNSC. We selected “h.all.v7.1.symbols.gmt” as the reference set, the number of permutations was set to 1000, and the permutation type set to the selected gene set. Gene sets with an FDR of <0.05 were significantly enriched, and the Normalized Enrichment Score (NES) was used to indicate the degree of enrichment.Protein–protein interaction network constructionWe obtained the protein–protein interactions (PPI) (score >0.4) of 223 immune‐related DEGs using STRING database21 (https://string-db.org/), and the main interaction network was visualized using Cytoscape.22 The top 10 key genes were obtained using CytoHubba, and the interaction network was drawn.Tumor immune infiltration analysisWe analyzed the gene expression data of HNSCC samples from TCGA using TIMER2.023 and analyzed the correlation between KRT4, KRT78, KRT13, and SPRR3 expression and immune cell (CD8+ T cells and macrophages) infiltration levels.ImmunohistochemistryThis study was approved by the Review Board Committee of Fujian Cancer Hospital, Fujian, China. The Review Board Committee granted a waiver of informed consent for this study. The informed consent Tissue samples from cases of pathologically diagnosed HNSCC (n = 70) treated at the Fujian Cancer Hospital from 2008 to 2017, and normal tissue samples (n = 10) were included in the analysis. The clinicopathological features of HNSCC and normal samples are shown in Table 3. The samples were fixed in formaldehyde and processed with heat‐mediated antigen retrieval in citrate buffer (pH = 6). The samples were then blocked and incubated with the following primary antibodies: rabbit polyclonal anti‐KRT13 (1:1000, Cat No. A116809, SIGMA Life Science, USA), rabbit polyclonal anti‐KRT78 (1:1000, 000015851, SIGMA Life Science, USA), and rabbit polyclonal anti‐SPRR3 (1:200, Cat No. ab218131, Abcam, USA) at 4°C overnight. The ElivisionTM plus Polyer HP (Mouse/Rabbit) immunohistochemistry (IHC) Kit (Cat. KIT‐9901, MXB biotechnologies, China) was used for IHC detection. Two independent pathologists, blinded to the clinicopathological data, evaluated the immunohistochemical score.Statistical analysisThe Wilcoxon rank‐sum test and Wilcoxon signed‐rank test were used to analyze the expression of key genes in tumor and normal tissues. The relationship between clinicopathological features and KRT4, KRT78, KRT13, and SPRR3 expression was assessed by the Kruskal–Wallis test, Wilcoxon signed‐rank test, and chi‐squared test. Survival curves were drawn using the Kaplan–Meier method, and differences between groups were assessed using the log‐rank test. OS was defined as the time of diagnosis until the date of death from any cause or last follow‐up. Progression‐free survival (PFS) was defined as the time of diagnosis until the date of disease progression, death or last follow‐up. A p‐value <.05 (two‐sided) was statistically significant. Statistical analyses were carried out using R (version 3.6.1) and SPSS (version 24.0).RESULTSIdentification and functional enrichment analysis of DEGs related to the occurrence of HNSCCTo explore the relationship between HNSCC development and immune infiltration, we conducted a comprehensive integrated analysis of transcriptomic and clinical data of HNSCC patients from TCGA (Figure 1, Table 1). First, differential expression analysis showed that 1869 and 1578 genes were significantly upregulated and downregulated in HNSCC, respectively (Figure 2A, Table S1).1FIGUREFlow chart of overall analysis1TABLEClinical characteristics of patients in the The Cancer Genome Atlas databaseClinical characteristicsNumber (%)GenderMale386 (73.11)Female142 (26.89)Prior malignancyYes27 (5.11)No501 (94.89)Vital statusAlive304 (57.58)Dead224 (42.42)AJCC clinical TT137 (7.01)T2152 (28.79)T3139 (26.33)T4184 (34.85)Not reported16 (3.03)AJCC clinical NN0246 (46.59)N185 (16.10)N2166 (31.44)N39 (1.70)Not reported22 (4.17)AJCC clinical MM0496 (93.94)M16 (1.14)Not reported26 (4.92)AJCC clinical stageI21 (3.98)II99 (18.75)III107 (20.27)IV287 (54.36)Not reported14 (2.65)2FIGUREIdentification and functional enrichment analysis of differentially expressed genes in head and neck squamous cell carcinoma (HNSCC)GO annotation demonstrated that these DEGs were mainly enriched in biological processes, including muscle contraction, myofibril assembly, and muscle structural development (Figure 2B). GO terms for molecular functions consisted of structural molecular activity, muscle structural components, and receptor ligand activity (Figure 2C), and the related cellular components included contractile fibers, extracellular matrix, and A‐bands (Figure 2D). Meanwhile, enrichment analysis of the immune characteristics of these DEGs showed that they were involved in the regulation of CD4+ T cells, macrophages, and CD8+ T cells (Figure 2E, Table S2). Furthermore, GSEA showed that biological processes that were significantly enriched in normal tissues included oxidative phosphorylation (Figure 3A), adipogenesis (Figure 3B), myogenesis (Figure 3C), bile acid metabolism (Figure 3D), fatty acid metabolism (Figure 3E), early estrogen response (Figure 3F), peroxisomes, (Figure 3G) and xenobiotic metabolism (Figure 3H). The biological processes significantly enriched in HNSCC encompassed E2F transcription factor targeting (Figure 3I), G2M checkpoint (Figure 3J), interferon α response (Figure 3K), epithelial‐mesenchymal transition (Figure 3L), interferon γ response (Figure 3M), MYC targeting (Figure 3N), DNA repair (Figure 3O) and mitotic spindle (Figure 3P; Table S3).3FIGUREGene set enrichment analysis (GSEA) of differentially expressed genes (DEGs) related to head and neck squamous cell carcinomaCo‐expression analysis of DEGs and PPI network in HNSCCTo identify gene sets related to the occurrence of HNSCC, we performed WGCNA based on transcriptomic data. Using a soft‐threshold of 4 (Figure 4A,B), WGCNA identified 13 modules (Figure 4C–E). Compared with the other modules, the brown module contained the most significant number of immune‐related genes (n = 7; Table S4). Two hundred and twenty‐three genes were enriched in the brown module, and many interactions were identified between the proteins encoded by these genes (Figure 5A,B). Proteins encoded by PPL, SCEL, KRT4, KRT24, KRT78, KRT13, SPRR3, TGM3, CRCT1, and CRNN were core components in the PPI network (Figure 5C).4FIGURECo‐expression analysis of differentially expressed genes (DEGs) in head and neck squamous cell carcinoma5FIGUREProtein–protein interaction (PPI) network of immune‐related differentially expressed genes (DEGs).Differential expression, immune cell infiltration, and survival analysis of key genesThe expression of key DEGs screened from normal and HNSCC tissues was analyzed. We found that their expression levels in the HNSCC group were significantly lower than in the normal group (Figure 6A–J, left). In the meantime, we found that except for KRT24 and CRCT1 (AUCs <0.75), the other eight genes yielded a high predictive performance for the occurrence of HNSCC, with AUC values of 0.828 (PPL), 0.804 (SECL), 0.854 (KRT4), 0.799 (KRT78), 0.824 (KRT13), 0.797 (SPRR3), 0.822 (TGM3), and 0.852 (CRNN; Figure 6A–J, right). Further analyses of immune cell infiltration of these eight genes showed that the expression levels of KRT4, KRT78, KRT13, and SPRR3 were significantly correlated with the infiltration levels of CD8+ T cells and macrophages in HNSCC (Figure 7). Then, survival analyses of these four genes were conducted, including OS and DFS, and results showed that except for KRT4, the other three genes were related to OS and DFS in HNSCC (Figure 8A). KRT78 expression was significantly related to OS and DFI (Figure 8B). KRT13 and SPRR3 expression levels were significantly related to DFI (Figure 8C) and OS (Figure 8D), respectively. Finally, KRT4, KRT78, and SPRR3 were selected as our candidate genes.6FIGUREThe RNA expression of key genes in head and neck squamous cell carcinoma7FIGUREAccurate prediction of the correlation between key genes and immune cell infiltration in head and neck squamous cell carcinoma8FIGURESurvival analysis of key genes in head and neck squamous cell carcinomaCorrelation analysis of key gene expression and clinical characteristicsA low correlation was found between KRT78 and SPRR3 expression levels and clinical T stage. However, the expression levels of these two genes in T4 HNSCC patients were significantly lower than in T3 HNSCC patients (Figure 9A,B), while no correlation was found between KRT13 expression and T stage (Figure 9C). The expression levels of these three candidate genes showed a decreasing trend as the N stage increased (Figure 9D–F). Similarly, KRT78, KRT13, and SPRR3 expression levels showed a downward trend with increasing HNSCC clinical stage (Figure 9G–I). These results suggested that the three candidate genes were closely related to HNSCC development (Table 2).9FIGURECorrelation between the expression level of candidate genes and the clinical stage of head and neck squamous cell carcinoma2TABLEDemographic and clinical characteristics of head and neck squamous cell carcinoma (HNSCC) patients in GSE41613 and GSE65858CharactersGSE41613GSE65858n (%)n (%)GenderFemale31 (32.0)43 (17.4)Male66 (68.0)223 (82.6)Age (year)<6050 (51.5)153 (56.7)≥6047 (48.5)117 (43.3)HPVNegative97 (100)197 (73.0)Positive0 (0)73 (27.0T classificationT1NA115(42.6)T3NA155 (57.4)N classificationN0‐N1NA126 (46.7)N2‐N3NA144 (53.3)Clinical stageI‐II41 (42.3)55 (20.4)III‐IV56 (57.7)215 (79.6)Comprehensive treatmentNo44 (45.4)Yes53 (54.6)External validation of candidate genesWe used Oncomine datasets (Ye Head‐Neck, and Peng Head‐Neck) to verify the expression levels of KRT13, KRT78, and SPRR3 in HNSCC tumor and normal tissues. Analysis of the Ye Head‐Neck (Figure 10A–C) and Peng Head‐Neck (Figure 10D–F) datasets showed that the expression levels of KRT13 (p < .001), KRT78 (p < .001), and SPRR3 (p < .001) in HNSCC were significantly lower than in normal tissues.10FIGUREVerification of candidate gene mRNA expression levels in Oncomine datasetsTo further analyze differences in the expression of KRT13, KRT78, and SPRR3 genes in HNSCC and normal tissues, we used a microarray of 10 normal tissues and 70 HNSCC tissues to verify the protein expression of the above genes. The results showed that the protein expression levels of KRT13 (p = .042, Figure 11A), KRT78 (p < .001, Figure 11B), and SPRR3 (p = .022, Figure 11C) in HNSCC tissues were significantly lower than in normal tissues (Table 3).11FIGUREValidation the protein expression of KRT13, KRT78, and SPRR3 in head and neck squamous cell carcinoma tissues and normal tissues by immunohistochemistry3TABLEKRT13, KRT78, and SPRR3 protein expression of head and neck squamous cell carcinoma (HNSCC) patients (n = 70) and normal tissue (n = 10) in validation cohort of protein expressionCharactersLevelHNSCCNormal tissuepGenderFemale5 (7.1%)2 (20.0%)0.210Male65 (92.9%)7 (80.0%)Age<60 years32 (45.7%)6 (60.0%)0.320≥60 years38 (44.4%)4 (40.0%)T classificationT1‐T230 (42.9%)T3‐T440 (57.1%)N classificationN0‐155 (78.6%)N2‐315 (21.4%)Clinical stageStage I–II26 (37.1%)Stage III–IV44 (62.9%)Histologic gradeG129 (59.3%)G221 (22.2%)G320 (18.5%)Anatomic siteOral cavity15 (29.6%)3 (30.0%)0.385Oropharynx4 (3.7%)2 (20.0%)Hypopharynx4 (3.7%)1 (10.0%)Larynx47 (63.0%)5 (50.0%)KRT13Negative29 (41.4)00.042Week12 (17.1)2 (20)Moderate16 (22.9)3 (30)Strong13 (18.6)5 (50)KRT78Negative40 (57.1)0 (0)<0.001Week17 (24.3)3 (30.0)Moderate8 (12.3)1 (10.0)Strong4 (5.7)6 (60.0)SPRR3Negative18 (25.7)00.022Week33(47.1)3 (30.0)Moderate12 (17.1)3 (30.0)Strong7 (10.0)4 (40.0)After verifying the expression differences of KRT13, KRT78, and SPRR3 in HNSCC and normal tissues, we used GEO datasets (GSE65858 and GSE41613) to analyze their prognostic value (Table 2). In the GSE65858 dataset, survival analysis demonstrated that patients with low KRT13 (p = .044), KRT78 (p = .0086) and SPRR3 (p = .017) expression had worse OS than those with high expression (Figure 12A–C). Furthermore, analysis of the GSE41613 dataset showed that patients with low KRT78 (p = .005) and SPRR3 (p = .02) expression had worse OS than those with low expression, while no significant difference in OS was found in patients with low and high expression of KRT13 (p = .96; Figure 11D–F).12FIGUREValidation the prognostic value of KRT13, KRT78, and SPRR3 mRNA expression in head and neck squamous cell carcinoma from GSE65858 and GSE41613DISCUSSIONThe present study explored the relationship between HNSCC development and immune infiltration via a comprehensive integrated analysis of transcriptomic and clinical data of HNSCC patients in public databases. Ten key genes, namely PPL, SCEL, KRT4, KRT24, KRT78, KRT13, SPRR3, TGM3, CRCT1, and CRNN, were screened by differential expression and PPI network analyses. We found that KRT4, KRT78, and SPRR3 were downregulated in HNSCC and negatively correlated with immune cell infiltration. Furthermore, survival analyses demonstrated that these genes were negatively correlated with OS and PFS in HNSCC. Analysis of the Oncomine and GEO datasets validated that the KRT13, KRT78, and SPRR3 expression levels were downregulated in HNSCC and negatively correlated with the prognosis of patients with HNSCC. Finally, IHC showed that KRT13, KRT78, and SPRR3 protein expression levels were downregulated in HNSCC compared to normal tissues.Moreover, we demonstrated that 1869 genes and 1578 genes were significantly upregulated and downregulated in HNSCC. Enrichment analysis showed that these DEGs regulate CD4+ T cells, macrophages, and CD8+ T cells. Furthermore, GSEA analysis showed that interferon α and γ reactions were also enriched in HNSCC, suggesting that the immune microenvironment of HNSCC was significantly different from that of normal tissues. It is widely acknowledged that the microenvironment is an essential factor leading to the occurrence and development of tumors,23 affecting the efficacy of immunotherapy.24,25 Importantly, WGCNA can harness the information of thousands or tens of thousands of the greatest gene expression changes or all genes to identify the gene set of interest and perform significant association analysis with the phenotype.26 In the present study, WGCNA showed that out of the 13 modules, the brown module contained the most immune‐related genes (n = 7).PPI network analysis showed that PPL, SCEL, KRT4, KRT24, KRT78, KRT13, SPRR3, TGM3, CRCT1, and CRNN were core components downregulated in HNSCC. KRT4, KRT13, KRT78, and SPRR3 were negatively correlated with CD8+ T cells and macrophage infiltration, suggesting that they mediated the formation of an inhibitory immune micro‐environment.Interestingly, these four key genes have been documented in keratinization, suggesting that keratinization not only reflects the cell differentiation level but also determines the immune state.Besides, we also revealed that the expression levels of KRT4, KRT78, KRT13, and SPRR3 in TCGA‐HNSCC negatively correlated with immune cell infiltration and influenced patient prognosis. Analysis of the two Oncomine datasets also confirmed that KRT13, KRT78, and SPRP3 expression levels were significantly lower in HNSCC than in normal tissues. In addition to the significant reduction in mRNA expression levels of KRT13, KRT78, and SPRR3 in HNSCC, we also documented a significant decline in the protein expression levels of KRT13, KRT78, and SPRR3 in HNSCC using IHC. This means that these 3 genes are related to the immune infiltration and prognosis of NPC, and they are potential immune‐related genes. Many previous studies had found that multiple immune‐related genes were associated with the prognosis of NPC. These genes were known immune‐related genes, such as CCR5, CD3E, CD4, SFRP4, SFRP4, CPXM1, and COL5A1 CCR6, CCL22, ROBO1, DKK1 and PDGFA, and so forth.11,27,28 In our study, WGCNA was used to lock immune‐related modules and then protein–protein interaction network was used to screen out potential immune‐related genes. The subsequent analysis also found that KRT13, KRT78, and SPRP3 were negatively correlated with immune cell infiltration. KRT13, KRT78, and SPRP3 are newly discovered immune‐related genes in HNSCC in our study.Herein, we found that patients with low expression of KRT13, KRT78, and SPRR3 in the TCGA‐HNSCC dataset had worse OS and PFI than those with high expression. This finding was also externally validated in two public GEO datasets (GSE65858 and GSE41613), where the low expression of KRT78 and SPRR3 in HNSCC patients was associated with poor OS. Our results indicated that the main biological function of Keratin 78 and SPRP3 is to participate in developmental biology and cell keratinization. Hence, these two key genes can inhibit the occurrence and progression of HNSCC acting as tumor suppressor genes. Although preliminary studies have shown that KRT7829 and SPRR330 are involved in cytokeratinization, their role in tumors remains unknown. Current evidence suggests that SPRR3 expression is low in HNSCC26 and esophageal cancer,31 while other studies have substantiated that SPRR3 is a prognostic factor of esophageal cancer and a sensitive marker of chemotherapy and radiotherapy.32–34 At present, few reports have shown that KRT78 is differentially expressed in melanoma,35 colon adenocarcinoma,36 cervical cancer,37 urinary bladder cancer,38 and HNSCC,39 and may be related to tumor prognosis.37,38 However, the biological function of KRT78 in tumors remains primarily understudied, warranting further in‐depth studies.Although this study enhanced our understanding of the involvement of KRT78 and SPRR3 in HNSCC, there were some limitations and shortcomings. First, the detailed mechanisms underlying the KRT78‐ and SPRR3‐mediated decrease in immune cell infiltration and biological functions in HNSCC were not explored. Indeed, in vivo and in vitro experiments are needed to explore and verify their biological functions. Although our study confirmed that KRT13, KRT78, and SPRR3 protein expression levels in HNSCC were significantly downregulated, the prognostic value of protein expression could not be analyzed due to the lack of survival data. It is necessary to conduct a perspective study to demonstrate the prognostic value of KRT13, KRT78, and SPRR3 protein expression.To conclude, we identified four hitherto unrecognized key genes, KRT4, KRT78, KRT13, and SPRR3, related to the occurrence and development of HNSCC and positively correlated with immune cell infiltration. KRT78 and SPRR3, novel immune‐related genes, might serve as diagnostic and prognostic biomarkers of HNSCC. Nonetheless, the detailed biological functions and clinical value of KRT78 and SPRR3 in HNSCC need further exploration.AUTHOR CONTRIBUTIONSQiaojuan Guo: Conceptualization (equal); investigation (equal); writing – original draft (equal); writing – review and editing (equal). Tianzhu Lu: Formal analysis (equal); writing – original draft (equal). Hanchuan Xu: Formal analysis (supporting); investigation (supporting); writing – original draft (supporting). Qingfeng Luo: Data curation (supporting); formal analysis (supporting). Zhilang Liu: Data curation (supporting); formal analysis (supporting). Sicong Jiang: Date curation(supporting); formal analysis (supporting). Shao jun Lin: Writing – original draft (supporting). Jianji Pan: Writing – original draft (supporting); writing – review and editing (supporting). Mengyao Lin: Formal analysis (supporting); writing – original draft (supporting); writing – review and editing (supporting). Fang Guo: Conceptualization (lead); writing – original draft (lead); writing – review and editing (lead).ACKNOWLEDGMENTSThis work was supported by Immuno Radiotherapy research fund project of Radiation Oncology Branch of Chinese Medical Association (No. Z‐2017‐24‐2020); National Natural Science Foundation of China (Nos. 81860664 and 81 660 453); Fujian Provincial Natural Science Foundation of China (No. 2019J01194); NHC Key Laboratory of Personalized Diagnosis and Treatment of Nasopharyngeal Carcinoma. The authors thank editage (http://www.fabiao@editage.cn) for editing this manuscript.CONFLICT OF INTEREST STATEMENTThe authors declare no conflict of interest.DATA AVAILABILITY STATEMENTThe datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.ETHICS STATEMENTThis study was approved by the Hospital Review Board of Fujian Cancer Hospital, Fujian, China. The Review Board Committee granted a waiver of informed consent for this study.REFERENCESSiegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. 2018;68(1):7‐30.Belcher R, Hayes K, Fedewa S, Chen AY. Current treatment of head and neck squamous cell cancer. J Surg Oncol. 2014;110(5):551‐574.Thariat J, Vignot S, Lapierre A, et al. Integrating genomics in head and neck cancer treatment: promises and pitfalls. Crit Rev Oncol Hematol. 2015;95(3):397‐406.Mehanna H, Wong WL, Dunn J. Management of Advanced Head and Neck Cancer. N Engl J Med. 2016;375(5):492‐493.de Miguel M, Calvo E. Clinical challenges of immune checkpoint inhibitors. Cancer Cell. 2020;38(3):326‐333.Pagès F, Mlecnik B, Marliot F, et al. International validation of the consensus Immunoscore for the classification of colon cancer: a prognostic and accuracy study. Lancet (London, England). 2018;391(10135):2128‐2139.Xu X, Tan Y, Qian Y, et al. Clinicopathologic and prognostic significance of tumor‐infiltrating CD8+ T cells in patients with hepatocellular carcinoma: a meta‐analysis. Medicine. 2019;98(2):e13923.Borel C, Jung AC, Burgy M. Immunotherapy breakthroughs in the treatment of recurrent or metastatic head and neck squamous cell carcinoma. Cancer. 2020;12(9):2691.Lim AR, Rathmell WK, Rathmell JC. The tumor microenvironment as a metabolic barrier to effector T cells and immunotherapy. Elife. 2020;9:e55185.Zhou D, Wang J, Wang J, Liu X. Profiles of immune cell infiltration and immune‐related genes in the tumor microenvironment of HNSCC with or without HPV infection. Am J Transl Res. 2021;13(4):2163‐2180.Chen Y, Li ZY, Zhou GQ, Sun Y. An immune‐related gene prognostic index for head and neck squamous cell carcinoma. Clin Cancer Res. 2021;27(1):330‐341.Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA‐seq data with DESeq2. Genome Biol. 2014;15(12):550.Ye H, Yu T, Temam S, et al. Transcriptomic dissection of tongue squamous cell carcinoma. BMC Genomics. 2008;9:69.Peng S, Creighton CJ, Zhang Y, et al. Tumor grafts derived from patients with head and neck squamous carcinoma authentically maintain the molecular and histologic characteristics of human cancers. J Transl Med. 2013;11:198.Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.Breuer K, Foroushani AK, Laird MR, et al. InnateDB: systems biology of innate immunity and beyond‐recent updates and continuing curation. Nucleic Acids Res. 2013;41:D1228‐33.Zhou Y, Zhou B, Pache L, et al. Metascape provides a biologist‐oriented resource for the analysis of systems‐level datasets. Nat Commun. 2019;10(1):1523.Tuo Z, Zhang Y, Wang X, et al. RUNX1 is a promising prognostic biomarker and related to immune infiltrates of cancer‐associated fibroblasts in human cancers. BMC Cancer. 2022;22(1):523.Subramanian A, Kuehn H, Gould J, Tamayo P, Mesirov JP. GSEA‐P: a desktop application for gene set enrichment analysis. Bioinformatics (Oxford, England). 2007;23(23):3251‐3253.von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258‐261.Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498‐2504.Li T, Fan J, Wang B, et al. TIMER: a web server for comprehensive analysis of tumor‐infiltrating immune cells. Cancer Res. 2017;77(21):e108‐e110.Giannone G, Ghisoni E, Genta S, et al. Immuno‐metabolism and microenvironment in cancer: key players for immunotherapy. Int J Mol Sci. 2020;21(12):4414.Bader JE, Voss K, Rathmell JC. Targeting metabolism to improve the tumor microenvironment for cancer immunotherapy. Mol Cell. 2020;78(6):1019‐1033.Yu L, Yang Z, Liu Y, et al. Identification of SPRR3 as a novel diagnostic/prognostic biomarker for oral squamous cell carcinoma via RNA sequencing and bioinformatic analyses. PeerJ. 2020;8:e9393.She Y, Kong X, Ge Y, et al. Immune‐related gene signature for predicting the prognosis of head and neck squamous cell carcinoma. Cancer Cell Int. 2020;20:22.Chengcheng L, Wenwen Q, Ningyue G, et al. Identification of the immune‐related genes in tumor microenvironment that associated with the recurrence of head and neck squamous cell carcinoma. Front Cell Dev Biol. 2021;9:723721.Rogers MA, Edler L, Winter H, Langbein L, Beckmann I, Schweizer J. Characterization of new members of the human type II keratin gene family and a general evaluation of the keratin gene domain on chromosome 12q13.13. J Invest Dermatol. 2005;124(3):536‐544.Koizumi H, Kartasova T, Tanaka H, Ohkawara A, Kuroki T. Differentiation‐associated localization of small proline‐rich protein in normal and diseased human skin. Br J Dermatol. 1996;134(4):686‐692.Souza‐Santos PT, de Oliveira DS, Bernardo V, Lima SC, Rapozo DC, et al. Quantitative evaluation of SPRR3 expression in esophageal squamous cell carcinoma by qPCR and its potential use as a biomarker. Exp Mol Pathol. 2011;91(2):584‐589.Luo A, Chen H, Ding F, et al. Small proline‐rich repeat protein 3 enhances the sensitivity of esophageal cancer cells in response to DNA damage‐induced apoptosis. Mol Oncol. 2013;7(5):955‐967.Luthra R, Wu TT, Luthra MG, et al. Gene expression profiling of localized esophageal carcinomas: association with pathologic response to preoperative chemoradiation. J Clin Oncol. 2006;24(2):259‐267.Zhang Y, Feng YB, Shen XM, et al. Exogenous expression of Esophagin/SPRR3 attenuates the tumorigenicity of esophageal squamous cell carcinoma cells via promoting apoptosis. Int J Cancer. 2008;122(2):260‐266.Wang LX, Li Y, Chen GZ. Network‐based co‐expression analysis for exploring the potential diagnostic biomarkers of metastatic melanoma. PloS One. 2018;13(1):e0190447.Wang H, Liu J, Li J, et al. Identification of gene modules and hub genes in colon adenocarcinoma associated with pathological stage based on WGCNA analysis. Cancer Genet. 2020;242:1‐7.Varga N, Mózes J, Keegan H, et al. The value of a novel panel of cervical cancer biomarkers for triage of HPV positive patients and for detecting disease progression. Pathol Oncol Res. 2017;23(2):295‐305.Chen S, Zhang N, Shao J, Wang T, Wang X. A novel gene signature combination improves the prediction of overall survival in urinary bladder cancer. J Cancer. 2019;10(23):5744‐5753.Zhang X, Feng H, Li Z, et al. Application of weighted gene co‐expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis. Onco Targets Ther. 2018;11:6001‐6021. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Cancer Reports Wiley

Identification of immune‐related genes contributing to head and neck squamous cell carcinoma development using weighted gene co‐expression network analysis

Download PDF
17 pages

Loading next page...
 
/lp/wiley/identification-of-immune-related-genes-contributing-to-head-and-neck-FZdPZ0R3Ey

References (79)

Publisher
Wiley
Copyright
© 2023 Wiley Periodicals LLC.
eISSN
2573-8348
DOI
10.1002/cnr2.1808
Publisher site
See Article on Publisher Site

Abstract

INTRODUCTIONHead and neck squamous cell carcinoma (HNSCC) is a common malignancy, with an annual incidence of 500 000 cases worldwide.1,2 It is well‐established that most patients with HNSCC are diagnosed with advanced or metastatic disease at diagnosis, with a 5‐year overall survival (OS) rate of only 50% and a recurrence or metastasis rate of 30%.3,4 Indeed, once recurrence and metastasis occur, the prognosis is poor, and the median survival time is only 6 months. These findings emphasize the need to identify novel approaches to reduce recurrence and metastasis rates and improve the survival of patients with recurrent and metastatic HNSCC.In recent years, tumor immunity research has gathered momentum given the efficacy of immune checkpoint inhibitor (ICI) treatment in various tumors.5 An increasing body of evidence suggests that immune cell infiltration is essential for the efficacy of ICIs6–8 and the cause and predictor of cancer recurrence and metastasis.6,7 It is well‐established that ICIs are efficient in HNSCC8; however, predicting the efficacy of ICIs and the mechanism underlying ICI resistance remain unclear. Current evidence suggests that the tumor immunosuppressive microenvironment can inhibit immune cell function and infiltration levels.9 Accordingly, it is necessary to explore the mechanisms underlying the inhibition of immune cell infiltration in HNSCC.The rapid advent of bioinformatics has enhanced our understanding of the biological functions of tumors, providing the foothold for the discovery and verification of tumor biomarkers. There is ample evidence showing that immune‐related genes in HNSCC can predict tumor prognosis and are potential biomarkers.10,11 However, these genes are known to be inversely related to immunity. Importantly, bioinformatic analysis provides a novel approach to uncover genes related to immune cell infiltration and explore their biological functions and their potential clinical translational value as diagnostic and prognostic markers.Herein, weighted correlation network analysis (WGCNA) was conducted to identify highly correlated gene clusters using module clustering.12 We sought to identify specific immune‐related genes associated with the development or prognosis of HNSCC from The Cancer Genome Atlas (TCGA). Furthermore, we explored key genes related to the occurrence of HNSCC, the biological process involved, and their prognostic value in this patient population.MATERIALS AND METHODSData collection and preprocessingThe TCGA (https://tcgadata.nci.nih.gov/tcga/) HNSCC transcriptomic and clinical data, and two Gene Expression Omnibus (GEO) datasets (GSE65858 and GSE41613) were downloaded. A total of 528 HNSCC patients were identified in the TCGA database, consisting predominantly of males (73.11%) with no previous history of cancer (94.89%). We used DESeq213 to analyze the differential expression of genes in the transcriptomic data of TCGA HNSCC database. The screening criteria for differentially expressed genes (DEGs) included: |log2 (fold change)|>2 and p < .05.In addition, the expression levels of key genes of two HNSCC datasets from the Oncomine database (Ye Head‐Neck,14 and Peng Head‐Neck),15 were included in our analyses.Weighted correlation network analysisWGCNA16 was used to identify co‐expressed gene modules in DEGs and convert the adjacency matrix into a topological overlap matrix (TOM). According to the TOM‐based dissimilarity, DEGs were finally divided into 13 different modules. Using the InnateDB17 database, a database of genes involved in innate immunity, we identified 824 immune‐related genes, of which 72 exhibited significant changes in HNSCC. Furthermore, genes from the brown module were selected as the gene set for subsequent analysis since it contained the most significant number of immune‐related DEGs (n = 223).Gene ontology function enrichment analysisMetascape18 (https://metascape.org/), an online analysis tool, was used to perform Gene ontology (GO) function enrichment analysis (Biological Process, Molecular Function, Cellular Component) on DEGs, with the parameters set to min overlap = 3, p‐value cutoff = .01, and min enrichment = 1.5.19 Specifically, the immunological signature (Immunologic Signature) was selected for functional enrichment analysis, and the parameters remain unchanged.Gene set enrichment analysisThe Gene set enrichment analysis (GSEA)20 method was used to screen the biological state or process related to the occurrence of HNSC. We selected “h.all.v7.1.symbols.gmt” as the reference set, the number of permutations was set to 1000, and the permutation type set to the selected gene set. Gene sets with an FDR of <0.05 were significantly enriched, and the Normalized Enrichment Score (NES) was used to indicate the degree of enrichment.Protein–protein interaction network constructionWe obtained the protein–protein interactions (PPI) (score >0.4) of 223 immune‐related DEGs using STRING database21 (https://string-db.org/), and the main interaction network was visualized using Cytoscape.22 The top 10 key genes were obtained using CytoHubba, and the interaction network was drawn.Tumor immune infiltration analysisWe analyzed the gene expression data of HNSCC samples from TCGA using TIMER2.023 and analyzed the correlation between KRT4, KRT78, KRT13, and SPRR3 expression and immune cell (CD8+ T cells and macrophages) infiltration levels.ImmunohistochemistryThis study was approved by the Review Board Committee of Fujian Cancer Hospital, Fujian, China. The Review Board Committee granted a waiver of informed consent for this study. The informed consent Tissue samples from cases of pathologically diagnosed HNSCC (n = 70) treated at the Fujian Cancer Hospital from 2008 to 2017, and normal tissue samples (n = 10) were included in the analysis. The clinicopathological features of HNSCC and normal samples are shown in Table 3. The samples were fixed in formaldehyde and processed with heat‐mediated antigen retrieval in citrate buffer (pH = 6). The samples were then blocked and incubated with the following primary antibodies: rabbit polyclonal anti‐KRT13 (1:1000, Cat No. A116809, SIGMA Life Science, USA), rabbit polyclonal anti‐KRT78 (1:1000, 000015851, SIGMA Life Science, USA), and rabbit polyclonal anti‐SPRR3 (1:200, Cat No. ab218131, Abcam, USA) at 4°C overnight. The ElivisionTM plus Polyer HP (Mouse/Rabbit) immunohistochemistry (IHC) Kit (Cat. KIT‐9901, MXB biotechnologies, China) was used for IHC detection. Two independent pathologists, blinded to the clinicopathological data, evaluated the immunohistochemical score.Statistical analysisThe Wilcoxon rank‐sum test and Wilcoxon signed‐rank test were used to analyze the expression of key genes in tumor and normal tissues. The relationship between clinicopathological features and KRT4, KRT78, KRT13, and SPRR3 expression was assessed by the Kruskal–Wallis test, Wilcoxon signed‐rank test, and chi‐squared test. Survival curves were drawn using the Kaplan–Meier method, and differences between groups were assessed using the log‐rank test. OS was defined as the time of diagnosis until the date of death from any cause or last follow‐up. Progression‐free survival (PFS) was defined as the time of diagnosis until the date of disease progression, death or last follow‐up. A p‐value <.05 (two‐sided) was statistically significant. Statistical analyses were carried out using R (version 3.6.1) and SPSS (version 24.0).RESULTSIdentification and functional enrichment analysis of DEGs related to the occurrence of HNSCCTo explore the relationship between HNSCC development and immune infiltration, we conducted a comprehensive integrated analysis of transcriptomic and clinical data of HNSCC patients from TCGA (Figure 1, Table 1). First, differential expression analysis showed that 1869 and 1578 genes were significantly upregulated and downregulated in HNSCC, respectively (Figure 2A, Table S1).1FIGUREFlow chart of overall analysis1TABLEClinical characteristics of patients in the The Cancer Genome Atlas databaseClinical characteristicsNumber (%)GenderMale386 (73.11)Female142 (26.89)Prior malignancyYes27 (5.11)No501 (94.89)Vital statusAlive304 (57.58)Dead224 (42.42)AJCC clinical TT137 (7.01)T2152 (28.79)T3139 (26.33)T4184 (34.85)Not reported16 (3.03)AJCC clinical NN0246 (46.59)N185 (16.10)N2166 (31.44)N39 (1.70)Not reported22 (4.17)AJCC clinical MM0496 (93.94)M16 (1.14)Not reported26 (4.92)AJCC clinical stageI21 (3.98)II99 (18.75)III107 (20.27)IV287 (54.36)Not reported14 (2.65)2FIGUREIdentification and functional enrichment analysis of differentially expressed genes in head and neck squamous cell carcinoma (HNSCC)GO annotation demonstrated that these DEGs were mainly enriched in biological processes, including muscle contraction, myofibril assembly, and muscle structural development (Figure 2B). GO terms for molecular functions consisted of structural molecular activity, muscle structural components, and receptor ligand activity (Figure 2C), and the related cellular components included contractile fibers, extracellular matrix, and A‐bands (Figure 2D). Meanwhile, enrichment analysis of the immune characteristics of these DEGs showed that they were involved in the regulation of CD4+ T cells, macrophages, and CD8+ T cells (Figure 2E, Table S2). Furthermore, GSEA showed that biological processes that were significantly enriched in normal tissues included oxidative phosphorylation (Figure 3A), adipogenesis (Figure 3B), myogenesis (Figure 3C), bile acid metabolism (Figure 3D), fatty acid metabolism (Figure 3E), early estrogen response (Figure 3F), peroxisomes, (Figure 3G) and xenobiotic metabolism (Figure 3H). The biological processes significantly enriched in HNSCC encompassed E2F transcription factor targeting (Figure 3I), G2M checkpoint (Figure 3J), interferon α response (Figure 3K), epithelial‐mesenchymal transition (Figure 3L), interferon γ response (Figure 3M), MYC targeting (Figure 3N), DNA repair (Figure 3O) and mitotic spindle (Figure 3P; Table S3).3FIGUREGene set enrichment analysis (GSEA) of differentially expressed genes (DEGs) related to head and neck squamous cell carcinomaCo‐expression analysis of DEGs and PPI network in HNSCCTo identify gene sets related to the occurrence of HNSCC, we performed WGCNA based on transcriptomic data. Using a soft‐threshold of 4 (Figure 4A,B), WGCNA identified 13 modules (Figure 4C–E). Compared with the other modules, the brown module contained the most significant number of immune‐related genes (n = 7; Table S4). Two hundred and twenty‐three genes were enriched in the brown module, and many interactions were identified between the proteins encoded by these genes (Figure 5A,B). Proteins encoded by PPL, SCEL, KRT4, KRT24, KRT78, KRT13, SPRR3, TGM3, CRCT1, and CRNN were core components in the PPI network (Figure 5C).4FIGURECo‐expression analysis of differentially expressed genes (DEGs) in head and neck squamous cell carcinoma5FIGUREProtein–protein interaction (PPI) network of immune‐related differentially expressed genes (DEGs).Differential expression, immune cell infiltration, and survival analysis of key genesThe expression of key DEGs screened from normal and HNSCC tissues was analyzed. We found that their expression levels in the HNSCC group were significantly lower than in the normal group (Figure 6A–J, left). In the meantime, we found that except for KRT24 and CRCT1 (AUCs <0.75), the other eight genes yielded a high predictive performance for the occurrence of HNSCC, with AUC values of 0.828 (PPL), 0.804 (SECL), 0.854 (KRT4), 0.799 (KRT78), 0.824 (KRT13), 0.797 (SPRR3), 0.822 (TGM3), and 0.852 (CRNN; Figure 6A–J, right). Further analyses of immune cell infiltration of these eight genes showed that the expression levels of KRT4, KRT78, KRT13, and SPRR3 were significantly correlated with the infiltration levels of CD8+ T cells and macrophages in HNSCC (Figure 7). Then, survival analyses of these four genes were conducted, including OS and DFS, and results showed that except for KRT4, the other three genes were related to OS and DFS in HNSCC (Figure 8A). KRT78 expression was significantly related to OS and DFI (Figure 8B). KRT13 and SPRR3 expression levels were significantly related to DFI (Figure 8C) and OS (Figure 8D), respectively. Finally, KRT4, KRT78, and SPRR3 were selected as our candidate genes.6FIGUREThe RNA expression of key genes in head and neck squamous cell carcinoma7FIGUREAccurate prediction of the correlation between key genes and immune cell infiltration in head and neck squamous cell carcinoma8FIGURESurvival analysis of key genes in head and neck squamous cell carcinomaCorrelation analysis of key gene expression and clinical characteristicsA low correlation was found between KRT78 and SPRR3 expression levels and clinical T stage. However, the expression levels of these two genes in T4 HNSCC patients were significantly lower than in T3 HNSCC patients (Figure 9A,B), while no correlation was found between KRT13 expression and T stage (Figure 9C). The expression levels of these three candidate genes showed a decreasing trend as the N stage increased (Figure 9D–F). Similarly, KRT78, KRT13, and SPRR3 expression levels showed a downward trend with increasing HNSCC clinical stage (Figure 9G–I). These results suggested that the three candidate genes were closely related to HNSCC development (Table 2).9FIGURECorrelation between the expression level of candidate genes and the clinical stage of head and neck squamous cell carcinoma2TABLEDemographic and clinical characteristics of head and neck squamous cell carcinoma (HNSCC) patients in GSE41613 and GSE65858CharactersGSE41613GSE65858n (%)n (%)GenderFemale31 (32.0)43 (17.4)Male66 (68.0)223 (82.6)Age (year)<6050 (51.5)153 (56.7)≥6047 (48.5)117 (43.3)HPVNegative97 (100)197 (73.0)Positive0 (0)73 (27.0T classificationT1NA115(42.6)T3NA155 (57.4)N classificationN0‐N1NA126 (46.7)N2‐N3NA144 (53.3)Clinical stageI‐II41 (42.3)55 (20.4)III‐IV56 (57.7)215 (79.6)Comprehensive treatmentNo44 (45.4)Yes53 (54.6)External validation of candidate genesWe used Oncomine datasets (Ye Head‐Neck, and Peng Head‐Neck) to verify the expression levels of KRT13, KRT78, and SPRR3 in HNSCC tumor and normal tissues. Analysis of the Ye Head‐Neck (Figure 10A–C) and Peng Head‐Neck (Figure 10D–F) datasets showed that the expression levels of KRT13 (p < .001), KRT78 (p < .001), and SPRR3 (p < .001) in HNSCC were significantly lower than in normal tissues.10FIGUREVerification of candidate gene mRNA expression levels in Oncomine datasetsTo further analyze differences in the expression of KRT13, KRT78, and SPRR3 genes in HNSCC and normal tissues, we used a microarray of 10 normal tissues and 70 HNSCC tissues to verify the protein expression of the above genes. The results showed that the protein expression levels of KRT13 (p = .042, Figure 11A), KRT78 (p < .001, Figure 11B), and SPRR3 (p = .022, Figure 11C) in HNSCC tissues were significantly lower than in normal tissues (Table 3).11FIGUREValidation the protein expression of KRT13, KRT78, and SPRR3 in head and neck squamous cell carcinoma tissues and normal tissues by immunohistochemistry3TABLEKRT13, KRT78, and SPRR3 protein expression of head and neck squamous cell carcinoma (HNSCC) patients (n = 70) and normal tissue (n = 10) in validation cohort of protein expressionCharactersLevelHNSCCNormal tissuepGenderFemale5 (7.1%)2 (20.0%)0.210Male65 (92.9%)7 (80.0%)Age<60 years32 (45.7%)6 (60.0%)0.320≥60 years38 (44.4%)4 (40.0%)T classificationT1‐T230 (42.9%)T3‐T440 (57.1%)N classificationN0‐155 (78.6%)N2‐315 (21.4%)Clinical stageStage I–II26 (37.1%)Stage III–IV44 (62.9%)Histologic gradeG129 (59.3%)G221 (22.2%)G320 (18.5%)Anatomic siteOral cavity15 (29.6%)3 (30.0%)0.385Oropharynx4 (3.7%)2 (20.0%)Hypopharynx4 (3.7%)1 (10.0%)Larynx47 (63.0%)5 (50.0%)KRT13Negative29 (41.4)00.042Week12 (17.1)2 (20)Moderate16 (22.9)3 (30)Strong13 (18.6)5 (50)KRT78Negative40 (57.1)0 (0)<0.001Week17 (24.3)3 (30.0)Moderate8 (12.3)1 (10.0)Strong4 (5.7)6 (60.0)SPRR3Negative18 (25.7)00.022Week33(47.1)3 (30.0)Moderate12 (17.1)3 (30.0)Strong7 (10.0)4 (40.0)After verifying the expression differences of KRT13, KRT78, and SPRR3 in HNSCC and normal tissues, we used GEO datasets (GSE65858 and GSE41613) to analyze their prognostic value (Table 2). In the GSE65858 dataset, survival analysis demonstrated that patients with low KRT13 (p = .044), KRT78 (p = .0086) and SPRR3 (p = .017) expression had worse OS than those with high expression (Figure 12A–C). Furthermore, analysis of the GSE41613 dataset showed that patients with low KRT78 (p = .005) and SPRR3 (p = .02) expression had worse OS than those with low expression, while no significant difference in OS was found in patients with low and high expression of KRT13 (p = .96; Figure 11D–F).12FIGUREValidation the prognostic value of KRT13, KRT78, and SPRR3 mRNA expression in head and neck squamous cell carcinoma from GSE65858 and GSE41613DISCUSSIONThe present study explored the relationship between HNSCC development and immune infiltration via a comprehensive integrated analysis of transcriptomic and clinical data of HNSCC patients in public databases. Ten key genes, namely PPL, SCEL, KRT4, KRT24, KRT78, KRT13, SPRR3, TGM3, CRCT1, and CRNN, were screened by differential expression and PPI network analyses. We found that KRT4, KRT78, and SPRR3 were downregulated in HNSCC and negatively correlated with immune cell infiltration. Furthermore, survival analyses demonstrated that these genes were negatively correlated with OS and PFS in HNSCC. Analysis of the Oncomine and GEO datasets validated that the KRT13, KRT78, and SPRR3 expression levels were downregulated in HNSCC and negatively correlated with the prognosis of patients with HNSCC. Finally, IHC showed that KRT13, KRT78, and SPRR3 protein expression levels were downregulated in HNSCC compared to normal tissues.Moreover, we demonstrated that 1869 genes and 1578 genes were significantly upregulated and downregulated in HNSCC. Enrichment analysis showed that these DEGs regulate CD4+ T cells, macrophages, and CD8+ T cells. Furthermore, GSEA analysis showed that interferon α and γ reactions were also enriched in HNSCC, suggesting that the immune microenvironment of HNSCC was significantly different from that of normal tissues. It is widely acknowledged that the microenvironment is an essential factor leading to the occurrence and development of tumors,23 affecting the efficacy of immunotherapy.24,25 Importantly, WGCNA can harness the information of thousands or tens of thousands of the greatest gene expression changes or all genes to identify the gene set of interest and perform significant association analysis with the phenotype.26 In the present study, WGCNA showed that out of the 13 modules, the brown module contained the most immune‐related genes (n = 7).PPI network analysis showed that PPL, SCEL, KRT4, KRT24, KRT78, KRT13, SPRR3, TGM3, CRCT1, and CRNN were core components downregulated in HNSCC. KRT4, KRT13, KRT78, and SPRR3 were negatively correlated with CD8+ T cells and macrophage infiltration, suggesting that they mediated the formation of an inhibitory immune micro‐environment.Interestingly, these four key genes have been documented in keratinization, suggesting that keratinization not only reflects the cell differentiation level but also determines the immune state.Besides, we also revealed that the expression levels of KRT4, KRT78, KRT13, and SPRR3 in TCGA‐HNSCC negatively correlated with immune cell infiltration and influenced patient prognosis. Analysis of the two Oncomine datasets also confirmed that KRT13, KRT78, and SPRP3 expression levels were significantly lower in HNSCC than in normal tissues. In addition to the significant reduction in mRNA expression levels of KRT13, KRT78, and SPRR3 in HNSCC, we also documented a significant decline in the protein expression levels of KRT13, KRT78, and SPRR3 in HNSCC using IHC. This means that these 3 genes are related to the immune infiltration and prognosis of NPC, and they are potential immune‐related genes. Many previous studies had found that multiple immune‐related genes were associated with the prognosis of NPC. These genes were known immune‐related genes, such as CCR5, CD3E, CD4, SFRP4, SFRP4, CPXM1, and COL5A1 CCR6, CCL22, ROBO1, DKK1 and PDGFA, and so forth.11,27,28 In our study, WGCNA was used to lock immune‐related modules and then protein–protein interaction network was used to screen out potential immune‐related genes. The subsequent analysis also found that KRT13, KRT78, and SPRP3 were negatively correlated with immune cell infiltration. KRT13, KRT78, and SPRP3 are newly discovered immune‐related genes in HNSCC in our study.Herein, we found that patients with low expression of KRT13, KRT78, and SPRR3 in the TCGA‐HNSCC dataset had worse OS and PFI than those with high expression. This finding was also externally validated in two public GEO datasets (GSE65858 and GSE41613), where the low expression of KRT78 and SPRR3 in HNSCC patients was associated with poor OS. Our results indicated that the main biological function of Keratin 78 and SPRP3 is to participate in developmental biology and cell keratinization. Hence, these two key genes can inhibit the occurrence and progression of HNSCC acting as tumor suppressor genes. Although preliminary studies have shown that KRT7829 and SPRR330 are involved in cytokeratinization, their role in tumors remains unknown. Current evidence suggests that SPRR3 expression is low in HNSCC26 and esophageal cancer,31 while other studies have substantiated that SPRR3 is a prognostic factor of esophageal cancer and a sensitive marker of chemotherapy and radiotherapy.32–34 At present, few reports have shown that KRT78 is differentially expressed in melanoma,35 colon adenocarcinoma,36 cervical cancer,37 urinary bladder cancer,38 and HNSCC,39 and may be related to tumor prognosis.37,38 However, the biological function of KRT78 in tumors remains primarily understudied, warranting further in‐depth studies.Although this study enhanced our understanding of the involvement of KRT78 and SPRR3 in HNSCC, there were some limitations and shortcomings. First, the detailed mechanisms underlying the KRT78‐ and SPRR3‐mediated decrease in immune cell infiltration and biological functions in HNSCC were not explored. Indeed, in vivo and in vitro experiments are needed to explore and verify their biological functions. Although our study confirmed that KRT13, KRT78, and SPRR3 protein expression levels in HNSCC were significantly downregulated, the prognostic value of protein expression could not be analyzed due to the lack of survival data. It is necessary to conduct a perspective study to demonstrate the prognostic value of KRT13, KRT78, and SPRR3 protein expression.To conclude, we identified four hitherto unrecognized key genes, KRT4, KRT78, KRT13, and SPRR3, related to the occurrence and development of HNSCC and positively correlated with immune cell infiltration. KRT78 and SPRR3, novel immune‐related genes, might serve as diagnostic and prognostic biomarkers of HNSCC. Nonetheless, the detailed biological functions and clinical value of KRT78 and SPRR3 in HNSCC need further exploration.AUTHOR CONTRIBUTIONSQiaojuan Guo: Conceptualization (equal); investigation (equal); writing – original draft (equal); writing – review and editing (equal). Tianzhu Lu: Formal analysis (equal); writing – original draft (equal). Hanchuan Xu: Formal analysis (supporting); investigation (supporting); writing – original draft (supporting). Qingfeng Luo: Data curation (supporting); formal analysis (supporting). Zhilang Liu: Data curation (supporting); formal analysis (supporting). Sicong Jiang: Date curation(supporting); formal analysis (supporting). Shao jun Lin: Writing – original draft (supporting). Jianji Pan: Writing – original draft (supporting); writing – review and editing (supporting). Mengyao Lin: Formal analysis (supporting); writing – original draft (supporting); writing – review and editing (supporting). Fang Guo: Conceptualization (lead); writing – original draft (lead); writing – review and editing (lead).ACKNOWLEDGMENTSThis work was supported by Immuno Radiotherapy research fund project of Radiation Oncology Branch of Chinese Medical Association (No. Z‐2017‐24‐2020); National Natural Science Foundation of China (Nos. 81860664 and 81 660 453); Fujian Provincial Natural Science Foundation of China (No. 2019J01194); NHC Key Laboratory of Personalized Diagnosis and Treatment of Nasopharyngeal Carcinoma. The authors thank editage (http://www.fabiao@editage.cn) for editing this manuscript.CONFLICT OF INTEREST STATEMENTThe authors declare no conflict of interest.DATA AVAILABILITY STATEMENTThe datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.ETHICS STATEMENTThis study was approved by the Hospital Review Board of Fujian Cancer Hospital, Fujian, China. The Review Board Committee granted a waiver of informed consent for this study.REFERENCESSiegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin. 2018;68(1):7‐30.Belcher R, Hayes K, Fedewa S, Chen AY. Current treatment of head and neck squamous cell cancer. J Surg Oncol. 2014;110(5):551‐574.Thariat J, Vignot S, Lapierre A, et al. Integrating genomics in head and neck cancer treatment: promises and pitfalls. Crit Rev Oncol Hematol. 2015;95(3):397‐406.Mehanna H, Wong WL, Dunn J. Management of Advanced Head and Neck Cancer. N Engl J Med. 2016;375(5):492‐493.de Miguel M, Calvo E. Clinical challenges of immune checkpoint inhibitors. Cancer Cell. 2020;38(3):326‐333.Pagès F, Mlecnik B, Marliot F, et al. International validation of the consensus Immunoscore for the classification of colon cancer: a prognostic and accuracy study. Lancet (London, England). 2018;391(10135):2128‐2139.Xu X, Tan Y, Qian Y, et al. Clinicopathologic and prognostic significance of tumor‐infiltrating CD8+ T cells in patients with hepatocellular carcinoma: a meta‐analysis. Medicine. 2019;98(2):e13923.Borel C, Jung AC, Burgy M. Immunotherapy breakthroughs in the treatment of recurrent or metastatic head and neck squamous cell carcinoma. Cancer. 2020;12(9):2691.Lim AR, Rathmell WK, Rathmell JC. The tumor microenvironment as a metabolic barrier to effector T cells and immunotherapy. Elife. 2020;9:e55185.Zhou D, Wang J, Wang J, Liu X. Profiles of immune cell infiltration and immune‐related genes in the tumor microenvironment of HNSCC with or without HPV infection. Am J Transl Res. 2021;13(4):2163‐2180.Chen Y, Li ZY, Zhou GQ, Sun Y. An immune‐related gene prognostic index for head and neck squamous cell carcinoma. Clin Cancer Res. 2021;27(1):330‐341.Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA‐seq data with DESeq2. Genome Biol. 2014;15(12):550.Ye H, Yu T, Temam S, et al. Transcriptomic dissection of tongue squamous cell carcinoma. BMC Genomics. 2008;9:69.Peng S, Creighton CJ, Zhang Y, et al. Tumor grafts derived from patients with head and neck squamous carcinoma authentically maintain the molecular and histologic characteristics of human cancers. J Transl Med. 2013;11:198.Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106.Breuer K, Foroushani AK, Laird MR, et al. InnateDB: systems biology of innate immunity and beyond‐recent updates and continuing curation. Nucleic Acids Res. 2013;41:D1228‐33.Zhou Y, Zhou B, Pache L, et al. Metascape provides a biologist‐oriented resource for the analysis of systems‐level datasets. Nat Commun. 2019;10(1):1523.Tuo Z, Zhang Y, Wang X, et al. RUNX1 is a promising prognostic biomarker and related to immune infiltrates of cancer‐associated fibroblasts in human cancers. BMC Cancer. 2022;22(1):523.Subramanian A, Kuehn H, Gould J, Tamayo P, Mesirov JP. GSEA‐P: a desktop application for gene set enrichment analysis. Bioinformatics (Oxford, England). 2007;23(23):3251‐3253.von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258‐261.Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498‐2504.Li T, Fan J, Wang B, et al. TIMER: a web server for comprehensive analysis of tumor‐infiltrating immune cells. Cancer Res. 2017;77(21):e108‐e110.Giannone G, Ghisoni E, Genta S, et al. Immuno‐metabolism and microenvironment in cancer: key players for immunotherapy. Int J Mol Sci. 2020;21(12):4414.Bader JE, Voss K, Rathmell JC. Targeting metabolism to improve the tumor microenvironment for cancer immunotherapy. Mol Cell. 2020;78(6):1019‐1033.Yu L, Yang Z, Liu Y, et al. Identification of SPRR3 as a novel diagnostic/prognostic biomarker for oral squamous cell carcinoma via RNA sequencing and bioinformatic analyses. PeerJ. 2020;8:e9393.She Y, Kong X, Ge Y, et al. Immune‐related gene signature for predicting the prognosis of head and neck squamous cell carcinoma. Cancer Cell Int. 2020;20:22.Chengcheng L, Wenwen Q, Ningyue G, et al. Identification of the immune‐related genes in tumor microenvironment that associated with the recurrence of head and neck squamous cell carcinoma. Front Cell Dev Biol. 2021;9:723721.Rogers MA, Edler L, Winter H, Langbein L, Beckmann I, Schweizer J. Characterization of new members of the human type II keratin gene family and a general evaluation of the keratin gene domain on chromosome 12q13.13. J Invest Dermatol. 2005;124(3):536‐544.Koizumi H, Kartasova T, Tanaka H, Ohkawara A, Kuroki T. Differentiation‐associated localization of small proline‐rich protein in normal and diseased human skin. Br J Dermatol. 1996;134(4):686‐692.Souza‐Santos PT, de Oliveira DS, Bernardo V, Lima SC, Rapozo DC, et al. Quantitative evaluation of SPRR3 expression in esophageal squamous cell carcinoma by qPCR and its potential use as a biomarker. Exp Mol Pathol. 2011;91(2):584‐589.Luo A, Chen H, Ding F, et al. Small proline‐rich repeat protein 3 enhances the sensitivity of esophageal cancer cells in response to DNA damage‐induced apoptosis. Mol Oncol. 2013;7(5):955‐967.Luthra R, Wu TT, Luthra MG, et al. Gene expression profiling of localized esophageal carcinomas: association with pathologic response to preoperative chemoradiation. J Clin Oncol. 2006;24(2):259‐267.Zhang Y, Feng YB, Shen XM, et al. Exogenous expression of Esophagin/SPRR3 attenuates the tumorigenicity of esophageal squamous cell carcinoma cells via promoting apoptosis. Int J Cancer. 2008;122(2):260‐266.Wang LX, Li Y, Chen GZ. Network‐based co‐expression analysis for exploring the potential diagnostic biomarkers of metastatic melanoma. PloS One. 2018;13(1):e0190447.Wang H, Liu J, Li J, et al. Identification of gene modules and hub genes in colon adenocarcinoma associated with pathological stage based on WGCNA analysis. Cancer Genet. 2020;242:1‐7.Varga N, Mózes J, Keegan H, et al. The value of a novel panel of cervical cancer biomarkers for triage of HPV positive patients and for detecting disease progression. Pathol Oncol Res. 2017;23(2):295‐305.Chen S, Zhang N, Shao J, Wang T, Wang X. A novel gene signature combination improves the prediction of overall survival in urinary bladder cancer. J Cancer. 2019;10(23):5744‐5753.Zhang X, Feng H, Li Z, et al. Application of weighted gene co‐expression network analysis to identify key modules and hub genes in oral squamous cell carcinoma tumorigenesis. Onco Targets Ther. 2018;11:6001‐6021.

Journal

Cancer Reports – Wiley

Published: May 1, 2023

Keywords: biomarkers; differentially expressed genes; head and neck squamous cell carcinoma; PPI network analysis; WGCNA

There are no references for this article.