Integrated transcriptomics identifies HIF1A and GSTP1 as biomarkers for cutaneous squamous cell carcinoma

Introduction

Cutaneous squamous cell carcinoma (cSCC) is a prevalent keratinocyte-derived malignancy and represents a major proportion of global skin cancer cases (1). Its incidence is increasing worldwide, imposing substantial morbidity and economic burden on healthcare systems, particularly among fair-skinned individuals and immunocompromised populations such as organ transplant recipients (2). Notably, cSCC accounts for nearly 20% of skin cancer-related deaths, with an annual mortality comparable to or exceeding that of melanoma (3). The pathogenesis of cSCC is strongly associated with environmental factors, particularly chronic exposure to ultraviolet radiation, which drives the progression from premalignant lesions such as actinic keratosis (4). Clinically, high-risk cSCC typically exhibits differentiation in early stages but may undergo epithelial-mesenchymal transition (EMT), leading to relapse, metastasis, and poor prognosis (5). Although primary cSCC is generally manageable with surgical excision or localized treatments, a subset of cases evolves aggressively, leading to local recurrence, nodal or distant metastasis, and disease-specific mortality (6). Current systemic therapies, including immune checkpoint inhibitors, chemotherapy, and combination regimens, often yield suboptimal outcomes in advanced or metastatic cSCC, contributing to elevated treatment costs and unfavorable prognoses (7). In light of these challenges, the identification of novel biomarkers via integrated transcriptomics holds significant promise for elucidating pathogenic pathways, refining risk stratification, and ultimately improving therapeutic efficacy in cSCC.

Post-translational modifications (PTMs) represent a crucial layer of protein regulation, involving the covalent modification of specific amino acid residues after translation. These alterations profoundly influence protein function, stability, and interactions, thereby modulating diverse cellular processes, including signal transduction to carcinogenesis (8-11). In squamous cell carcinoma (SCC), accumulating experimental evidence underscores the role of specific PTMs in driving tumor progression. For example, genomic analyses have revealed frequent mutations in the lysine-specific histone methyltransferase KMT2D in metastatic cSCC, underscoring the role of histone modifications in aggressive disease phenotypes (12,13). Emerging PTMs such as histone lysine lactylation (Kla) have also been implicated in tumor immune escape. In head and neck squamous cell carcinoma (HNSCC), cleavage under targets and tagmentation (CUT&Tag) analysis identified H3K9la as a modifier of interleukin expression [e.g., interleukin (IL)-6 and IL-10], leading to poor immunotherapy response (14). Consistently, lactylation profiling in oral SCC using immunohistochemical and proteomic analyses has linked Kla to higher malignancy grade, increased tumor proliferation, and microenvironment remodeling (15,16). Moreover, METTL1-mediated m7G RNA modification has been implicated in the development of multiple cancers, including cSCC (17). Complement system components such as C1r may also promote cSCC cell invasion via PTM pathways, whilst C1r knockout inhibits its invasiveness (18). In summary, PTMs participate in cSCC development through multidimensional mechanisms, including metabolic reprogramming, epigenetic regulation, and modulation of the immune microenvironment.

In this study, we bridge this knowledge gap through an integrated transcriptomic analysis of cSCC, focusing on PTM-related processes. Our multi-step approach synergized single-cell and bulk transcriptome data to systematically pinpoint a concise gene signature robustly associated with disease diagnosis and progression. Beyond establishing a validated biomarker model, we delineated the functional implications of these genes in tumor immunity and identified potential targeting agents. This work provides novel, reproducible biomarkers and insights into cSCC pathogenesis, offering a valuable resource for improved diagnosis and the future development of targeted therapies. The flowchart for this study is shown in Figure 1. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0281/rc).

Figure 1 The flowchart of this study. AUC, area under the curve; DEG, differentially expressed gene; GSEA, gene set enrichment analysis; LASSO, least absolute shrinkage and selection operator; PCR, polymerase chain reaction; RF, random forest; ROC, receiver operating characteristic; sc-DEG, single-cell DEG; SVM-RFE, support vector machine-recursive feature elimination; WGCNA, weighted gene co‑expression network analysis; XGBoost, extreme gradient boosting.

Methods

Data sources

Transcriptomic and single-cell RNA sequencing (scRNA-seq) data of cSCC were obtained from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov). The GSE42677 dataset comprised 10 cSCC samples and 10 normal skin samples, generated using the GPL571 platform (HG-U133A_2, Affymetrix Human Genome U133A 2.0 Array). The GSE45164 dataset included 10 cSCC samples and 3 normal samples, annotated on the same platform (GPL571). These two datasets (GSE42677 and GSE45164) were merged and batch effects were corrected as described in reference (19), resulting in a combined training cohort consisting of 20 cSCC samples and 13 normal controls. The validation set was derived from the GSE191334 dataset, which contained 8 normal skin samples and 8 cSCC samples for expression verification of key genes. GSE149614 comprises single-cell transcriptomic data, encompassing 10 control samples and 10 cSCC samples. In addition, a collection of 822 post-translational modification-related genes (PTMRGs) was retrieved from the literature (20), according to the definition and selection criteria specified in that study. These genes are involved in the writing, erasing, or recognition of the 20 PTM types examined in our study.

scRNA-seq analysis

Single-cell transcriptomic analysis was performed using the Seurat package (v4.3.0.1). The Seurat object was constructed with the CreateSeuratObject function under default parameters. Low-quality cells and low-expression genes were filtered by setting thresholds of 200–7,500 for the number of features per cell and 200–50,000 for counts per cell. Additionally, mitochondrial gene content was calculated using the PercentageFeatureSet function, and cells with >10% mitochondrial gene proportion were excluded to ensure data quality. Batch correction among samples was conducted using the Harmony package (v1.2.0) with the RunHarmony function to eliminate potential batch effects. Data normalization was performed with NormalizeData, and the top 2,000 highly variable genes were selected using the “vst” method for downstream analysis. Subsequently, a linear transformation was applied through ScaleData, and principal component analysis (PCA) was conducted for dimensionality reduction. Unsupervised clustering was performed using FindNeighbors and FindClusters functions (resolution =0.8), followed by visualization via uniform manifold approximation and projection (UMAP). Marker genes of each cell cluster were identified using FindAllMarkers, and cluster identities were annotated by comparing these markers with previously reported cell-type-specific markers from reference literature (21). UMAP plots were used to illustrate the spatial distribution of cell clusters across groups, while dot plots generated by Seurat visualized marker gene expression patterns. The relative proportions of different cell types in each sample were visualized as stacked bar charts using the ggpubr package (v0.6.0). Directional scores for each cell were calculated using the AddModuleScore function, and Wilcoxon rank-sum tests were applied to identify cell types showing significant differences in directional scores or abundance between groups (P<0.05). These significantly altered cell populations were defined as key cells. Finally, differential expression analysis was performed between the control and disease groups for key cells (|log₂fold change (FC)| >0.25 and P<0.05) to identify single-cell-differentially expressed genes (sc-DEGs).

Differential expression analysis and weighted gene co-expression network analysis

Based on the training set, differential expression analysis was first performed using the limma package. A threshold of |log₂FC| >0.5 and P<0.05 was set to identify DEGs. To identify co-expression gene modules and their core genes associated with disease phenotypes, weighted gene co‑expression network analysis (WGCNA) was employed. The pickSoftThreshold function evaluated the scale-free topology fit index and average connectivity under different soft threshold powers β to determine the optimal soft threshold. The selected β was used to construct an adjacency matrix, which was then transformed into a topology overlap matrix (TOM) to measure co-expression similarity between genes. Hierarchical clustering was performed based on TOM’s dissimilarity matrix, and gene modules were partitioned using the dynamic tree cutting algorithm [minimum module size =30, module merging threshold (cut height) set to 0.25]. The module eigengene (ME) is calculated for each module, and Pearson correlation analysis is performed to assess the association between the ME and the clinical phenotype (disease/control). Modules significantly associated with disease (|cor| >0.5 and P<0.05) are defined as key modules. Genes within key modules were defined as key module genes.

Functional enrichment analysis of candidate genes

An intersection analysis was performed on sc-DEGs, PTMRGs, key module genes, and DEGs to obtain a candidate gene. To further explore the potential biological functions of these candidate genes in disease progression, functional enrichment analysis was conducted using the clusterProfiler package (22). This specifically included annotation analysis based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO). GO annotation encompasses three major categories: biological processes (BPs), cellular components (CCs), and molecular functions (MFs). By calculating the enrichment levels of candidate genes across KEGG pathways and GO entries, and employing Fisher’s exact test for statistical significance assessment, pathways and entries showing significant enrichment (P<0.05) were identified. Plot chromosome circles using the RCircos package to visualise gene locations on chromosomes.

Identification of key genes

To further identify key genes with potential diagnostic value in disease classification, we applied five machine learning algorithms for feature selection to the candidate genes identified in the preliminary screening. These included least absolute shrinkage and selection operator (LASSO), support vector machine-recursive feature elimination (SVM-RFE), Boruta, extreme gradient boosting (XGBoost), and random forest (RF). Among these, LASSO logistic regression analysis was implemented using the glmnet package (23), effectively identifying important features and collapsing the coefficients of insignificant variables to zero. The Boruta algorithm, based on the RF model, performed automatic feature selection via the Boruta package; the RF method utilised the randomForest package to calculate feature importance and screen significant variables (P<0.05); the XGBoost algorithm, implemented through the XGBoost package, screened candidate genes based on feature importance values. Ultimately, the intersection genes jointly identified by all five algorithms were defined as key candidate genes closely associated with disease diagnosis and progression. To further refine the key genes, the Wilcoxon signed-rank test was first employed to compare expression level differences between the key candidate genes in the training and validation datasets. Subsequently, the pROC package was utilised for receiver operating characteristic (ROC) curve analysis to evaluate the diagnostic accuracy of each key gene in both the training and validation sets. Predictive performance was measured by calculating the area under the curve (AUC), where an AUC close to 1 indicates high diagnostic efficacy. Ultimately, genes exhibiting consistent expression trends across both training and validation sets, with significant differences (P<0.05) and an AUC >0.7, were selected as the final key genes.

Predictive modeling and model interpretation

To comprehensively assess disease risk, a predictive model incorporating 12 machine learning algorithms was constructed, utilising key genes as predictive features. All modeling procedures were performed using the caret package (v 6.0-94). These algorithms encompassed RF, gradient boosting, support vector machine Kernel (SVM_Kernel), logistic model, neighbour method, partial least squares regression (PLS Model), decision tree, neural net, Bayesian method, discriminant model, LASSO regression, and adaptive boosting. To ensure robust parameter estimation and avoid overfitting, we applied 10-fold cross-validation repeated three times on the training set for each algorithm. Hyperparameters were tuned using default settings as implemented in caret. Model discriminatory performance was evaluated by the AUC using the pROC package (v 1.18.5). The model with the highest AUC value was selected as the final predictive model. Considering the prevalent ‘black box’ issue in machine learning models, the SHapley Additive Explanations (SHAP) method was further applied to conduct interpretative analysis on the optimal model. SHAP values quantify each feature’s contribution to the model output in a game‑theoretic manner, where positive values indicate that higher expression of the gene increases the predicted probability of disease, and negative values suggest a protective or risk‑reducing effect. SHAP values were calculated using the Permshap function, and visualizations were generated with the shapviz package.

Functional enrichment analysis

To identify biological pathways associated with key genes, the training dataset was integrated with KEGG gene sets from the Molecular Signatures Database (MSigDB). Correlation analysis between key genes and all other genes was performed, followed by single-sample gene set enrichment analysis (ssGSEA) using the clusterProfiler package (version 4.14.4) to determine potential KEGG pathways related to each key gene. Statistical significance was defined as P<0.05. Furthermore, to predict the subcellular localization of key proteins, the FASTA sequences of each gene were retrieved from the UniProt database (www.uniprot.org) and analyzed using the WoLF PSORT tool (wolfpsort.hgc.jp). The location with the highest prediction score was considered as the specific subcellular localization for each key gene, providing insights into their functional distribution within cells.

Immune microenvironment analyses

To characterise differences in immune cell and stromal cell infiltration between cSCC and control samples, ssGSEA was performed using the gene set variation analysis (GSVA) software package (v1.48.2). The Wilcoxon signed-rank test assessed differences in immune cell enrichment scores between the cSCC and control groups, with P<0.05 indicating statistical significance; cell types exhibiting significant differences were defined as differentially enriched cells. Spearman’s correlation test examined associations between key genes and immune cells; correlations with |cor| >0.3 and P<0.05 were deemed clinically significant. To further deepen the immune analysis, we investigated the correlations between the key genes and immune checkpoint molecules (24) using Spearman’s rank correlation test. A correlation coefficient |cor| >0.3 and P<0.05 was considered statistically significant, and the results were visualized using a heatmap.

Construction of regulatory network and molecular docking

To further explore the functional roles of key genes at the protein level, the STRING database (https://string-db.org) was used to predict protein-protein interactions (PPI) associated with the key genes. Subsequently, the GeneMANIA platform (https://genemania.org) was employed to identify functionally similar genes or proteins based on extensive genomic and proteomic datasets, thereby constructing a gene-gene interaction (GGI) network. In terms of drug discovery, potential compounds targeting key genes were predicted using the Drug-Gene Interaction Database (DGIdb). Drugs that were commonly associated with multiple key genes were selected for subsequent molecular docking analyses. The three-dimensional (3D) structures of the proteins encoded by the key genes were obtained from UniProt (https://www.uniprot.org/), while the corresponding 3D structures of candidate drugs were downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/). Protein structures were preprocessed using AutoDock software by removing water molecules and unrelated ligands, adding hydrogen atoms, and saving the processed files in Protein Data Bank, partial charges and atom types (PDBQT) format to obtain optimized protein conformations suitable for docking. Drug ligands were also hydrogenated using AutoDock and defined as docking ligands. Molecular docking between each key gene–encoded protein and its corresponding drug was performed using AutoDock Vina. The complex with the lowest binding energy was considered to represent the most stable interaction. Finally, PyMOL and LigPlot software were used to visualize the optimal docking results, providing an intuitive depiction of the binding patterns and interaction interfaces between proteins and ligands.

Cell-cell communication analysis and Pseudotime trajectory analysis

Cell-cell communication analysis was performed using the CellChat package (v2.1.2, https://github.com/sqjin/CellChat). The “Secreted Signaling” interaction database was utilized to integrate single-cell transcriptomic data from both normal and disease groups. Comparative analyses were conducted between the two groups to identify alterations in intercellular communication patterns, including differences in signaling pathway interaction strength and the major ligand-receptor pairs between key cell types and other cell populations. To evaluate the relationship between enhanced intercellular signaling and key gene expression, PPI networks were constructed using the STRING database to screen for ligand-receptor genes that interact with HIF1A or GSTP1. Subsequently, Spearman correlation analysis was performed between the expression levels of HIF1A/GSTP1 and these ligand-receptor-related genes within key cell types; |cor| >0.3 and P<0.05 were considered statistically significant.

Furthermore, pseudotime trajectory analysis was conducted using the Monocle package (v2.32.0). DEGs between normal and disease states within key cell types (criteria: |avg_log₂FC| >1 and P<0.05) were used to infer cellular developmental trajectories.

Quantitative real-time polymerase chain reaction (qRT-PCR)

A total of 10 paired cSCC and adjacent non-tumor tissue samples were collected from Shanghai Fourth People’s Hospital. The study was approved by the Ethics Committee of Shanghai Fourth People’s Hospital (No. TJBH07224101), and written informed consent was obtained from all participants. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Total RNA was extracted using the Animal Total RNA Isolation Kit (Foregene; RE-03014) according to the manufacturer’s instructions, and RNA concentration and purity were assessed with a nucleic acid analyzer (K2800, Beijing Kaiao Technology). Quantitative PCR was conducted on a CFX Connect Real-Time PCR System (Bio-Rad) using EvaGreen Express 2× qPCR MasterMix (ABM; MasterMix-ES). All reactions were performed in triplicate. Relative gene expression levels were calculated using the 2^−ΔΔCT method. Primer sequences are provided in Table S1.

Statistical analysis

Statistical computations were executed in R (V 4.3.3) software. Analyses of scRNA-seq data leaned on the Seurat (V 5.3.0) package. Comparisons across groups were conducted via the Wilcoxon test, where statistical significance was set at P<0.05.

Results

Single-cell profiling and cell type identification

Single-cell RNA-seq data were processed using the Seurat pipeline to filter out low-quality cells (Figure S1A). Following normalization, highly variable genes were selected for downstream analysis (Figure S1B). PCA and batch correction were performed, with the elbow plot indicating the top 30 principal components to be retained. No significant batch effects were observed among samples (Figure S1C,S1D). Unsupervised clustering based on these components identified 29 distinct cell clusters (Figure 2A). Cell type annotation was conducted and validated by visualizing canonical marker gene expression via bubble plots. Major cell types identified included B cells, dendritic cells (DCs), endothelial cells, fibroblasts, keratinocytes, Langerhans cells, macrophages, melanocytes, and T cells (Figure 2B,2C). Notably, the distribution of annotated cell types differed significantly between the experimental and control groups (Figure 2D-2F). Integrating PTMRGs, GSVA was applied to calculate directionality scores across cell types, assessing their distribution patterns and differences between groups. Directionality scores were significantly enriched in DCs and other specific populations (P<0.05, Figure 2G). Given that Keratinocytes exhibited significant differences in both abundance and directionality scores between groups, they were designated as key cells for subsequent analyses. Finally, differential expression analysis was performed between the normal and cSCC groups of key cells, yielding 999 sc-DEGs (Figure 2H).

Figure 2 Single-cell profiling and cell-type identification in the study samples. (A) UMAP visualization of 29 distinct cell clusters identified from single-cell RNA-seq data after Seurat-based preprocessing. (B,C) Marker-based annotation of major cell lineages, with representative marker expression shown as bubble plots and feature plots. (D-F) Comparison of annotated cell-type distributions between experimental and control groups. (G) GSVA-derived directionality scores across cell types incorporating PTMRGs; directionality was significantly enriched in DCs and select other populations (P<0.05). (H) Differential expression analysis within the key cell type. *, P<0.05; **, P<0.01; ns, not significant. DC, dendritic cell; GSVA, gene set variation analysis; PTMRG, post-translational modification-related gene; UMAP, uniform manifold approximation and projection.

Identification of DEGs and key module genes

In the training dataset, differential expression analysis identified a total of 5,771 genes that were significantly differentially expressed between the disease and control groups, including 3,047 upregulated and 2,724 downregulated genes (Figure S2A,S2B). The pickSoftThreshold function was employed to test various soft-threshold powers (β), and the optimal value of 9 was selected based on scale-free topology criteria. Using this threshold, an adjacency matrix was constructed and subsequently transformed into a TOM to quantify gene co-expression similarity (Figure S2C,S2D). Hierarchical clustering combined with dynamic tree cutting was performed, a total of 17 distinct co-expression modules were identified (Figure S2E). Analysis of ME connectivity revealed that inter-module distances exceeded 0.25, indicating good independence among modules. Correlation analysis between MEs and disease status showed that six modules (turquoise, brown, pink, blue, green, and black) were significantly associated with the disease group. These modules were defined as key modules, and the 10,551 genes contained within them were designated as key module genes (Figure S2F).

A total of 19 candidate genes were screened

An intersection analysis was performed among sc-DEGs, PTMRGs, bulk DEGs, and key module genes, yielding 19 candidate genes (Figure 3A). Functional enrichment analyses for these candidate genes were conducted using GO and KEGG annotations. The GO enrichment identified 401 enriched terms, comprising 319 BPs, 15 CCs, and 67 MFs. Notably enriched GO terms include the hypoxia-inducible factor-1 alpha signaling pathway, regulation of cytokine-mediated signaling pathway, and regulation of apoptotic signaling pathway (Figure 3B). The KEGG analysis revealed 14 enriched pathways, including ubiquitin-mediated proteolysis and glycolysis/gluconeogenesis (Figure 3C). Chromosomal localization was shown in Figure 3D.

Figure 3 Identification and functional enrichment of candidate genes. (A) Venn plot of candidate genes. (B) GO enrichment overview with 401 terms (BP 319, CC 15, MF 67). (C) KEGG enrichment. (D) Chromosomal localization of the candidate genes. BP, biological process; CC, cellular component; DEG, differentially expressed gene; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; PTMRG, post-translational modification-related gene.

HIF1A and GSTP1 were identified as key genes

Based on the previously identified candidate genes, five machine learning algorithms were applied to further select diagnostic gene targets with significant features. LASSO logistic regression analysis performed identified six genes (Figure 4A). The RF algorithm selected the top five genes (Figure 4B), Boruta identified 13 significant feature genes (Figure 4C), SVM-RFE selected 17 key genes (Figure 4D), and XGBoost identified nine feature genes (Figure 4E). Through the combined screening of gene intersections using five methods, three key candidate genes were ultimately identified: HIF1A, UBC, and GSTP1. These genes were designated as candidate key genes for cSCC (Figure 4F).

Figure 4 Machine learning-based feature selection of candidate genes. (A) Genes selected by LASSO logistic regression analysis. (B) Top features identified by RF algorithm. (C) Feature genes selected by Boruta algorithm. (D) Key genes selected by SVM-RFE algorithm. (E) Features identified by XGBoost algorithm. (F) Intersection of genes selected by all five methods. CV, cross-validation; LASSO, least absolute shrinkage and selection operator; RF, random forest; SVM-RFE, support vector machine-recursive feature elimination; XGBoost, extreme gradient boosting.

To further identify key genes, the expression levels of candidate genes were compared between the disease and control groups in both the training dataset (Figure 5A) and the validation dataset (Figure 5B). The results showed that HIF1A and GSTP1 were significantly upregulated in the disease group compared with the controls (P<0.05), displaying consistent trends across both datasets. To evaluate the diagnostic accuracy of each key gene, ROC curve analysis was performed in both the training and independent validation cohorts. The AUC was calculated for each gene, where an AUC value approaching 1 indicates superior diagnostic performance. In both the training set (Figure 5C) and validation set (Figure 5D), HIF1A and GSTP1 achieved AUC values greater than 0.7 for distinguishing disease samples from controls, demonstrating strong diagnostic potential. qRT-PCR analysis further revealed that HIF1A and GSTP1 were significantly upregulated in the cSCC samples (Figure 5E). Consequently, HIF1A and GSTP1 were identified as the final key genes.

Figure 5 Validation and diagnostic performance of key genes HIF1A and GSTP1. (A,B) Expression levels of candidate key genes in training and validation datasets, respectively. (C,D) ROC curve analyses evaluating the diagnostic accuracy of HIF1A and GSTP1 in training and validation cohorts. (E) Relative mRNA expression levels of HIF1A and GSTP1. Data are presented as mean ± SD (n=10 per group). Statistical significance was assessed using the Student’s t-test (P<0.05). *, P<0.05; **, P<0.01; ****, P<0.0001; ns, not significant. AUC, area under the curve; cSCC, cutaneous squamous cell carcinoma; ROC, receiver operating characteristic; SCC, squamous cell carcinoma; SD, standard deviation.

Diagnostic model construction and SHAP-based interpretability

To construct a diagnostic model, ten machine learning algorithms were applied to the training dataset to evaluate the diagnostic value of key genes. Residual distributions were compared, and ROC curves from cross-validation assessed each model’s performance. The LASSO model achieved the highest AUC, demonstrating superior discriminative ability compared to other models (Figure S3A). SHAP analysis revealed the contributions of key genes HIF1A and GSTP1 in the predictive model. Specifically, SHAP values were computed for each feature to quantify their marginal contributions to the model’s output, with positive values indicating promotion of disease prediction. HIF1A consistently exhibited high SHAP values, indicating its strong influence on prediction accuracy, while high expression levels of both genes promoted disease prediction (Figure S3B,S3C). Interaction analysis showed that high HIF1A expression often followed elevated GSTP1 levels and jointly increased disease prediction probability (Figure S3D). Key predictive features HIF1A and GSTP1 showed strong contributions to model performance. In sample 3, f(x) = 0.549 and E[f(x)] = 0.606, confirming the robustness of the predictive model (Figure S3E). In summary, HIF1A and GSTP1 were identified as key drivers of model performance.

GSEA, subcellular localization and immune infiltration analysis

GSEA enrichment analysis of HIF1A and GSTP1 revealed significant enrichment in several KEGG pathways, including the Neurotrophin signaling pathway, Chemokine signaling pathway, Toll-like receptor signaling pathway, Wnt signaling pathway, and NOD-like receptor signaling pathway (Figure 6A). Furthermore, subcellular localization analysis indicated that HIF1A was predominantly located in the nucleus, whereas GSTP1 was mainly localized in mitochondria (Figure 6B).

Figure 6 GSEA, subcellular localization and immune infiltration analyses. (A) GSEA showing enrichment of HIF1A and GSTP1 in several KEGG pathways. (B) Subcellular localization. (C) Heatmap of immune cell infiltration across samples derived from ssGSEA. (D) Differences in 18 immune cell subsets between disease and control groups. (E) Correlations between HIF1A and GSTP1 expression and infiltration levels of immune cells across subsets. (F) Heatmap of the correlation between key genes and immune checkpoints. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001; ns, P>0.05. GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; SCC, squamous cell carcinoma; ssGSEA, single-sample gene set enrichment analysis.

Enrichment analysis indicated that immune regulation plays a crucial role in disease progression. The ssGSEA method was applied to evaluate immune cell infiltration in control and disease samples, with a heatmap illustrating the types and abundance of infiltrating immune cells in each sample (Figure 6C). Comparative analysis revealed 18 immune cell subsets showing significant differences between groups; notably, neutrophils, plasmacytoid DCs, and macrophages were markedly increased in the disease group (P<0.05, Figure 6D). Furthermore, correlation analysis indicated that the expression levels of HIF1A and GSTP1 were positively correlated with the infiltration of most immune cells, including macrophages, NK cells, and neutrophils, while exhibiting a negative correlation with immature DCs and immature B cells (P<0.05) (Figure 6E). The correlation analysis between immune checkpoint molecules and the two key genes (GSTP1 and HIF1A) revealed that both GSTP1 and HIF1A were significantly negatively correlated with immune activation‑related molecules, including CD28, CD27, and TNFSF4 (P<0.001). Meanwhile, both genes showed significant positive correlations with immune regulatory molecules ICOS and CD80 (P<0.001). Additionally, HIF1A was positively correlated with CD48 and CD40 (P<0.001 and P<0.01, respectively) and negatively correlated with CD40LG (P<0.001). NRP1 exhibited a moderate positive correlation only with GSTP1 (P<0.01) but showed no significant association with HIF1A (Figure 6F).

PPI and GGI network construction, drug prediction and molecular docking

PPI network analysis identified 4 nodes and 3 edges (Figure S4A). Genemania revealed interactions between 2 hub genes and a total of 20 highly correlated genes. GO/KEGG enrichment analyses showed enrichment primarily in pathways related to cellular response to hypoxia and oxygen levels, regulation of glycolytic processes, cellular response to oxidative stress, negative regulation of apoptotic signaling, glycolytic process, and regulation of vascular endothelial growth factor (VEGF) production (Figure S4B). Drug prediction identified three compounds that target the key genes (Figure S4C), including camptothecin, hydroquinone, and resveratrol. Docking results (Table S2) indicate interactions such as HIF1A with camptothecin (−8.7 kcal/mol) and GSTP1 with camptothecin (−8.6 kcal/mol), demonstrating favorable binding affinities (Figure S4D).

Cell-cell communication analysis and keratinocyte pseudotime trajectory

Based on inter-sample differences, cell-cell communication analysis between the two groups was performed using the Secreted Signaling database. The results showed that the disease group exhibited stronger interaction intensity and a greater number of intercellular communications compared with the control group (Figure 7A), primarily involving signaling pathways such as VEGF and endothelin (EDN) (Figure 7B).

Figure 7 Cell-cell communication analysis and keratinocyte pseudotime trajectory. (A) Comparison of intercellular communication strength between disease and control groups. (B) Major signaling pathways involved in altered cell communication. (C) Interactions between keratinocytes as target cells and various immune and stromal cell types. (D) Significant ligand-receptor pairs between keratinocytes as source cells and B cells with increased communication in the disease group. (F) The PPI network between the ligand and GSTP1 and HIF1A. (G) Heatmap of the correlation between the ligand and GSTP1 and HIF1A. (H) Pseudotime trajectory of keratinocyte differentiation from normal to metastatic states, showing seven distinct states. (I,J) Temporal gene clusters enriched in keratinocyte development, proliferation, and differentiation pathways. (K) Expression dynamics of HIF1A and GSTP1 along keratinocyte differentiation trajectories. PPI, protein-protein interaction.

Furthermore, keratinocytes acting as target cells were found to interact with multiple cell types, including B cells, DCs, endothelial cells, fibroblasts, Langerhans cells, macrophages, melanocytes, and T cells (Figure 7C). Notably, when keratinocytes functioned as source cells interacting with B cells, significant ligand-receptor pairs were identified, including MIF-(CD74+CXCR4) and MIF-(CD74+CD44). The communication probabilities of these interactions were markedly higher in the disease group than in the control group (Figure 7D,7E). PPI analysis showed that these enhanced ligand-receptor pairs have known protein interactions with HIF1A, and GSTP1 also interacts with HIF1A (Figure 7F). Furthermore, the expression levels of HIF1A and GSTP1 were significantly positively correlated with nine related genes (IFNG, ITGAV, WNT5A, PLAUR, PLAU, IGFBP3, SDC1, SPP1, MDK) within key cell types (|cor| >0.3, P<0.05, Figure 7G).

Pseudotime trajectory analysis of the key cell type, keratinocytes, revealed a differentiation progression primarily from normal to metastatic states. Keratinocytes were categorized into seven distinct differentiation states, with cells from the normal group predominantly occupying the initial states, while those from the metastatic group were enriched in later stages (Figure 7H). Three gene clusters exhibiting temporal expression patterns were identified along the trajectory, which were significantly enriched in pathways related to keratinocyte biology, including keratinocyte development, regulation of keratinocyte proliferation, keratinocyte differentiation, and regulation of keratinocyte differentiation (Figure 7I,7J). Finally, Pseudotime trajectory analysis revealed that during the transition from normal to malignant states, HIF1A expression progressively increased along the pseudotime axis, and GSTP1 showed a similar upward trend, indicating that both genes are upregulated during malignant transformation (Figure 7K).

Discussion

Leveraging integrative transcriptome and single-cell analyses, this study systematically identifies and characterizes the pivotal roles of HIF1A and GSTP1 in the pathogenesis of cSCC. PTMs contribute to the pathogenesis of cSCC through multiple mechanisms, including epigenetic regulation, aberrant signaling pathway activation, and remodeling of the tumor immune microenvironment (14,25,26). Our single-cell analysis pinpointed keratinocytes as a central compartment undergoing substantial reprogramming. Through a consensus of machine learning algorithms applied to PTM-related gene sets, HIF1A and GSTP1 emerged as consistently upregulated genes with robust diagnostic power (AUC >0.7), underscoring their robust diagnostic potential.

HIF1A, a master regulator of cellular adaptation to hypoxia, is a well-established driver of angiogenesis and metastasis in various cancers, including SCCs (27,28). Its significance in SCC is supported by transcriptomic analyses identifying HIF1A as a core upregulated transcription factor (29). Mechanistically, HIF1A promotes angiogenesis by directly transactivating VEGF, a relationship observed in HNSCCs where their expression is synchronized (30,31). In various SCC subtypes, including HNSCC and esophageal SCC, HIF1A overexpression has been linked to poor prognosis via activation of glycolysis, Wnt/β-catenin signaling, and EMT, thereby enhancing stemness and proliferative capacity (32-35). Additional mechanisms involve interactions with Akt signaling, regulation of lactate dehydrogenase A to support glycolytic metabolism (34), and modulation of chemotherapy sensitivity through mitochondrial fission factor (Mff) (36). In tumor cells, preventing HIF1A ubiquitination to stabilize the protein promotes tumor angiogenesis within the hypoxic microenvironment (37). HIF1A also enhances PGK1 succinylation levels in glioblastoma, facilitating tumor growth (38). Through GSEA and cell-cell communication analysis, the present study confirms HIF1A upregulation in cSCC and reveals concomitant activation of the VEGF signaling pathway. These findings consolidate the role of HIF1A in driving angiogenesis in cSCC and further suggest its involvement in a broader pro-tumorigenic network within the tumor microenvironment. Collectively, this positions HIF1A as both a prognostic marker and a potential therapeutic target in cSCC.

GSTP1, a key glutathione S-transferase, plays a dual role in cancer by managing oxidative stress and promoting malignant progression (39). It facilitates intracellular reactive oxygen species (ROS) clearance and is central to maintaining redox homeostasis (40,41). Paradoxically, its aberrant overexpression is associated with aggressive tumor phenotypes. In lung cancer, high GSTP1 expression correlates with stemness, metastasis, and therapy resistance (42). Similarly, in esophageal SCC, GSTP1 inhibits apoptosis by modulating oxidative stress (43). In HNSCC, GSTP1 overexpression has been documented, and silencing GSTP1 in HNSCC cell lines significantly reduces proliferation and increases apoptosis (44). Genetic polymorphisms in GSTP1 also influence susceptibility to SCC and therapeutic outcomes. For example, the Arg187Trp variant is associated with increased oral SCC risk (45), and the c.341C>T variant confers a higher risk, particularly among tobacco users (46). Moreover, the Ile105Val polymorphism has been linked to greater toxicity and poorer prognosis in HNSCC patients receiving cisplatin-based chemoradiotherapy (47). Regarding PTMs, ubiquitin-mediated proteasome degradation suppresses the pro-tumorigenic function of GSTP1 in various cancers (48,49). In breast cancer, GSTP1 promotes tumorigenesis by regulating glucose-6-phosphate dehydrogenase (G6PD) phosphorylation (50). In glioma, serine phosphorylation of GSTP1 enhances its ability to metabolize cisplatin, thereby inducing resistance to cisplatin chemotherapy (51). These findings underscore GSTP1’s relevance as both a biomarker and a potential therapeutic target. Collectively, these insights suggest that in cSCC, GSTP1 may enhance tumor cell survival and chemoresistance by suppressing apoptosis, thereby facilitating adaptation to microenvironmental stress.

Although neither HIF1A nor GSTP1 is a classical PTM enzyme, both are critical substrates and functional effectors of PTMs under hypoxic and oxidative stress. Their activities are finely tuned by ubiquitination, phosphorylation, and succinylation, which directly impact cSCC progression. Therefore, they are legitimately categorized as “PTM-related genes” from a functional perspective. Based on convergent evidence from transcriptomic, cellular, and network-level analyses, we propose that HIF1A and GSTP1 may constitute a synergistic regulatory axis driving cSCC progression. This hypothesis is supported by their complementary biological functions and coordinated association with pathological features in our data. A potential model involves initial HIF1A activation by tumoral hypoxia, leading to the upregulation of pro-angiogenic and remodeling factors (e.g., VEGF, EDN), a process supported by our GSEA and cell-cell communication analyses. The consequent rapid proliferation and aberrant vasculature likely elevate metabolic and oxidative stress. In this context, mitochondrial GSTP1 could facilitate tumor cell adaptation by scavenging excess ROS and maintaining redox homeostasis, thereby promoting survival. Furthermore, both genes may cooperatively shape an immunosuppressive microenvironment. Their expression correlated positively with infiltrating macrophages and neutrophils, cell types often bolstered by HIF1A-driven metabolic reprogramming (52). Notably, our cell-cell communication analysis revealed enhanced MIF signaling from keratinocytes to B cells in cSCC. Given the established role of MIF in promoting immune escape and immunotherapy resistance in other cancers (53,54), this pathway may represent a shared mechanism through which the HIF1A-GSTP1 axis influences the tumor immune landscape. Finally, PPI and genetic interaction network analyses, which show co-enrichment of these genes in hypoxia response, glycolysis, and oxidative stress pathways, provide further conceptual support for their functional interconnection. This integrative hypothesis, derived from transcriptomic, cellular, and network-level evidence, delineates a plausible mechanistic framework for future experimental validation.

Therapeutically, our drug prediction and molecular docking analyses provide a computational foundation for exploring potential treatment strategies in cSCC. The results indicate that small-molecule compounds, including camptothecin and resveratrol, may form stable bindings with HIF1A or GSTP1 proteins. Camptothecin, a classic topoisomerase I inhibitor, induces apoptosis through DNA strand breaks (55). Notably, prior studies indicate it may influence HIF1α mRNA levels (56), and its toxicity profile could involve GSTP1-mediated oxidative stress pathways (49). Separately, the natural polyphenol resveratrol possesses established anti-inflammatory and antioxidant properties (57). In pancreatic cancer models, molecular docking and dynamics simulations have shown that resveratrol can directly bind to and inhibit HIF1α (58). In colorectal cancer, resveratrol influences chemotherapy responses through the β1-integrin/HIF-1α axis (59), while also enhancing cellular antioxidant capacity via activation of the Nrf2 pathway (60), which may indirectly regulate GSTP1 expression or activity. Integrating these pharmacological insights with our proposed HIF1A-GSTP1 axis model, we speculate that dual targeting—disrupting hypoxia adaptation via HIF1A inhibition while impairing oxidative damage repair and detoxification via GSTP1—could synergistically sensitize cSCC cells to chemotherapy and oxidative stress, thereby offering a promising combinatorial therapeutic approach.

Beyond the central roles of HIF1A and GSTP1, our study reveals additional features of the cSCC tumor microenvironment that merit attention. We observed significant alterations in the PTM-related gene expression profile in DCs, suggesting a potential impairment of their antigen-presenting function (61). Given the decisive role of DCs in initiating anti-tumor immunity, such dysfunction may contribute to immune escape in cSCC, which could be particularly relevant for understanding the high incidence of this malignancy in immunocompromised individuals (62,63). Furthermore, our GSEA and singlecell communication analyses identified coordinated activation of the VEGF and EDN signaling pathways. VEGF acts not only as a key regulator of tumor angiogenesis (64,65) but also participates in establishing a hypoxiarelated immunosuppressive microenvironment, forming a positive feedback loop that exacerbates vascular abnormality and immune evasion (66). The EDN pathway, which regulates EMT and cellular responses to hypoxia, can further upregulate VEGF expression, thereby cooperatively promoting angiogenesis and EMT (67,68). These findings suggest that in cSCC, alongside the canonical HIF1A-driven axis, a microenvironment remodeling program coordinated by VEGF and EDN signaling may actively support tumor progression and metastasis.

Despite the robust bioinformatic evidence, this study has limitations. First, although the sample size was adequate for exploratory discovery, it may not provide sufficient statistical power to detect associations with specific clinical subtypes. Second, the lack of prognostic data precludes evaluating the clinical value of HIF1A and GSTP1 in predicting patient outcomes. Finally, the functional roles of these genes, inferred computationally, require experimental validation in vitro and in vivo to confirm their necessity and explore their potential synergy in cSCC pathogenesis.

Conclusions

In conclusion, this study defines HIF1A and GSTP1 as pivotal drivers and diagnostic biomarkers in cSCC through a compelling multi-omics narrative that spans from cellular reprogramming to druggable targets. Our findings provide a substantive resource for understanding PTM-mediated oncogenesis and a concrete foundation for developing diagnostic and therapeutic strategies. Future studies focusing on the experimental dissection of the HIF1A-GSTP1 axis and its clinical prognostic value are warranted to advance these findings toward precision oncology applications.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0281/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0281/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0281/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0281/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Shanghai Fourth People’s Hospital (No. TJBH07224101) and informed consent was obtained from all participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Nassir S, Yousif M, Li X, et al. Multiomic Sequencing of Intermediate- to High-Risk Cutaneous Squamous Cell Carcinoma Identifies Critical Genes and Expression Patterns Associated with Disease and Poor Outcomes. J Invest Dermatol 2025;145:2060-2070.e5. [Crossref] [PubMed]
2. Klein JC, Shahwan KT, Petric UB, et al. Impact of immunosuppression on cutaneous squamous cell carcinoma outcomes. J Am Acad Dermatol 2026;94:151-60. [Crossref] [PubMed]
3. Cheraghlou S, Stevenson ML, Christensen SR, et al. Standardizing Retrospective Observational Research in Cutaneous Squamous Cell Carcinoma: Expert Panel Guidelines from ITSCC. JAMA Dermatol 2024;160:989-92. [Crossref] [PubMed]
4. Wang Y, Xu X, Jiang G. Microplastics exposure promotes the proliferation of skin cancer cells but inhibits the growth of normal skin cells by regulating the inflammatory process. Ecotoxicol Environ Saf 2023;267:115636. [Crossref] [PubMed]
5. Bone M, Schreyer D, Treanor-Taylor M, et al. The landscape of long noncoding RNA during cutaneous squamous cell carcinoma progression. Br J Dermatol 2025;193:490-501. [Crossref] [PubMed]
6. Zakhem GA, Pulavarty AN, Carucci J, et al. Association of Patient Risk Factors, Tumor Characteristics, and Treatment Modality With Poor Outcomes in Primary Cutaneous Squamous Cell Carcinoma: A Systematic Review and Meta-analysis. JAMA Dermatol 2023;159:160-71. [Crossref] [PubMed]
7. Petzold A, Steeb T, Wessely A, et al. Comparative efficacy analysis identifies immune checkpoint blockade as a new survival benchmark in advanced cutaneous squamous cell carcinoma. Eur J Cancer 2022;170:42-53. [Crossref] [PubMed]
8. Ding LJ, Jiang X, Li T, et al. Role of UFMylation in tumorigenesis and cancer immunotherapy. Front Immunol 2024;15:1454823. [Crossref] [PubMed]
9. Pan S, Chen R. Pathological implication of protein post-translational modifications in cancer. Mol Aspects Med 2022;86:101097. [Crossref] [PubMed]
10. Zhang H, Yan Q, Jiang S, et al. Protein post-translational modifications and tumor immunity: A pan-cancer perspective. Phys Life Rev 2025;55:142-209. [Crossref] [PubMed]
11. Hermann J, Schurgers L, Jankowski V. Identification and characterization of post-translational modifications: Clinical implications. Mol Aspects Med 2022;86:101066. [Crossref] [PubMed]
12. Ji YZ, Jia LL, Liu SR. Inflammation and epigenetics of sporotrichosis disease. Semin Cell Dev Biol 2024;154:193-8. [Crossref] [PubMed]
13. Dauch C, Shim S, Cole MW, et al. KMT2D loss drives aggressive tumor phenotypes in cutaneous squamous cell carcinoma. Am J Cancer Res 2022;12:1309-22.
14. Wang R, Li C, Cheng Z, et al. H3K9 lactylation in malignant cells facilitates CD8(+) T cell dysfunction and poor immunotherapy response. Cell Rep 2024;43:114686. [Crossref] [PubMed]
15. Jing F, Zhu L, Zhang J, et al. Multi-omics reveals lactylation-driven regulatory mechanisms promoting tumor progression in oral squamous cell carcinoma. Genome Biol 2024;25:272. [Crossref] [PubMed]
16. Zhu W, Fan C, Hou Y, et al. Lactylation in tumor microenvironment and immunotherapy resistance: New mechanisms and challenges. Cancer Lett 2025;627:217835. [Crossref] [PubMed]
17. Viiklepp K, Nissinen L, Ojalill M, et al. C1r Upregulates Production of Matrix Metalloproteinase-13 and Promotes Invasion of Cutaneous Squamous Cell Carcinoma. J Invest Dermatol 2022;142:1478-1488.e9. [Crossref] [PubMed]
18. Dousset L, Khalife F, Chalopin-Fillot D, et al. Metabolism-based scoring of cutaneous squamous cell carcinoma to predict tumour features and responses to treatment with dihydroorotate dehydrogenase inhibitors. Br J Dermatol 2025;193:936-47. [Crossref] [PubMed]
19. Xing J, Chen M, Han Y. Multiple datasets to explore the tumor microenvironment of cutaneous squamous cell carcinoma. Math Biosci Eng 2022;19:5905-24. [Crossref] [PubMed]
20. Zhang P, Wang D, Zhou G, et al. Novel post-translational modification learning signature reveals B4GALT2 as an immune exclusion regulator in lung adenocarcinoma. J Immunother Cancer 2025;13:e010787. [Crossref] [PubMed]
21. Ji AL, Rubin AJ, Thrane K, et al. Multimodal Analysis of Composition and Spatial Architecture in Human Squamous Cell Carcinoma. Cell 2020;182:497-514.e22. [Crossref] [PubMed]
22. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
23. Simon N, Friedman J, Hastie T, et al. Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent. J Stat Softw 2011;39:1-13. [Crossref] [PubMed]
24. Hu J, Nie K, Wang X. Novel ubiquitination-related biomarkers for Crohn's disease identified by multi-omics study and experimental validation. Front Immunol 2025;16:1687606. [Crossref] [PubMed]
25. South AP, Laimer M, Gueye M, et al. Type VII Collagen Deficiency in the Oncogenesis of Cutaneous Squamous Cell Carcinoma in Dystrophic Epidermolysis Bullosa. J Invest Dermatol 2023;143:2108-19. [Crossref] [PubMed]
26. Zhang X, Chen T, Zhang F, et al. METTL1 coordinates cutaneous squamous cell carcinoma progression via the m7G modification of the ATF4 mRNA. Cell Death Discov 2025;11:27. [Crossref] [PubMed]
27. Semenza GL. HIF-1: mediator of physiological and pathophysiological responses to hypoxia. J Appl Physiol 1985;2000:1474-80. [Crossref] [PubMed]
28. Semenza GL. Targeting HIF-1 for cancer therapy. Nat Rev Cancer 2003;3:721-32. [Crossref] [PubMed]
29. Shen LI, Liu L, Yang Z, et al. Identification of genes and signaling pathways associated with squamous cell carcinoma by bioinformatics analysis. Oncol Lett 2016;11:1382-90. [Crossref] [PubMed]
30. Mukhopadhyay D, Chakraborty B, Sarkar S, et al. Clinical implications of activation of the LIMD1-VHL-HIF1α pathway during head-&-neck squamous cell carcinoma development. Indian J Med Res 2024;159:479-93. [Crossref] [PubMed]
31. Forsythe JA, Jiang BH, Iyer NV, et al. Activation of vascular endothelial growth factor gene transcription by hypoxia-inducible factor 1. Mol Cell Biol 1996;16:4604-13. [Crossref] [PubMed]
32. Ji J, Wang Y, Jing A, et al. HIF1A-dependent overexpression of MTFP1 promotes lung squamous cell carcinoma development by activating the glycolysis pathway. Heliyon 2024;10:e28440. [Crossref] [PubMed]
33. Lv Z, Liu RD, Chen XQ, et al. HIF-1α promotes the stemness of oesophageal squamous cell carcinoma by activating the Wnt/β-catenin pathway. Oncol Rep 2019;42:726-34. [Crossref] [PubMed]
34. Chen X, Liu HY, Zhou WB, et al. Hypoxia-inducible factor 1-alpha and lactate dehydrogenase-A axis in metabolic changes and aggression in esophageal squamous-cell carcinoma. World J Gastrointest Oncol 2025;17:103450. [Crossref] [PubMed]
35. Yang Z, Zhang C, Feng Y, et al. Tenascin-C is involved in promotion of cancer stemness via the Akt/HIF1ɑ axis in esophageal squamous cell carcinoma. Exp Mol Pathol 2019;109:104239. [Crossref] [PubMed]
36. Wu K, Mao YY, Chen Q, et al. Hypoxia-induced ROS promotes mitochondrial fission and cisplatin chemosensitivity via HIF-1α/Mff regulation in head and neck squamous cell carcinoma. Cell Oncol (Dordr) 2021;44:1167-81. [Crossref] [PubMed]
37. Bai M, Xu P, Cheng R, et al. ROS-ATM-CHK2 axis stabilizes HIF-1α and promotes tumor angiogenesis in hypoxic microenvironment. Oncogene 2025;44:1609-19. [Crossref] [PubMed]
38. Yang S, Zhan Q, Su D, et al. HIF1α/ATF3 partake in PGK1 K191/K192 succinylation by modulating P4HA1/succinate signaling in glioblastoma. Neuro Oncol 2024;26:1405-20. [Crossref] [PubMed]
39. Cui J, Li G, Yin J, et al. GSTP1 and cancer: Expression, methylation, polymorphisms and signaling Int J Oncol 2020;56:867-78. (Review). [Crossref] [PubMed]
40. Kanwal R, Pandey M, Bhaskaran N, et al. Protection against oxidative DNA damage and stress in human prostate by glutathione S-transferase P1. Mol Carcinog 2014;53:8-18. [Crossref] [PubMed]
41. Sun X, Guo C, Huang C, et al. GSTP alleviates acute lung injury by S-glutathionylation of KEAP1 and subsequent activation of NRF2 pathway. Redox Biol 2024;71:103116. [Crossref] [PubMed]
42. Wang SQ, Chen JJ, Jiang Y, et al. Targeting GSTP1 as Therapeutic Strategy against Lung Adenocarcinoma Stemness and Resistance to Tyrosine Kinase Inhibitors. Adv Sci (Weinh) 2023;10:e2205262. [Crossref] [PubMed]
43. Yin S, Zhao S, Li J, et al. NUMA1 modulates apoptosis of esophageal squamous cell carcinoma cells through regulating ASK1-JNK signaling pathway. Cell Mol Life Sci 2023;80:211. [Crossref] [PubMed]
44. Mutallip M, Nohata N, Hanazawa T, et al. Glutathione S-transferase P1 (GSTP1) suppresses cell apoptosis and its regulation by miR-133α in head and neck squamous cell carcinoma (HNSCC). Int J Mol Med 2011;27:345-52. [Crossref] [PubMed]
45. Rajesh D, Balakrishna S, Azeem Mohiyuddin SM, et al. Novel association of oral squamous cell carcinoma with GSTP1 Arg187Trp gene polymorphism. J Cell Biochem 2019;120:5906-12. [Crossref] [PubMed]
46. Rajesh D, Balakrishna S, Azeem Mohiyuddin SM, et al. GSTP1 c.341C>T gene polymorphism increases the risk of oral squamous cell carcinoma. Mutat Res Genet Toxicol Environ Mutagen 2018;831:45-9. [Crossref] [PubMed]
47. Pincinato EC, Costa EFD, Lopes-Aguiar L, et al. GSTM1, GSTT1 and GSTP1 Ile105Val polymorphisms in outcomes of head and neck squamous cell carcinoma patients treated with cisplatin chemoradiation. Sci Rep 2019;9:9312. [Crossref] [PubMed]
48. FeiFei W. FBX8 degrades GSTP1 through ubiquitination to suppress colorectal cancer progression. Cell Death Dis 2019;10:351. [Crossref] [PubMed]
49. Zhang W, Dai J, Hou G, et al. SMURF2 predisposes cancer cell toward ferroptosis in GPX4-independent manners by promoting GSTP1 degradation. Mol Cell 2023;83:4352-4369.e8. [Crossref] [PubMed]
50. Sun Y, He Q, Li J, et al. A GSTP1-mediated lactic acid signaling promotes tumorigenesis through the PPP oxidative branch. Cell Death Dis 2023;14:463. [Crossref] [PubMed]
51. Singh S, Okamura T, Ali-Osman F. Serine phosphorylation of glutathione S-transferase P1 (GSTP1) by PKCα enhances GSTP1-dependent cisplatin metabolism and resistance in human glioma cells. Biochem Pharmacol 2010;80:1343-55. [Crossref] [PubMed]
52. Corcoran SE, O'Neill LA. HIF1α and metabolic reprogramming in inflammation. J Clin Invest 2016;126:3699-707. [Crossref] [PubMed]
53. Balakrishnan CK, Tye GJ, Balasubramaniam SD, et al. CD74 and HLA-DRA in Cervical Carcinogenesis: Potential Targets for Antitumour Therapy. Medicina (Kaunas) 2022;58:190. [Crossref] [PubMed]
54. Tran TT, Sánchez-Zuno GA, Osmani L, et al. Improving immunotherapy responses by dual inhibition of macrophage migration inhibitory factor and PD-1. JCI Insight 2025;10:e191539. [Crossref] [PubMed]
55. Chen Z, Liu M, Wang N, et al. Unleashing the Potential of Camptothecin: Exploring Innovative Strategies for Structural Modification and Therapeutic Advancements. J Med Chem 2024;67:3244-73. [Crossref] [PubMed]
56. Bertozzi D, Marinello J, Manzo SG, et al. The natural inhibitor of DNA topoisomerase I, camptothecin, modulates HIF-1α activity by changing miR expression patterns in human cancer cells. Mol Cancer Ther 2014;13:239-48. [Crossref] [PubMed]
57. Ren B, Kwah MX, Liu C, et al. Resveratrol for cancer therapy: Challenges and future perspectives. Cancer Lett 2021;515:63-72. [Crossref] [PubMed]
58. Srivani G, Behera SK, Dariya B, et al. Resveratrol binds and inhibits transcription factor HIF-1α in pancreatic cancer. Exp Cell Res 2020;394:112126. [Crossref] [PubMed]
59. Brockmueller A, Girisa S, Kunnumakkara AB, et al. Resveratrol Modulates Chemosensitisation to 5-FU via β1-Integrin/HIF-1α Axis in CRC Tumor Microenvironment. Int J Mol Sci 2023;24:4988. [Crossref] [PubMed]
60. Alavi M, Farkhondeh T, Aschner M, et al. Resveratrol mediates its anti-cancer effects by Nrf2 signaling pathway activation. Cancer Cell Int 2021;21:579. [Crossref] [PubMed]
61. Verneau J, Sautés-Fridman C, Sun CM. Dendritic cells in the tumor microenvironment: prognostic and theranostic impact. Semin Immunol 2020;48:101410. [Crossref] [PubMed]
62. Li R, Fang F, Jiang M, et al. STAT3 and NF-κB are Simultaneously Suppressed in Dendritic Cells in Lung Cancer. Sci Rep 2017;7:45395. [Crossref] [PubMed]
63. Chen J, Duan Y, Che J, et al. Dysfunction of dendritic cells in tumor microenvironment and immunotherapy. Cancer Commun (Lond) 2024;44:1047-70. [Crossref] [PubMed]
64. Yang Y, Cao Y. The impact of VEGF on cancer metastasis and systemic disease. Semin Cancer Biol 2022;86:251-61. [Crossref] [PubMed]
65. Ma S, Lu CC, Yang LY, et al. ANXA2 promotes esophageal cancer progression by activating MYC-HIF1A-VEGF axis. J Exp Clin Cancer Res 2018;37:183. [Crossref] [PubMed]
66. Zhao H, Wu X, Wang Y, et al. Histone variant H2AZ1 drives lung cancer progression through the RELA-HIF1A-EGFR signaling pathway. Cell Commun Signal 2024;22:453. [Crossref] [PubMed]
67. Arndt PF, Turkowski K, Cekay MJ, et al. Endothelin and the tumor microenvironment: a finger in every pie. Clin Sci (Lond) 2024;138:617-34. [Crossref] [PubMed]
68. Sestito R, Tocci P, Roman C, et al. Functional interaction between endothelin-1 and ZEB1/YAP signaling regulates cellular plasticity and metastasis in high-grade serous ovarian cancer. J Exp Clin Cancer Res 2022;41:157. [Crossref] [PubMed]