Online ISSN: 1945-4589
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 15, Issue 19 pp 10347—10369

Mitochondrial-related genes markers that predict survival in patients with head and neck squamous cell carcinoma affect immunomodulation through hypoxia, glycolysis, and angiogenesis pathways

Zhonghua Li2, , Haoxi Cai3, , Jinyang Zheng4, , Xun Chen5, , Guancheng Liu6, , Yunxia Lv7, , Hui Ye1, , Gengming Cai1,8,9, ,

  • 1 Haicang Hospital Affiliated of Xiamen Medical College, Xiamen 361026, China
  • 2 Department of Otolaryngology Head and Neck Surgery, Quanzhou First Hospital Affiliated to Fujian Medical University, Quanzhou 362000, China
  • 3 School of Stomatology, Ningxia Medical University, Yinchuan 750004, China
  • 4 Department of Pathology, Quanzhou First Hospital Affiliated to Fujian Medical University, Quanzhou 362000, China
  • 5 Department of Oral Surgery, Quanzhou First Hospital Affiliated to Fujian Medical University, Quanzhou 362000, China
  • 6 Department of Otolaryngology Head and Neck Surgery, The Hospital Affiliated of Guilin Medical College, Guilin 541000, China
  • 7 Department of Thyroid Surgery, The Second Affiliated Hospital to Nanchang University, Nanchang 330006, China
  • 8 The School of Clinical Medicine, Fujian Medical University, Fuzhou 361026, China
  • 9 The Graduate School of Fujian Medical University, Fuzhou 361026, China

Received: April 14, 2023       Accepted: September 8, 2023       Published: October 4, 2023      

https://doi.org/10.18632/aging.205081
How to Cite

Copyright: © 2023 Li et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Mitochondria play a crucial role in the occurrence and development of tumors. We used mitochondria-related genes for consistent clustering to identify three stable molecular subtypes of head and neck squamous cell carcinoma (HNSCC) with different prognoses, mutations, and immune characteristics. Significant differences were observed in clinical characteristics, immune microenvironment, immune cell infiltration, and immune cell scores. TP53 was the most significantly mutated; cell cycle-related pathways and tumorigenesis-related pathways were activated in different subtypes. Risk modeling was conducted using a multifactor stepwise regression method, and nine genes were identified as mitochondria-related genes affecting prognosis (DKK1, EFNB2, ITGA5, AREG, EPHX3, CHGB, P4HA1, CCND1, and JCHAIN). Risk score calculations revealed significant differences in prognosis, immune cell scores, immune cell infiltration, and responses to conventional chemotherapy drugs. Glycolysis, angiogenesis, hypoxia, and tumor-related pathways were positively correlated with the RiskScore. Clinical samples were subjected to qPCR to validate the results. In this work, we constructed a prognostic model based on the mitochondrial correlation score, which well reflects the risk and positive factors for the prognosis of patients with HNSCC. This model can be used to guide individualized adjuvant and immunotherapy in patients with HNSCC.

Introduction

Predominantly, head and neck squamous cell carcinoma (HNSCC) manifests as the primary histological variant of head and neck cancers, emerging as the world’s sixth most ubiquitous malignancy, encompassing over 90% of cases. Factors such as genetic predispositions, interaction with tobacco carcinogens, immoderate alcohol indulgence, and HPV infections are acknowledged determinants for HNSCC. Regrettably, the vast majority of HNSCC cases are discerned in their advanced phases. Even with significant progress in diagnostic and therapeutic strategies, the quintessential five-year survival trajectory remains unyielding at 50%. Consequently, the pressing imperative is to unearth steadfast molecular markers to optimize HNSCC clinical interventions [1].

Mitochondria are organelles with notable, metabolic activities and are considered signaling hubs with biosynthetic, bioenergetic, and signaling functions responsible for coordinating key biological pathways [2]. Mitochondria are associated with various diseases, namely cardiovascular, neurological, and metabolic disorders. Furthermore, mitochondria can influence all the processes involved in tumor formation and progression. As a result, they can affect all tumorigenic processes by regulating the metabolic, oxidative, and apoptotic processes in cancer cells. Mitochondrial DNA defects in mitochondrial functional defects have been reported in several cancers. Beyond their pivotal bioenergetic roles, mitochondria underpin tumor anabolism, orchestrate redox and calcium equilibrium, guide transcriptional governance, and regulate cellular demise. Moreover, mitochondrial dynamics hold significance in stress signaling. Cancer formation and progression are closely related to mitochondria; however, much regarding this complex relationship remains unclear [3, 4].

In this research, we harnessed mitochondria-associated genes to delineate enduring molecular subclasses through uniform clustering, juxtaposing the clinical, pathway, and immune attributes across these subcategories. Vital prognostic influencers rooted in mitochondria related genes were discerned via multi-factorial stepwise regression, culminating in the formulation of a clinical risk paradigm for such genes. To enhance the prognostic framework and fortify survival forecasts, we amalgamated the RiskScore with clinicopathological nuances, crafting a nomogram to appraise the peril faced by HNC patients, thereby paving the way for tailored therapeutic approaches.

Results

Establishment of molecular subtypes

Differential analysis was performed using the limma package on tumor and normal samples next to cancer in TCGA (|FC|>1.2 and FDR <0.05, Supplementary Table 1). Next, a one-way Cox analysis of mitochondria-related genes was performed in GSE41613 in TCGA-HNSC (P < 0.05), and the differential genes were intersected with mitochondrial genes with a significant prognosis in both datasets. The final results included 39 mitochondria-related differentially expressed genes with a significant prognosis (Figure 1A). Tumor samples expressed most of these genes (Figure 1B). Patients were then classified based on the coherent grouping of expression patterns connected to 39 key mitochondria-related genes with significant prognostic implications. The ideal cluster count was determined using the cumulative distribution function (CDF), and the CDF Delta area curve showed that selecting three clusters increased stability (Figure 1C, 1D). This final choice of k = 3 culminated in three distinct molecular subtypes: C1, C2, and C3 (Figure 1E). Delving deeper into the prognostic attributes of these subtypes revealed marked disparities in their prognostic outcomes. As depicted in Figure 1F, C1 boasted the most favorable prognosis, succeeded by C2, while C3 manifested the least favorable outcome. Employing an analogous methodology for the GSE41613 dataset, these molecular subtypes again exhibited pronounced prognostic variations (e.g., Figure 1G), resonating with the findings from the TCGA dataset (Supplementary Tables 2 and 3 detail the clustering outcomes for TCGA and GSE41613 subtypes).

Molecular typing based on mitochondrial-related genes. (A) Identification of differential mitochondria-related genes with significant prognosis; (B) Expression of 39 genes in tumor and paraneoplastic normal samples; (C) CDF curves of TCGA cohort samples; (D) CDF Delta area curves of TCGA cohort samples, Delta area curve of consensus clustering, indicating the relative change in area under the cumulative distribution function (CDF) curve for each category number k compared with k – 1. The horizontal axis represents the category number k, and the vertical axis represents the relative change in area under the CDF curve; (E) Heat map of sample clustering at consensus k = 3; (F) KM curves of the relationship among the prognosis of the three subtypes of TCGA; (G) KM curves of the prognosis of the three subtypes in the GSE41613 cohort.

Figure 1. Molecular typing based on mitochondrial-related genes. (A) Identification of differential mitochondria-related genes with significant prognosis; (B) Expression of 39 genes in tumor and paraneoplastic normal samples; (C) CDF curves of TCGA cohort samples; (D) CDF Delta area curves of TCGA cohort samples, Delta area curve of consensus clustering, indicating the relative change in area under the cumulative distribution function (CDF) curve for each category number k compared with k – 1. The horizontal axis represents the category number k, and the vertical axis represents the relative change in area under the CDF curve; (E) Heat map of sample clustering at consensus k = 3; (F) KM curves of the relationship among the prognosis of the three subtypes of TCGA; (G) KM curves of the prognosis of the three subtypes in the GSE41613 cohort.

Clinical characteristics among subtypes

We assessed the clinicopathological attributes across TCGA-HNSC molecular subtypes, uncovering pronounced distinctions in gender, clinical progression, staging, gradation, and survival outcomes among the three categories (Figure 2). Furthermore, in examining the expression of 39 mitochondria-related genes bearing prognostic relevance across distinct molecular subtypes, we discerned that the “Risk” gene was preeminent in C3, whereas the “Protective” gene exhibited pronounced elevation in C1.

Relationship between gene expression profiles and clinical characteristics among molecular subtypes.

Figure 2. Relationship between gene expression profiles and clinical characteristics among molecular subtypes.

Inter-subtype immune profile

To delineate the immune milieu across varied molecular subtypes, we scrutinized immune cell penetration employing ESTIMATE, revealing a markedly diminished “ImmuneScore” for C3 compared to its counterparts (Figure 3A), characterized by attenuated immune engagement. Utilizing the CIBERSORT methodology, we approximated the distribution of immune cell variants, identifying disparities among the molecular classifications (Figure 3B). Immune engagement evaluations were further determined via MCP-count and TIIMER, with a majority of immune cell metrics manifesting significant variations across subtypes (Figure 3C, 3D). Subsequently, the TIDE software provided insights into the prospective therapeutic implications of immunotherapy for the identified molecular categories. Elevated TIDE prognostications signify an increased propensity for immune evasion, hinting at a diminished therapeutic responsiveness to immunotherapy. Evidently, as depicted in Figure 3E, C3 from the TCGA cohort, marked by an adverse prognosis, displayed a considerably amplified TIDE metric, insinuating an escalated tendency for immune evasion and diminished therapeutic prospects from immunotherapy relative to the alternate subtypes. This is further underscored by the minimalistic response rates to immune checkpoint inhibitors in C3 (Figure 3F).

Immunological characteristics between molecular subtypes in TCGA cohort. (A) Differences in TCGA cohort ESTIMATE immune scores between molecular subtypes; (B) Differences in TCGA cohort TIDE scores and immune response between molecular subtypes CIBERSORT immune scores between molecular subtypes; (C) Differences in TCGA cohort TIMER calculated immune cell scores between molecular subtypes; (D) Differences in TCGA cohort MCP-Count immune scores between molecular subtypes; (E, F): TCGA cohort TIDE scores and immune response between molecular subtypes.

Figure 3. Immunological characteristics between molecular subtypes in TCGA cohort. (A) Differences in TCGA cohort ESTIMATE immune scores between molecular subtypes; (B) Differences in TCGA cohort TIDE scores and immune response between molecular subtypes CIBERSORT immune scores between molecular subtypes; (C) Differences in TCGA cohort TIMER calculated immune cell scores between molecular subtypes; (D) Differences in TCGA cohort MCP-Count immune scores between molecular subtypes; (E, F): TCGA cohort TIDE scores and immune response between molecular subtypes.

Inter-subtype mutation characterization/pathway characterization

Disparities in genomic modifications between these two molecular subtypes within the TCGA cohort underwent meticulous analysis. Insights into the molecular attributes of TCGA-HNSCs were extrapolated from referenced literature [5, 6]. C3 exhibited an elevated aneuploidy index, proportion altered, segment count, surpassing both C1 and C2. Conversely, C1 demonstrated notably diminished tumor purity relative to C2 and C3 (Figure 4A). Subsequently, the mutation dataset, refined by TCGA mutect2 software, was acquired. Genes surpassing a mutation frequency of three were shortlisted for salient high-frequency mutations within each subtype via the Fisher test, adhering to a selection benchmark of P < 0.05. The mutational idiosyncrasies of the paramount 15 genes within each subtype are delineated in in Figure 4B. Subsequently, the manifestation of variably activated pathways across distinct molecular subtypes was scrutinized. For pathway identification, we executed a GSEA, leveraging the h.all.v7.5.1.entrez.gmt gene set from the MSigDB database [7], with FDR <0.05 denoting noteworthy enrichment. Figure 4C portrays the outcomes of the TCGA assessment, revealing an inhibition of cell cycle-associated pathways in C1, juxtaposed with the activation of both cell cycle and tumorigenesis-associated pathways.

Genomic alterations in molecular subtypes of TCGA cohort. (A) Comparison of homologous recombination defects, aneuploidy score, fraction altered, number of segments, tumor purity in molecular subtypes of TCGA cohort differences; (B) Somatic mutations in the three molecular subtypes; (C) GSEA results among molecular subtypes of TCGA cohort.

Figure 4. Genomic alterations in molecular subtypes of TCGA cohort. (A) Comparison of homologous recombination defects, aneuploidy score, fraction altered, number of segments, tumor purity in molecular subtypes of TCGA cohort differences; (B) Somatic mutations in the three molecular subtypes; (C) GSEA results among molecular subtypes of TCGA cohort.

Identification of differential genes

In our prior analysis in this paper, three molecular subtypes were established based on prognostically significant differential mitochondria-associated genes that differed in clinical, immune, and pathway features. Subsequently, genes exhibiting differential expression amongst C1, C2, and C3 relative to other subtypes were discerned utilizing the limma package, adhering to criteria of FDR <0.05 and |log2FC| > 1. This analysis culminated in the identification of 224, 443, and 226 distinctively expressed genes in C1, C2, and C3, respectively, aggregating to an ultimate count of 662 unique genes (Supplementary Tables 4–6). Additionally, a comprehensive functional enrichment assessment of these differentially expressed genes was undertaken via the R software’s clusterProfiler package (FDR <0.05). The enriched findings from both GO and KEGG pathways pertinent to genes predominantly elevated in C3 isoforms are depicted in Figure 5 (Supplementary Table 7).

Enrichment results of GO and KEGG pathways of differentially upregulated genes in C3 subtype. (A) Results of KEGG analysis of upregulated genes in TCGA-HNSC cohort C3; (B–D) Results of GO analysis of upregulated genes in TCGA-HNSC cohort C3.

Figure 5. Enrichment results of GO and KEGG pathways of differentially upregulated genes in C3 subtype. (A) Results of KEGG analysis of upregulated genes in TCGA-HNSC cohort C3; (B–D) Results of GO analysis of upregulated genes in TCGA-HNSC cohort C3.

Risk modeling

Utilizing the coxph function from the SURVIVAL package, a univariate Cox analysis was conducted on the 662 differential genes, pinpointing 20 genes of profound prognostic significance (P < 0.001) as illustrated in Supplementary Table 8. Subsequent to this, a multifactorial stepwise regression analysis was undertaken. This analytical approach employs the AIC information criterion—a metric balancing the statistical fit against the parameter count in the model. Using the AIC method within the MASS package, the analysis commences with the most intricate model, systematically excluding variables to minimize the AIC value. A diminished value signifies an optimal model, suggesting a satisfactory fit achieved with minimal parameters. Ultimately, nine pivotal genes influencing prognosis were discerned, as exemplified in Figure 6A. The definitive model equation is presented below: RiskScore = 0.073 × DKK1 + 0.154 × EFNB2 + (−0.232 × ITGA5) + 0.064 × AREG + (−0.071 × EPHX3) + 0.069 × CHGB + 0.242 × P4HA1 + 0.103 × CCND1 + (−0.061 × JCHAIN) Subsequently, risk scores for each specimen within the TCGA cohort were derived utilizing the expression metrics of the nine distinct genes. Based on the median RiskScore, TCGA samples were segregated into high-risk and low-risk contingents, and the Kaplan-Meier curves demonstrated a pronounced disparity between these cohorts (Figure 6B). The timeROC R package facilitated an ROC analysis on the prognostic stratification by RiskScore, examining the predictive efficacy over intervals from 1 to 5 years. The model exhibited an impressive area beneath the AUC curve (Figure 6C). In an endeavor to reinforce the model’s robustness validation, analogous methodology was applied to authenticate the GSE41613 and GSE65858 datasets. These validations confirmed a commendable AUC area, with conspicuous contrasts between high-risk and low-risk classifications (Figure 6D–6G).

Determination of risk model and its KM, ROC curves. (A) Multifactor forest plot of prognostic key genes; (B) KM curve of risk model constructed for 9 genes in TCGA dataset; (C) ROC curve of the risk model in TCGA dataset; (D) KM curve of risk model constructed for 9 genes in GSE41613 dataset; (E) ROC curve of risk model in GSE41613 dataset; (F) KM curve of risk model constructed for 9 genes in GSE65858 dataset; (G) ROC curve of risk model in GSE65858 dataset.

Figure 6. Determination of risk model and its KM, ROC curves. (A) Multifactor forest plot of prognostic key genes; (B) KM curve of risk model constructed for 9 genes in TCGA dataset; (C) ROC curve of the risk model in TCGA dataset; (D) KM curve of risk model constructed for 9 genes in GSE41613 dataset; (E) ROC curve of risk model in GSE41613 dataset; (F) KM curve of risk model constructed for 9 genes in GSE65858 dataset; (G) ROC curve of risk model in GSE65858 dataset.

Nomogram

Univariate and multifactorial Cox regression assessments of RiskScore, juxtaposed with clinical attributes, ascertained that RiskScore, age, and stage stood as paramount prognostic determinants (Figure 7A, 7B). To crystallize the interplay between risk evaluation and patient survival trajectories, we melded the RiskScore with pertinent clinicopathological indicators, giving rise to an illustrative columnar depiction (Figure 7C); here, RiskScore emerged as the preeminent force in forecasting survival. The model’s prognostic fidelity underwent scrutiny via the calibration trajectory depicted in Figure 7D. Calibration curves at the 1-, 3-, and 5-year milestones were in near-congruence with the benchmark trajectory, underscoring the model’s adept prognostic prowess. The model’s robustness garnered further validation through DCA (Decisioncurve); wherein both RiskScore and the nomogram exhibited a palpable ascendancy over extreme curves (Figure 7E). As illuminated in Figure 7F, the ROC curve analysis elucidated that the nomogram, coupled with RiskScore, manifested unparalleled sensitivity and specificity in forecasting the overarching survival of head and neck carcinoma patients, eclipsing other pertinent clinical markers such as age, gender, histological caliber, and TNM clinical echelon, thus asserting their superlative prognostic merit. To elucidate the nexus between the RiskScore and TCGA-HNSC clinical attributes, we probed the variances in risk stratifications and scores vis-à-vis clinical gradations utilizing the TCGA-HNSC compendium. The representation of C3 preponderated within the high-risk cohort (Figure 8A), whilst the RiskScore ascended concomitantly with the clinical tier (Figure 8B). Additionally, juxtapositions of the RiskScore across molecular subtypes revealed C3, bearing the most ominous prognosis, also brandished the zenith RiskScore. We further juxtaposed the prognostic disparities between high and low RiskScore strata across diverse clinicopathological categorizations. The results demonstrated that the risk groupings were equally good in different clinical groupings, proving the reliability of the risk groupings (Figure 8C–8F).

Determination and survival prediction ability of nomogram. (A) Single-factor forest plot of RiskScore and clinical features; (B) Multifactor forest plot of RiskScore and clinical features; (C) RiskScore combined with clinical features column line plot; (D) Calibration curves of column line plot at 1, 3, and 5 years; (E) Decision curve of column line plot; (F) ROC curves of various clinical features for overall survival (OS) at 1, 3, and 5 years.

Figure 7. Determination and survival prediction ability of nomogram. (A) Single-factor forest plot of RiskScore and clinical features; (B) Multifactor forest plot of RiskScore and clinical features; (C) RiskScore combined with clinical features column line plot; (D) Calibration curves of column line plot at 1, 3, and 5 years; (E) Decision curve of column line plot; (F) ROC curves of various clinical features for overall survival (OS) at 1, 3, and 5 years.

Clinical characteristics between risk groups. (A) Key gene expression in relation to clinical characteristics and RiskScore; (B) Differences between RiskScore between clinicopathological subgroups in TCGA-HNSC cohort; (C) KM curves between RiskScore of the high- and low-risk groups for different stages in TCGA-HNSC cohort; (D) KM curves between RiskScore high- and low-risk groups for grades in TCGA-HNSC cohort grade; (E) KM curves between RiskScore high- and low-risk groups of age groups in TCGA-HNSC cohort; (F) KM curves between the RiskScore of high- and low-risk of different genders in TCGA-HNSC cohort.

Figure 8. Clinical characteristics between risk groups. (A) Key gene expression in relation to clinical characteristics and RiskScore; (B) Differences between RiskScore between clinicopathological subgroups in TCGA-HNSC cohort; (C) KM curves between RiskScore of the high- and low-risk groups for different stages in TCGA-HNSC cohort; (D) KM curves between RiskScore high- and low-risk groups for grades in TCGA-HNSC cohort grade; (E) KM curves between RiskScore high- and low-risk groups of age groups in TCGA-HNSC cohort; (F) KM curves between the RiskScore of high- and low-risk of different genders in TCGA-HNSC cohort.

Immune/pathway characteristics among risk subgroups

To discern variances in the immune microenvironment across risk strata, we employed the CIBERSORT algorithm to quantify immune cell prevalence within TCGA-HNSC’s high- and low-risk factions (Figure 9A). Utilizing both MCP-count and TIMER methodologies, we derived the immune infiltration indices for the TCGA-HNSC cohort, revealing pronounced disparities in most immune cell indices between risk strata (Figure 9B). The “ImmuneScore” was markedly elevated in the low-risk segment compared to its high-risk counterpart. Analyzing immunotherapeutic divergences between risk sectors within the TCGA cohort, we harnessed the TIDE software to gauge potential clinical immunotherapy outcomes for both risk sectors. Notably, the TIDE score for the high-risk faction surpassed that of the low-risk, implying an augmented propensity for immune evasion in the former, yet paradoxically suggesting a greater therapeutic benefit potential. Conversely, the amplified TIDE score in the high-risk cohort indicated increased susceptibility to immune subterfuge and a diminished probability of immunotherapeutic success (Figure 9D). We also employed the Pearson approach to deduce the correlation between RiskScore and immune score, unearthing a pronounced inverse association between RiskScore and a majority of immune cell indices (Figure 9C). Furthermore, we assessed the responsiveness to conventional chemotherapeutics across TCGA cohort subgroups, observing heightened sensitivity within the high-risk cohort compared to its low-risk counterpart (Figure 9E). To delve deeper into the nexus between RiskScore and biological functions across diverse samples, we employed the TCGA dataset, designating h.all.v7.5.1.symbols.gmt as our gene set. Thereafter, we executed single-sample GSEA (ssGSEA) analyses via the R package GSVA. For each distinct function, the ssGSEA enrichment scores were derived for every sample. Employing the rank sum test, we discerned distinct pathways between risk groups (P < 0.05) and further evaluated the interrelation between these pathways and the RiskScore. As illustrated in Figure 10A, pathways integral to glycolysis, angiogenesis, hypoxia, and specific tumorigenic processes manifested a robust positive association with the RiskScore. Through GSEA, we also analyzed prominently enriched pathways within both high- and low-risk cohorts. Adopting a threshold of NP <0.05, we pinpointed the saliently enriched pathways (Figure 10B).

Immune and pathway characteristics between different risk groups. (A) difference in CIBERSORT immune infiltration between risk subgroups in TCGA cohort; (B) difference in MCP-count, TIMER immune score between risk subgroups in TCGA cohort; (C) correlation between immune score and RiskScore in TCGA cohort; (D) difference in TIDE score and ESTIMATE immune score between risk subgroups in TCGA cohort; (E) difference in drug sensitivity (IC50) between risk subgroups in TCGA cohort.

Figure 9. Immune and pathway characteristics between different risk groups. (A) difference in CIBERSORT immune infiltration between risk subgroups in TCGA cohort; (B) difference in MCP-count, TIMER immune score between risk subgroups in TCGA cohort; (C) correlation between immune score and RiskScore in TCGA cohort; (D) difference in TIDE score and ESTIMATE immune score between risk subgroups in TCGA cohort; (E) difference in drug sensitivity (IC50) between risk subgroups in TCGA cohort.

Correlation between differential pathways and RiskScore in risk groups and their GSEA results in TCGA dataset. (A) Correlation between risk intergroup differential pathways and RiskScore in TCGA dataset; (B) GSEA results for high- and low-risk groups in TCGA dataset.

Figure 10. Correlation between differential pathways and RiskScore in risk groups and their GSEA results in TCGA dataset. (A) Correlation between risk intergroup differential pathways and RiskScore in TCGA dataset; (B) GSEA results for high- and low-risk groups in TCGA dataset.

Expression of mitochondrial prognostic-related genes in HNSCC

Prior findings indicated that CCL22, CTSG, and FGD3 expressions were auspiciously linked with HNSCC prognosis, while TPP1 exhibited an inverse correlation. Using RT-qPCR, we probed the mRNA expression levels of these four genes in both tumor and adjacent non-tumorous tissues from 10 HNSCC patients. Notably, expressions of (A) AREG (observed in 8 of 10 samples, or 80%; Figure 11A), (B) DKK1 (observed in 6 of 10 samples, or 60%; Figure 11B), and (C) EFNB2 (observed in 7 of 10 samples, or 70%; Figure 11C) in tumorous specimens were markedly elevated compared to their non-tumorous counterparts. These data underscore that the expression patterns of certain prognostic genes in HNSCC align seamlessly with our anticipations.

RT-qPCR. Expression of prognostic-specific genes in HNSCC were consistent with the predicted trend. RT-qPCR detected mRNA expression levels of (A) AREG, (B) DKK1, and (C) EFNB2 in the tumor tissue and adjacent tissue of patients with HNSCC p

Figure 11. RT-qPCR. Expression of prognostic-specific genes in HNSCC were consistent with the predicted trend. RT-qPCR detected mRNA expression levels of (A) AREG, (B) DKK1, and (C) EFNB2 in the tumor tissue and adjacent tissue of patients with HNSCC p < 0.05.

Discussion

While aberrations in mitochondrial genes frequently manifest in cancer cells, they seldom incapacitate mitochondrial energy metabolism. Instead, they recalibrate the mitochondria’s bioenergetic and biosynthetic equilibrium. Such altered states liaise with the nucleus via ‘retrograde signaling’, orchestrating a gamut of signal transduction pathways, transcriptional networks, and chromatin configurations, adeptly catering to the mitochondrial and nuclear exigencies of malignancies. Subsequently, malignant cells reconfigure neighboring stromal cells, refining the tumoral milieu. Mitochondria, pivotal in tumor anabolism, are custodians of redox balance, calcium homeostasis, transcriptional oversight, and cellular demise. Tumorigenesis, progression, and therapeutic responses are intricately intertwined with host immunological processes, many of which pivot on unimpaired mitochondrial metabolism [2, 4]. Reactive oxygen species (ROS) don pivotal roles—spanning cell death regulation, DNA repair, stem cell sustenance, metabolic shifts, and sculpting the tumor microenvironment. Intriguingly, such ROS also influence T-cell dynamics. In HNC, modulating ROS levels during chemotherapy/radiotherapy is of paramount clinical import. An array of non-oncologic pharmaceuticals, molecular compounds, traditional herbal remedies, and avant-garde interdisciplinary methodologies—encompassing photodynamics, nanotech systems, and Bio Electro-Magnetic-Energy-Regulation (BEMER) therapies—can adeptly modulate cellular ROS in HNC, potentiating the efficacy of conventional treatment modalities [8–11]. Rencelj reported posited that mitochondrial receptors incisively target pivotal metabolic entities, influencing myriad oncogenic signaling cascades and holding sway over the Warburg effect—a phenomenon vital for cancer cell proliferation. A nexus between mitochondria and hypoxia has been delineated in breast cancer and HNCs [12]. Lee et al. elucidated that the majority of HNCs, rooted in mucosal epithelial cells, retain epithelial characteristics. However, as malignancy advances, it often morphs into mesenchymal or hypofractionated phenotypes, paving the way for invasion, metastasis, and treatment resistance. Erosion of epithelial traits, courtesy of epithelial-mesenchymal transition, might predispose such tenacious malignancies to ferroptosis. Steering mitochondrial or iron metabolism can amplify intracellular ferrous iron and lipid peroxidation, rendering drug-resistant neoplasms more amenable to ferroptosis [13]. Huang and associates postulated that metabolic shifts originating from nasopharyngeal carcinoma cells might foster tumor advancement and immunosuppression via intercellular discourse with adjacent immune cells. Crafting therapeutic targets at this metabolic-immune interface could herald novel treatment avenues for nasopharyngeal carcinoma [14]. Ergo, strategic manipulation of mitochondrial metabolism emerges as a promising therapeutic frontier for HNC. In our investigation, we discerned a differential expression of mitochondrial-associated genes between tumoral and adjacent healthy tissues, pinpointed molecular isoforms through robust clustering, and spotlighted nine pivotal mitochondrial-associated genes. Among them, CCND1, located on chromosome band 11q13, codes for cell cycle protein D1, pivotal in orchestrating G1-S phase transition. Beyond its pro-proliferative prowess, it augments cellular migratory potential, curtails differentiation, impedes DNA reparative processes, and is intricately woven into the mitochondrial metabolism governing the cell cycle. Consequently, its role as an oncogene in specific neoplasms is undeniable [15]. CCND1, perceived as a mitotic cellular sentinel, can, when dysregulated, plunge the cell into cycles of disarray and dysfunction [16]. Given its frequent amplification in malignancies [17], CCND1 remains a cynosure in oncologic pursuits [18]. Dickkopf-related protein 1 (Dkk1), a protein consisting of 266 amino acids, manifests predominantly within mature tissues, including bone, placenta, prostate, spleen, and colon. It operates as a formidable adversary to the Wnt signaling cascade. DKK1 obstructs the intricate interplay between Wnt, FZD, and LRP6, culminating in the degradation of β-catenin and the subsequent dormancy of the β-catenin/T-cell-specific factor (TCF) transcriptional ensemble, thereby attenuating genes steered by T-cell factors. Clinical evaluations across a spectrum of malignancies have detected amplified concentrations of this Wnt antagonist, Dickkopf-1(DKK1), within patient sera and neoplasms, frequently serving as harbingers of an unfavorable prognosis. DKK1 exhibits prowess in modulating immune cellular dynamics and curbing the tumoral immunosuppressive milieu. Its involvement in T-cell differentiation and fostering tumor subterfuge from immune scrutiny, primarily through the proliferation of MDSCs, has thrust DKK1 into the spotlight as a prospective linchpin in cancer immunotherapy [19]. IIt acts as a catalyst in the tumorigenesis of pancreatic, esophageal, and hepatocellular carcinomas, stands as a prognostic marker in breast, gastric, and colorectal malignancies, and spurs prostate neoplastic cell growth and migration [20–26]. Additionally, DKK1 orchestrates the accrual and function of MDSCs in malignancies, is ubiquitously acknowledged as a potent tumor-affiliated antigen in multiple myeloma, and has been anointed as an innovative immunotherapy target for myeloma [27–29]. Contrarily, its role in mitigating breast neoplastic cell movement and infiltration is achieved through suppression of the β-catenin/MMP7 signaling conduit [30]. Ephrins serve as the ligands for the Eph receptor tyrosine kinase subset, with EFNB2 donning the dual hats of a receptor ligand and a signaling receptor. Its omnipresence within tumoral vasculature, coupled with its role in catalyzing tumor angiogenesis and neovascularization by anchoring vascular endothelial forerunners to neoplastic sites and its influence on lymphangiogenesis, underscores its significance. EFNB2 gene expression, fostering cellular proliferation, migration, and invasive tendencies particularly in pancreatic ductal adenocarcinoma, emerges as both a biological compass and a prognostic touchstone in a plethora of malignancies. Its expression intensity resonates with the malignancy’s severity [31–34]. Eph signaling’s intricate nexus with tumoral migration, growth, angiogenesis, apoptosis, and metastasis heralds its prospective therapeutic potential [35]. ITGA5, an affiliate of the integrin chain consortium, is an ever-present adhesive entity. These cellular sentinels are pivotal in adhesive dynamics and intracellular communication. They sculpt neovascularization patterns and wield influence over neoplastic growth, invasion, and metastatic spread via conduits like FAK and PI3K [36]. Its association spans from hepatocellular carcinomas to gliomas, including mesenchymal shifts in tongue squamous cell carcinomas [37–39]. Its prognostic worth extends to non-small cell lung neoplasms and pediatric acute myeloid leukemia [40, 41]. Augmenting this, curbing ITGA5 has been shown to bolster certain chemotherapeutic regimens’ efficacy [42]. Amphiregulin (AREG), a ligand specific to the epidermal growth factor receptor (EGFR), plays a pivotal role in the intricate orchestration of juxtacrine signaling amongst neighboring cells. In a distinct capacity, AREG is excreted, functioning as an autocrine or paracrine mediator. Its gene expression and subsequent discharge are precipitated by a myriad of agents, encompassing inflammatory lipids, cytokines, hormones, growth factors, and foreign substances [43]. By binding to EGFR, AREG triggers a slew of salient intracellular signaling pathways that delineate cellular viability, proliferation, and dynamism [44]. Exogenously secreted AREG, upon interfacing with EGFR, instigates a self-propagating feedback loop that enhances AREG transcription. Multiple intracellular conduits, such as MAPK, PI3K/Akt, STAT, and mTOR, are activated in the wake of AREG interaction, serving as the fulcrum for AREG/EGFR-directed cellular operations [45]. AREG has also been associated with oxidative stress and ferroptosis pathways, and HIF2-α regulates AREG, which are consistent with the results in the literature on AREG-induced transcriptional pathways under hypoxic conditions [46]. Studies have also shown that HIF-2α can enhance oxidative death in colon cancer cells through ferroptosis activators and DMF [47]. Areg is involved in immune regulation; can be an autocrine factor for tissue Treg, Treg expresses EGF receptors; and enhances Treg function [48]. It also enhances the differentiation of Th9 cells [49]. Mast cell-derived AREG enhances the immunosuppressive capacity of regulatory T (Treg) cells, except for AREG expressed by cancer cells in specific cases, which produce acquired resistance to anti-EGFR therapies (e.g., cetuximab), and sensitivity to cetuximab depends on high expression of both EREG and AREG [50, 51]. AREG is also associated with hepatocellular carcinoma, cholangiocarcinoma, pancreatic cancer, lung cancer, and breast cancer [52–59]. Studies have shown that AREG can promote multiple cancer models in the invasion and in regulating tumorigenesis [60, 61] and that targeted stromal-derived AREG can eliminate the AREG generated by stromal and cancer cells at the TME ecological site [50] and can be used in combination with anti-PD-1 antibodies [62].

EPHX3, an adept epoxide hydrolase, proficiently metabolizes volatile xenobiotic epoxides and modulates endogenous epoxides integral to cellular signaling. With pronounced expression in the proximal digestive tract, bone marrow, lymphatic structures, and dermal layers, EPHX3 emerges as a paramount arbiter of tumorigenesis across 13 malignancies. In the context of HNSCC, it not only epitomizes a prognostic indicator but also manifests its antineoplastic prowess by curbing tumor immune checkpoint articulation and immune cellular permeation [63]. EPHX3 hypermethylation may contribute to the development of OSCC and is associated with adenoid cystic carcinogenesis and progression [64]. Chromophobe granule B (CHGB) is one of the two major soluble proteins in the chromophobe granules of the adrenal medulla [65]. CHGB, instrumental in immune modulation, exhibits deviant gene expression across myriad tumor varieties, with its augmented expression being intrinsically linked to metastatic events. It may be a prognostic marker in neuroendocrine tumors, correlates with survival in colon cancer, correlates with malignant behavior in pheochromocytomas and abdominal paragangliomas, and correlates with poor prognosis in several squamous cell carcinomas (SCCs), namely HNSCC [66–69]. P4HA1 expression, catalyzed by HIF-1 under hypoxic conditions, orchestrates various tumorigenic pathways, including EMT, angiogenesis, invasion, inflammation, and notably, the glycolytic process. P4HA1 has an important role in the HIF-1 signaling pathway [70] and is associated with various tumors, namely HNSCC, breast cancer, colorectal cancer, and B-cell lymphoma [71–75]. JCHAIN is a polypeptide chain specific to polymeric immunoglobulins (Igs) and is thought to be responsible for linking monomeric subunits into a polymeric form [76]. It is required for the multimerization of IgM and IgA, a small polypeptide required for the transport of these Ig classes across the mucosal epithelium in a multi-Ig receptor-mediated process [77] associated with various cancers, namely ovarian, gastric, and breast cancers [78–80]. JCHAIN is linked to both inherent and adaptive resistance to radiation therapy in nasopharyngeal carcinoma patients, influencing their prognostic outcomes [81].

We established three molecular subtypes, C1-C3, based on prognostically significant differential mitochondria-related genes, which were associated with cell cycle-related and tumorigenesis-related pathways, with cell cycle-related pathways being inhibited in C1 and cell cycle-related and tumorigenesis-related pathways being activated in C3. The differential genes of the subtypes were further analyzed by one-way Cox analysis and multifactor stepwise regression analysis. In the final analysis, nine pivotal genes were discerned as instrumental in influencing prognosis. Notably, DKK1, EFNB2, AREG, CHGB, P4HA1, and CCND1 manifested as risk determinants, while ITGA5, EPHX3, and JCHAIN emerged as protective indicators. The RiskScore affirmed the model’s robustness, evident from the substantial area beneath the AUC curve, delineating marked disparities between high-risk and low-risk cohorts. The mRNA expression corroborated in clinical specimens mirrored largely congruent patterns; predominantly, there was a pronounced diminution in tumorous tissues, juxtaposed with an amplification in normal samples. The RiskScore showed a significant negative correlation with most immune cell scores. The RiskScore was positively correlated with glycolysis, angiogenesis, hypoxia, and tumor-related pathways (WNT, TFGβ, EMT) and further influenced immune regulation. These phenomena were verified using the GSE41613 and GSE65858 datasets.

In summation, we devised a risk paradigm informed by nine mitochondria-centric genes, synergizing the RiskScore with clinicopathological nuances to elevate the precision of our prognostic framework and survival foresight, which boasted notable predictive acuity. Moreover, we scrutinized variances within the immune milieu across risk categories, discerning that the RiskScore harbored a marked inverse association with the majority of immune cell evaluations. We further gauged the responsiveness of distinct subgroups to traditional chemotherapy regimens. This model stands as a salient tool in HNC, facilitating therapeutic guidance and tailoring patient-centric interventions.

Materials and Methods

Data collection and processing

Utilizing the TCGA GDC API, we procured RNA-seq data from TCGA-HNSC. Post meticulous selection, we incorporated 499 primary tumor specimens and 44 normative samples. From the GEO repository, we garnered expression metrics for GSE41613, leading to the inclusion of 97 distinguished samples. Similarly, for GSE65858, post-evaluation, 270 samples were integrated.

Source of mitochondria-related genes

The mitochondria-associated genes were obtained from the literature [82]. Mitochondria-associated genes encode mitochondria-localized proteins, namely all proteins in the mitochondrial membrane, matrix, cristae, and mitochondria-associated endoplasmic reticulum membrane. On the basis of subcellular localization, all genes were obtained from the Molecular Signature Database (MSigDB) database (http://software.broadinstitute.org/gsea/msigdb) for a set of 23 mitochondria-associated cellular component genes; 1576 genes were identified as mitochondria-associated genes (Supplementary Table 9).

Data preprocessing

From TCGA, RNA-seq data underwent a meticulous four-phase refinement:

  1. Exclusion of specimens lacking clinical follow-up data.

  2. Preservation of samples boasting survival durations exceeding 0 days.

  3. Elimination of samples devoid of status.

  4. Retain coding protein genes.

For the GEO dataset, the subsequent procedures were executed:

Upon accessing the pertinent microarray platform’s annotation details, probes were aligned with genes in accordance with this annotated data. Probes correlating with multiple genes were excised. Conversely, when several probes corresponded to a singular gene, their mean expression served as the definitive gene expression value.

Molecular subtyping of mitochondria-related genes

We constructed consistency matrices by using ConsensusClusterPlus for the clustering samples [83]. We obtained the molecular subtypes of the model by using mitochondria-related gene expression data with a significant prognosis. Utilizing the “km” algorithm with “Euclidean” as the distance metric, we executed five hundred bootstraps, each encompassing 80% of the patients from the training set. The cluster range was established between 2 and 10. The optimal classification emerged from evaluating the consistency matrix and the CDF, thereby discerning the molecular subtypes of the specimen.

Risk model

  1. Identify the differentially expressed genes between subtypes via the molecular subtypes identified previously (|log2FC|>1 and FDR <0.05).

  2. Select genes with a significant prognosis (P < 0.001).

  3. Further reduce the number of genes by using multifactor stepwise regression to obtain the prognostically significant genes associated with the mitochondrial phenotype.

  4. Conduct risk modeling. The risk score for each patient was determined using the equation: RiskScore = Σβi × Expi. In this equation, ‘I’ denotes the expression level of genes associated with mitochondrial phenotype prognostic features, and β’ represents the pertinent gene’s Cox regression coefficient. Patients were stratified into two categories—high risk and low risk based on the median risk score. Survival curves were crafted using the Kaplan-Meier approach, with the log-rank test assessing the significance of observed disparities.

GSEA

To investigate the pathways of different biological processes in different molecular subtypes, we used “GSEA” for pathway analysis: we performed gene set enrichment analysis using candidate gene sets from the HALLMARK database.

Calculation of TME cell invasion abundance

The CIBERSORT algorithm ascertained the proportional abundances of 22 immune cells within the tumor tissue. Immune infiltration was gauged via the ESTIMATE software, while the MCP Count function assessed the metrics of 10 immune cells. Conversely, the TIMER function evaluated the metrics of six distinct immune cells.

Prediction of responsiveness to immunotherapy

The TIDE algorithm, a computational stratagem forecasting ICB responsiveness via gene expression profiling, authenticated the impact of IMS on predicting clinical receptivity to ICIs [84]. This algorithm scrutinized three cellular entities impeding T-cell penetration in tumors: the M2 variant of CAF, MDSCs, and TAMs. Moreover, it appraised two distinct tumor immune evasion methodologies: the compromised efficacy of tumor-infiltrating CTLs and the counteraction of CTLs due to immunosuppressive elements.

qRT-PCR

TRIzol facilitated the extraction of total RNA from pristine human HNSCC specimens and adjacent tissues, subsequently undergoing reverse transcription to cDNA. For the scope of this investigation, the patient consented to the surgical utilization of human tissue at the Quanzhou First Hospital, affiliated with Fujian Medical University, spanning October 2021 to November 2022. This endeavor received the endorsement of the Ethics Committee of the said institution, affirming the projected blueprint Quantitative real-time PCR, standardized to actin, was executed on the ABI 7900 apparatus employing the SYBR Green RT-PCR assay. The following primers were used for PCR: AREG: 5′-CTGTCGCTCTTGATACTCG-3′ (sense), 5′-CAGAAAATGGTTCACGCT-3′ (antisense); DKK1: 5′-AACTGGGAGAAGATGGCT-3′ (sense), 5′-TCCTGGGGTGAAAGTATG-3′ (antisense); and EFNB2: 5′-ACATTCGGGGAACAACAT-3′ (sense), 5′-TTCAGCAAGAGGACCACC-3′ (antisense).

Statistical analysis

Analyses were executed with GraphPad Prism 5 and R (version 3.6.3), representing data as mean ± SD. The Student’s t-test assessed differential expression across tissues. The Cox function in R facilitated the Univariate Cox analysis, elucidating associations with HNSCC prognosis. Prognostic evaluations employed Kaplan-Meier survival curves, with the log-rank test discerning the significance of disparities. A threshold of P < 0.05 demarcated statistical significance.

Data availability

The data used to support the findings of this study have been included in this article.

Abbreviations

HNSC: Head and Neck Squamous Cancer; HPV: Human papillomavirus; TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus database; MSigDB: Molecular Signatures Database; GSEA: Gene Set Enrichment Analysis; CAF: tumor-associated fibroblasts; MDSCs: myeloid-derived suppressor cells; TAM: tumor-associated macrophages; CTLs: tumor-infiltrating cytotoxic T lymphocytes; TIDE: Tumor Immune Dysfunction and Exclusion; DCA: Decisioncurve; CNV: Copy Number Variation; SNV: Simple Nucleotide Variation; RNA-seq: Transcriptome expression data; ROC: receiver operating characteristic curve; OS: Overall Survival; ICB: Immune Checkpoint Blockade.

Author Contributions

(I) Conception and design: Gengming Cai; (II) Administrative support: Hui Ye; (III) Provision of study materials or patients: Zhonghua Li; (IV) Experiments designed and performed: Zhonghua Li, Haoxi Cai; (V) Collection and assembly of data: Zhonghua Li; (VI) Data analysis and interpretation: Haoxi Cai, Jinyang Zheng, Xun Chen, Guancheng Liu, Yunxia Lv, Zhonghua Li; (VII) Manuscript writing: All authors; (VIII) Final approval of manuscript: All authors.

Acknowledgments

We would like to thank Editage (http://www.editage.cn/) for English language editing.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Ethical Statement and Consent

The studies involving human participants were reviewed and approved by the Ethics Committee of Quanzhou First Hospital Affiliated to Fujian Medical University (No. (2021) 228). All the patients provided written informed consent to participate in the study.

Funding

This work was funded by The Natural Science Foundation of Fujian Provincial Department of Science and Technology (No. 2020J011281).

References

  • 1. Zhao X, Chen H, Qiu Y, Cui L. FAM64A promotes HNSCC tumorigenesis by mediating transcriptional autoregulation of FOXM1. Int J Oral Sci. 2022; 14:25. https://doi.org/10.1038/s41368-022-00174-4 [PubMed]
  • 2. Porporato PE, Filigheddu N, Pedro JMB, Kroemer G, Galluzzi L. Mitochondrial metabolism and cancer. Cell Res. 2018; 28:265–80. https://doi.org/10.1038/cr.2017.155 [PubMed]
  • 3. Klein K, He K, Younes AI, Barsoumian HB, Chen D, Ozgen T, Mosaffa S, Patel RR, Gu M, Novaes J, Narayanan A, Cortez MA, Welsh JW. Role of Mitochondria in Cancer Immune Evasion and Potential Therapeutic Approaches. Front Immunol. 2020; 11:573326. https://doi.org/10.3389/fimmu.2020.573326 [PubMed]
  • 4. Srinivasan S, Guha M, Kashina A, Avadhani NG. Mitochondrial dysfunction and mitochondrial dynamics-The cancer connection. Biochim Biophys Acta Bioenerg. 2017; 1858:602–14. https://doi.org/10.1016/j.bbabio.2017.01.004 [PubMed]
  • 5. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, Porta-Pardo E, Gao GF, Plaisier CL, Eddy JA, Ziv E, Culhane AC, Paull EO, et al, and Cancer Genome Atlas Research Network. The Immune Landscape of Cancer. Immunity. 2018; 48:812–30.e14. https://doi.org/10.1016/j.immuni.2018.03.023 [PubMed]
  • 6. Bagaev A, Kotlov N, Nomie K, Svekolkin V, Gafurov A, Isaeva O, Osokin N, Kozlov I, Frenkel F, Gancharova O, Almog N, Tsiper M, Ataullakhanov R, Fowler N. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell. 2021; 39:845–65.e7. https://doi.org/10.1016/j.ccell.2021.04.014 [PubMed]
  • 7. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015; 1:417–25. https://doi.org/10.1016/j.cels.2015.12.004 [PubMed]
  • 8. Chávez MD, Tse HM. Targeting Mitochondrial-Derived Reactive Oxygen Species in T Cell-Mediated Autoimmune Diseases. Front Immunol. 2021; 12:703972. https://doi.org/10.3389/fimmu.2021.703972 [PubMed]
  • 9. Huang G, Pan ST. ROS-Mediated Therapeutic Strategy in Chemo-/Radiotherapy of Head and Neck Cancer. Oxid Med Cell Longev. 2020; 2020:5047987. https://doi.org/10.1155/2020/5047987 [PubMed]
  • 10. Jing YZ, Li SJ, Sun ZJ. Gas and gas-generating nanoplatforms in cancer therapy. J Mater Chem B. 2021; 9:8541–57. https://doi.org/10.1039/d1tb01661j [PubMed]
  • 11. Maharjan PS, Bhattarai HK. Singlet Oxygen, Photodynamic Therapy, and Mechanisms of Cancer Cell Death. J Oncol. 2022; 2022:7211485. https://doi.org/10.1155/2022/7211485 [PubMed]
  • 12. Rencelj A, Gvozdenovic N, Cemazar M. MitomiRs: their roles in mitochondria and importance in cancer cell metabolism. Radiol Oncol. 2021; 55:379–92. https://doi.org/10.2478/raon-2021-0042 [PubMed]
  • 13. Lee J, Roh JL. Induction of ferroptosis in head and neck cancer: A novel bridgehead for fighting cancer resilience. Cancer Lett. 2022; 546:215854. https://doi.org/10.1016/j.canlet.2022.215854 [PubMed]
  • 14. Huang H, Li S, Tang Q, Zhu G. Metabolic Reprogramming and Immune Evasion in Nasopharyngeal Carcinoma. Front Immunol. 2021; 12:680955. https://doi.org/10.3389/fimmu.2021.680955 [PubMed]
  • 15. González-Ruiz L, González-Moles MÁ, González-Ruiz I, Ruiz-Ávila I, Ramos-García P. Prognostic and Clinicopathological Significance of CCND1/Cyclin D1 Upregulation in Melanomas: A Systematic Review and Comprehensive Meta-Analysis. Cancers (Basel). 2021; 13:1314. https://doi.org/10.3390/cancers13061314 [PubMed]
  • 16. Strzyz P. Cell signalling: Signalling to cell cycle arrest. Nat Rev Mol Cell Biol. 2016; 17:536. https://doi.org/10.1038/nrm.2016.108 [PubMed]
  • 17. Zhu D, Huang J, Liu N, Li W, Yan L. PSMC2/CCND1 axis promotes development of ovarian cancer through regulating cell growth, apoptosis and migration. Cell Death Dis. 2021; 12:730. https://doi.org/10.1038/s41419-021-03981-5 [PubMed]
  • 18. Farah CS. Molecular landscape of head and neck cancer and implications for therapy. Ann Transl Med. 2021; 9:915. https://doi.org/10.21037/atm-20-6264 [PubMed]
  • 19. Chu HY, Chen Z, Wang L, Zhang ZK, Tan X, Liu S, Zhang BT, Lu A, Yu Y, Zhang G. Dickkopf-1: A Promising Target for Cancer Immunotherapy. Front Immunol. 2021; 12:658097. https://doi.org/10.3389/fimmu.2021.658097 [PubMed]
  • 20. Dahlmann M, Monks A, Harris ED, Kobelt D, Osterland M, Khaireddine F, Herrmann P, Kemmner W, Burock S, Walther W, Shoemaker RH, Stein U. Combination of Wnt/β-Catenin Targets S100A4 and DKK1 Improves Prognosis of Human Colorectal Cancer. Cancers (Basel). 2021; 14:37. https://doi.org/10.3390/cancers14010037 [PubMed]
  • 21. Forget MA, Turcotte S, Beauseigle D, Godin-Ethier J, Pelletier S, Martin J, Tanguay S, Lapointe R. The Wnt pathway regulator DKK1 is preferentially expressed in hormone-resistant breast tumours and in some common cancer types. Br J Cancer. 2007; 96:646–53. https://doi.org/10.1038/sj.bjc.6603579 [PubMed]
  • 22. Wang Y, Cui D, Li D, Li D, Wang H, Wu Y. Clinical Value on Combined Detection of Serum CA724, DKK1, and TK1 in Diagnosis of Gastric Cancer. J Oncol. 2022; 2022:6941748. https://doi.org/10.1155/2022/6941748 [PubMed]
  • 23. Zheng Y, Zhou Z, Wei R, Xiao C, Zhang H, Fan T, Zheng B, Li C, He J. The RNA-binding protein PCBP1 represses lung adenocarcinoma progression by stabilizing DKK1 mRNA and subsequently downregulating β-catenin. J Transl Med. 2022; 20:343. https://doi.org/10.1186/s12967-022-03552-y [PubMed]
  • 24. Ding G, Lu W, Zhang Q, Li K, Zhou H, Wang F, Zhao C, Fan C, Wang J. ZBTB38 suppresses prostate cancer cell proliferation and migration via directly promoting DKK1 expression. Cell Death Dis. 2021; 12:998. https://doi.org/10.1038/s41419-021-04278-3 [PubMed]
  • 25. Zhang H, Yu C, Dai J, Keller JM, Hua A, Sottnik JL, Shelley G, Hall CL, Park SI, Yao Z, Zhang J, McCauley LK, Keller ET. Parathyroid hormone-related protein inhibits DKK1 expression through c-Jun-mediated inhibition of β-catenin activation of the DKK1 promoter in prostate cancer. Oncogene. 2014; 33:2464–77. https://doi.org/10.1038/onc.2013.203 [PubMed]
  • 26. Zhao Y, Wang C, Goel A. Andrographis overcomes 5-fluorouracil-associated chemoresistance through inhibition of DKK1 in colorectal cancer. Carcinogenesis. 2021; 42:814–25. https://doi.org/10.1093/carcin/bgab027 [PubMed]
  • 27. Qian J, Yi Q. DKK1 as a novel target for myeloma immunotherapy. Oncoimmunology. 2012; 1:756–8. https://doi.org/10.4161/onci.19655 [PubMed]
  • 28. Qian J, Xie J, Hong S, Yang J, Zhang L, Han X, Wang M, Zhan F, Shaughnessy JD Jr, Epstein J, Kwak LW, Yi Q. Dickkopf-1 (DKK1) is a widely expressed and potent tumor-associated antigen in multiple myeloma. Blood. 2007; 110:1587–94. https://doi.org/10.1182/blood-2007-03-082529 [PubMed]
  • 29. D'Amico L, Mahajan S, Capietto AH, Yang Z, Zamani A, Ricci B, Bumpass DB, Meyer M, Su X, Wang-Gillam A, Weilbaecher K, Stewart SA, DeNardo DG, Faccio R. Dickkopf-related protein 1 (Dkk1) regulates the accumulation and function of myeloid derived suppressor cells in cancer. J Exp Med. 2016; 213:827–40. https://doi.org/10.1084/jem.20150950 [PubMed]
  • 30. Niu J, Li XM, Wang X, Liang C, Zhang YD, Li HY, Liu FY, Sun H, Xie SQ, Fang D. DKK1 inhibits breast cancer cell migration and invasion through suppression of β-catenin/MMP7 signaling pathway. Cancer Cell Int. 2019; 19:168. https://doi.org/10.1186/s12935-019-0883-1 [PubMed]
  • 31. Zhu F, Dai SN, Xu DL, Hou CQ, Liu TT, Chen QY, Wu JL, Miao Y. EFNB2 facilitates cell proliferation, migration, and invasion in pancreatic ductal adenocarcinoma via the p53/p21 pathway and EMT. Biomed Pharmacother. 2020; 125:109972. https://doi.org/10.1016/j.biopha.2020.109972 [PubMed]
  • 32. Scherer D, Deutelmoser H, Balavarca Y, Toth R, Habermann N, Buck K, Kap EJ, Botma A, Seibold P, Jansen L, Lorenzo Bermejo J, Weigl K, Benner A, et al. Polymorphisms in the Angiogenesis-Related Genes EFNB2, MMP2 and JAG1 Are Associated with Survival of Colorectal Cancer Patients. Int J Mol Sci. 2020; 21:5395. https://doi.org/10.3390/ijms21155395 [PubMed]
  • 33. Tang XX, Zhao H, Robinson ME, Cohen B, Cnaan A, London W, Cohn SL, Cheung NK, Brodeur GM, Evans AE, Ikegaki N. Implications of EPHB6, EFNB2, and EFNB3 expressions in human neuroblastoma. Proc Natl Acad Sci U S A. 2000; 97:10936–41. https://doi.org/10.1073/pnas.190123297 [PubMed]
  • 34. Zhang Y, Zhang R, Ding X, Ai K. EFNB2 acts as the target of miR-557 to facilitate cell proliferation, migration and invasion in pancreatic ductal adenocarcinoma by bioinformatics analysis and verification. Am J Transl Res. 2018; 10:3514–28. [PubMed]
  • 35. Piffko A, Uhl C, Vajkoczy P, Czabanka M, Broggini T. EphrinB2-EphB4 Signaling in Neurooncological Disease. Int J Mol Sci. 2022; 23:1679. https://doi.org/10.3390/ijms23031679 [PubMed]
  • 36. Shi Y, Wu M, Liu Y, Hu L, Wu H, Xie L, Liu Z, Wu A, Chen L, Xu C. ITGA5 Predicts Dual-Drug Resistance to Temozolomide and Bevacizumab in Glioma. Front Oncol. 2021; 11:769592. https://doi.org/10.3389/fonc.2021.769592 [PubMed]
  • 37. Zhang X, Chen F, Huang P, Wang X, Zhou K, Zhou C, Yu L, Peng Y, Fan J, Zhou J, Lu Z, Hu J, Wang Z. Exosome-depleted MiR-148a-3p derived from Hepatic Stellate Cells Promotes Tumor Progression via ITGA5/PI3K/Akt Axis in Hepatocellular Carcinoma. Int J Biol Sci. 2022; 18:2249–60. https://doi.org/10.7150/ijbs.66184 [PubMed]
  • 38. Chen J, Wang H, Wang J, Niu W, Deng C, Zhou M. LncRNA NEAT1 Enhances Glioma Progression via Regulating the miR-128-3p/ITGA5 Axis. Mol Neurobiol. 2021; 58:5163–77. https://doi.org/10.1007/s12035-021-02474-y [PubMed]
  • 39. Li T, Wu Q, Liu D, Wang X. miR-27b Suppresses Tongue Squamous Cell Carcinoma Epithelial-Mesenchymal Transition by Targeting ITGA5. Onco Targets Ther. 2020; 13:11855–67. https://doi.org/10.2147/OTT.S281211 [PubMed]
  • 40. Zheng W, Jiang C, Li R. Integrin and gene network analysis reveals that ITGA5 and ITGB1 are prognostic in non-small-cell lung cancer. Onco Targets Ther. 2016; 9:2317–27. https://doi.org/10.2147/OTT.S91796 [PubMed]
  • 41. Walter RB, Laszlo GS, Alonzo TA, Gerbing RB, Levy S, Fitzgibbon MP, Gudgeon CJ, Ries RE, Harrington KH, Raimondi SC, Hirsch BA, Gamis AS, W McIntosh M, Meshinchi S. Significance of expression of ITGA5 and its splice variants in acute myeloid leukemia: a report from the Children's Oncology Group. Am J Hematol. 2013; 88:694–702. https://doi.org/10.1002/ajh.23486 [PubMed]
  • 42. Kuninty PR, Bansal R, De Geus SWL, Mardhian DF, Schnittert J, van Baarlen J, Storm G, Bijlsma MF, van Laarhoven HW, Metselaar JM, Kuppen PJK, Vahrmeijer AL, Östman A, et al. ITGA5 inhibition in pancreatic stellate cells attenuates desmoplasia and potentiates efficacy of chemotherapy in pancreatic cancer. Sci Adv. 2019; 5:eaax2770. https://doi.org/10.1126/sciadv.aax2770 [PubMed]
  • 43. Berasain C, Avila MA. Amphiregulin. Semin Cell Dev Biol. 2014; 28:31–41. https://doi.org/10.1016/j.semcdb.2014.01.005 [PubMed]
  • 44. Taniguchi H, Takeuchi S, Fukuda K, Nakagawa T, Arai S, Nanjo S, Yamada T, Yamaguchi H, Mukae H, Yano S. Amphiregulin triggered epidermal growth factor receptor activation confers in vivo crizotinib-resistance of EML4-ALK lung cancer and circumvention by epidermal growth factor receptor inhibitors. Cancer Sci. 2017; 108:53–60. https://doi.org/10.1111/cas.13111 [PubMed]
  • 45. Fang L, Yu Y, Zhang R, He J, Sun YP. Amphiregulin mediates hCG-induced StAR expression and progesterone production in human granulosa cells. Sci Rep. 2016; 6:24917. https://doi.org/10.1038/srep24917 [PubMed]
  • 46. Koeppen M, Lee JW, Seo SW, Brodsky KS, Kreth S, Yang IV, Buttrick PM, Eckle T, Eltzschig HK. Hypoxia-inducible factor 2-alpha-dependent induction of amphiregulin dampens myocardial ischemia-reperfusion injury. Nat Commun. 2018; 9:816. https://doi.org/10.1038/s41467-018-03105-2 [PubMed]
  • 47. Singhal R, Mitta SR, Das NK, Kerk SA, Sajjakulnukit P, Solanki S, Andren A, Kumar R, Olive KP, Banerjee R, Lyssiotis CA, Shah YM. HIF-2α activation potentiates oxidative cell death in colorectal cancers by increasing cellular iron. J Clin Invest. 2021; 131:143691. https://doi.org/10.1172/JCI143691 [PubMed]
  • 48. Sakai R, Ito M, Komai K, Iizuka-Koga M, Matsuo K, Nakayama T, Yoshie O, Amano K, Nishimasu H, Nureki O, Kubo M, Yoshimura A. Kidney GATA3+ regulatory T cells play roles in the convalescence stage after antibody-mediated renal injury. Cell Mol Immunol. 2021; 18:1249–61. https://doi.org/10.1038/s41423-020-00547-x [PubMed]
  • 49. Roy S, Rizvi ZA, Clarke AJ, Macdonald F, Pandey A, Zaiss DMW, Simon AK, Awasthi A. EGFR-HIF1α signaling positively regulates the differentiation of IL-9 producing T helper cells. Nat Commun. 2021; 12:3182. https://doi.org/10.1038/s41467-021-23042-x [PubMed]
  • 50. Xu Q, Long Q, Zhu D, Fu D, Zhang B, Han L, Qian M, Guo J, Xu J, Cao L, Chin YE, Coppé JP, Lam EW, et al. Targeting amphiregulin (AREG) derived from senescent stromal cells diminishes cancer resistance and averts programmed cell death 1 ligand (PD-L1)-mediated immunosuppression. Aging Cell. 2019; 18:e13027. https://doi.org/10.1111/acel.13027 [PubMed]
  • 51. Job S, Reyniès A, Heller B, Weiss A, Guérin E, Macabre C, Ledrappier S, Bour C, Wasylyk C, Etienne-Selloum N, Brino L, Gaiddon C, Wasylyk B, Jung AC. Preferential Response of Basal-Like Head and Neck Squamous Cell Carcinoma Cell Lines to EGFR-Targeted Therapy Depending on EREG-Driven Oncogenic Addiction. Cancers (Basel). 2019; 11:795. https://doi.org/10.3390/cancers11060795 [PubMed]
  • 52. Han SX, Bai E, Jin GH, He CC, Guo XJ, Wang LJ, Li M, Ying X, Zhu Q. Expression and clinical significance of YAP, TAZ, and AREG in hepatocellular carcinoma. J Immunol Res. 2014; 2014:261365. https://doi.org/10.1155/2014/261365 [PubMed]
  • 53. Chang YC, Li CH, Chan MH, Chen MH, Yeh CN, Hsiao M. Regorafenib inhibits epithelial-mesenchymal transition and suppresses cholangiocarcinoma metastasis via YAP1-AREG axis. Cell Death Dis. 2022; 13:391. https://doi.org/10.1038/s41419-022-04816-7 [PubMed]
  • 54. Wang L, Wang L, Zhang H, Lu J, Zhang Z, Wu H, Liang Z. AREG mediates the epithelial-mesenchymal transition in pancreatic cancer cells via the EGFR/ERK/NF-κB signalling pathway. Oncol Rep. 2020; 43:1558–68. https://doi.org/10.3892/or.2020.7523 [PubMed]
  • 55. Elangovan IM, Vaz M, Tamatam CR, Potteti HR, Reddy NM, Reddy SP. FOSL1 Promotes Kras-induced Lung Cancer through Amphiregulin and Cell Survival Gene Regulation. Am J Respir Cell Mol Biol. 2018; 58:625–35. https://doi.org/10.1165/rcmb.2017-0164OC [PubMed]
  • 56. Casimiro MC, Wang C, Li Z, Di Sante G, Willmart NE, Addya S, Chen L, Liu Y, Lisanti MP, Pestell RG. Cyclin D1 determines estrogen signaling in the mammary gland in vivo. Mol Endocrinol. 2013; 27:1415–28. https://doi.org/10.1210/me.2013-1065 [PubMed]
  • 57. Willmarth NE, Ethier SP. Autocrine and juxtacrine effects of amphiregulin on the proliferative, invasive, and migratory properties of normal and neoplastic human mammary epithelial cells. J Biol Chem. 2006; 281:37728–37. https://doi.org/10.1074/jbc.M606532200 [PubMed]
  • 58. Baillo A, Giroux C, Ethier SP. Knock-down of amphiregulin inhibits cellular invasion in inflammatory breast cancer. J Cell Physiol. 2011; 226:2691–701. https://doi.org/10.1002/jcp.22620 [PubMed]
  • 59. Sato T, Tran TH, Peck AR, Liu C, Ertel A, Lin J, Neilson LM, Rui H. Global profiling of prolactin-modulated transcripts in breast cancer in vivo. Mol Cancer. 2013; 12:59. https://doi.org/10.1186/1476-4598-12-59 [PubMed]
  • 60. Naffar-Abu Amara S, Kuiken HJ, Selfors LM, Butler T, Leung ML, Leung CT, Kuhn EP, Kolarova T, Hage C, Ganesh K, Panayiotou R, Foster R, Rueda BR, et al. Transient commensal clonal interactions can drive tumor metastasis. Nat Commun. 2020; 11:5799. https://doi.org/10.1038/s41467-020-19584-1 [PubMed]
  • 61. Goupille C, Vibet S, Frank PG, Mahéo K. EPA and DHA Fatty Acids Induce a Remodeling of Tumor Vasculature and Potentiate Docetaxel Activity. Int J Mol Sci. 2020; 21:4965. https://doi.org/10.3390/ijms21144965 [PubMed]
  • 62. Jiang YQ, Wang ZX, Zhong M, Shen LJ, Han X, Zou X, Liu XY, Deng YN, Yang Y, Chen GH, Deng W, Huang JH. Investigating Mechanisms of Response or Resistance to Immune Checkpoint Inhibitors by Analyzing Cell-Cell Communications in Tumors Before and After Programmed Cell Death-1 (PD-1) Targeted Therapy: An Integrative Analysis Using Single-cell RNA and Bulk-RNA Sequencing Data. Oncoimmunology. 2021; 10:1908010. https://doi.org/10.1080/2162402X.2021.1908010 [PubMed]
  • 63. Ding S, Hong Q, Duan T, Xu Z, He Q, Qiu D, Li L, Yan J, Zhang Q, Mu Z. miRNA-Mediated Low Expression of EPHX3 Is Associated with Poor Prognosis and Tumor Immune Infiltration in Head and Neck Squamous Cell Carcinomas. J Oncol. 2022; 2022:7633720. https://doi.org/10.1155/2022/7633720 [PubMed]
  • 64. Bai G, Song J, Yuan Y, Chen Z, Tian Y, Yin X, Niu Y, Liu J. Systematic analysis of differentially methylated expressed genes and site-specific methylation as potential prognostic markers in head and neck cancer. J Cell Physiol. 2019; 234:22687–702. https://doi.org/10.1002/jcp.28835 [PubMed]
  • 65. Zheng HT, Zhuang ZX, Chen CJ, Liao HY, Chen HL, Hsueh HC, Chen CF, Chen SE, Huang SY. Effects of acute heat stress on protein expression and histone modification in the adrenal gland of male layer-type country chickens. Sci Rep. 2021; 11:6499. https://doi.org/10.1038/s41598-021-85868-1 [PubMed]
  • 66. Ti W, Wei T, Wang J, Cheng Y. Comparative Analysis of Mutation Status and Immune Landscape for Squamous Cell Carcinomas at Different Anatomical sites. Front Immunol. 2022; 13:947712. https://doi.org/10.3389/fimmu.2022.947712 [PubMed]
  • 67. Zhou Z, Xie X, Wang X, Zhang X, Li W, Sun T, Cai Y, Wu J, Dang C, Zhang H. Correlations Between Tumor Mutation Burden and Immunocyte Infiltration and Their Prognostic Value in Colon Cancer. Front Genet. 2021; 12:623424. https://doi.org/10.3389/fgene.2021.623424 [PubMed]
  • 68. Stenman A, Svahn F, Hojjat-Farsangi M, Zedenius J, Söderkvist P, Gimm O, Larsson C, Juhlin CC. Molecular Profiling of Pheochromocytoma and Abdominal Paraganglioma Stratified by the PASS Algorithm Reveals Chromogranin B as Associated With Histologic Prediction of Malignant Behavior. Am J Surg Pathol. 2019; 43:409–21. https://doi.org/10.1097/PAS.0000000000001190 [PubMed]
  • 69. Zhang Y, Chen P, Zhou Q, Wang H, Hua Q, Wang J, Zhong H. A Novel Immune-Related Prognostic Signature in Head and Neck Squamous Cell Carcinoma. Front Genet. 2021; 12:570336. https://doi.org/10.3389/fgene.2021.570336 [PubMed]
  • 70. Eriksson J, Le Joncour V, Jahkola T, Juteau S, Laakkonen P, Saksela O, Hölttä E. Prolyl 4-hydroxylase subunit alpha 1 (P4HA1) is a biomarker of poor prognosis in primary melanomas, and its depletion inhibits melanoma cell invasion and disrupts tumor blood vessel walls. Mol Oncol. 2020; 14:742–62. https://doi.org/10.1002/1878-0261.12649 [PubMed]
  • 71. Wu ZH, Zhong Y, Zhou T, Xiao HJ. miRNA biomarkers for predicting overall survival outcomes for head and neck squamous cell carcinoma. Genomics. 2021; 113:135–41. https://doi.org/10.1016/j.ygeno.2020.12.002 [PubMed]
  • 72. Xiong G, Stewart RL, Chen J, Gao T, Scott TL, Samayoa LM, O'Connor K, Lane AN, Xu R. Collagen prolyl 4-hydroxylase 1 is essential for HIF-1α stabilization and TNBC chemoresistance. Nat Commun. 2018; 9:4456. https://doi.org/10.1038/s41467-018-06893-9 [PubMed]
  • 73. Huang W, Liu X, Zhang Y, Deng M, Li G, Chen G, Yu L, Jin L, Liu T, Wang Y, Chen Y. USP5 promotes breast cancer cell proliferation and metastasis by stabilizing HIF2α. J Cell Physiol. 2022; 237:2211–9. https://doi.org/10.1002/jcp.30686 [PubMed]
  • 74. Afik R, Zigmond E, Vugman M, Klepfish M, Shimshoni E, Pasmanik-Chor M, Shenoy A, Bassat E, Halpern Z, Geiger T, Sagi I, Varol C. Tumor macrophages are pivotal constructors of tumor collagenous matrix. J Exp Med. 2016; 213:2315–31. https://doi.org/10.1084/jem.20151193 [PubMed]
  • 75. Jiang W, Zhou X, Li Z, Liu K, Wang W, Tan R, Cong X, Shan J, Zhan Y, Cui Z, Jiang L, Li Q, Shen S, et al. Prolyl 4-hydroxylase 2 promotes B-cell lymphoma progression via hydroxylation of Carabin. Blood. 2018; 131:1325–36. https://doi.org/10.1182/blood-2017-07-794875 [PubMed]
  • 76. Kang YS, Calvanico NJ, Tomasi TB Jr. Human J-chain: isolation and molecular weight studies. J Immunol. 1974; 112:162–7. [PubMed]
  • 77. Xu AQ, Barbosa RR, Calado DP. Genetic timestamping of plasma cells in vivo reveals tissue-specific homeostatic population turnover. Elife. 2020; 9:e59850. https://doi.org/10.7554/eLife.59850 [PubMed]
  • 78. Zou J, Li Y, Liao N, Liu J, Zhang Q, Luo M, Xiao J, Chen Y, Wang M, Chen K, Zeng J, Mo Z. Identification of key genes associated with polycystic ovary syndrome (PCOS) and ovarian cancer using an integrated bioinformatics analysis. J Ovarian Res. 2022; 15:30. https://doi.org/10.1186/s13048-022-00962-w [PubMed]
  • 79. Pan S, Gao Q, Chen Q, Liu P, Tan Y, Liu F, Xu H. Integrative analysis-based identification and validation of a prognostic immune cell infiltration-based model for patients with advanced gastric cancer. Int Immunopharmacol. 2021; 101:108258. https://doi.org/10.1016/j.intimp.2021.108258 [PubMed]
  • 80. Wang Y, Zhu M, Guo F, Song Y, Fan X, Qin G. Identification of Tumor Microenvironment-Related Prognostic Biomarkers in Luminal Breast Cancer. Front Genet. 2020; 11:555865. https://doi.org/10.3389/fgene.2020.555865 [PubMed]
  • 81. Liu G, Zeng X, Wu B, Zhao J, Pan Y. RNA-Seq analysis of peripheral blood mononuclear cells reveals unique transcriptional signatures associated with radiotherapy response of nasopharyngeal carcinoma and prognosis of head and neck cancer. Cancer Biol Ther. 2020; 21:139–46. https://doi.org/10.1080/15384047.2019.1670521 [PubMed]
  • 82. Zhang T, Nie Y, Gu J, Cai K, Chen X, Li H, Wang J. Identification of Mitochondrial-Related Prognostic Biomarkers Associated With Primary Bile Acid Biosynthesis and Tumor Microenvironment of Hepatocellular Carcinoma. Front Oncol. 2021; 11:587479. https://doi.org/10.3389/fonc.2021.587479 [PubMed]
  • 83. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. 2010; 26:1572–3. https://doi.org/10.1093/bioinformatics/btq170 [PubMed]
  • 84. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018; 24:1550–8. https://doi.org/10.1038/s41591-018-0136-1 [PubMed]
Download PDF
Corresponding Authors

Gengming Cai
cgmkgx@hotmail.com
https://orcid.org/0000-0002-5144-1637

Hui Ye
entyh@126.com
https://orcid.org/0009-0002-8108-6047

Keywords
mitochondrial-related genes head and neck squamous cell carcinoma immunomodulation risk score prognostic model predicting survival
Paper Sections