Crosstalk Between Immunity and Oncogenes Within the Tumor Microenvironment of HPV-Associated Cervical Squamous Cell Carcinoma

Abstract

Introduction: Every two minutes, a woman dies from cervical cancer, which is considered the fourth most common cancer among women worldwide. The dynamic interplay between tumor inflammation, immune crosstalk, oncogenes, and tumor suppressor genes plays a crucial role in tumor development and progression.

Methods: Using clinical and integrated bioinformatics, the mRNA expression pattern of 168 immune and tumor-related genes in the tumor microenvironment (TME) of HPV-positive cervical squamous cell carcinoma (CSCC) was analyzed.

Results: The study identified 94 DEGs, of which 55 genes were remarkably upregulated, including CASP8, ZHX2, BCL2L1, CTNNB1, RB1, BAX, CD274, CCL20, FOXP3, and CCL18. The top three-fold changes were associated with CASP8, ZHX2, and BCL2L1, respectively. In contrast, downregulation was discovered for 39 genes associated with immunity, regulation of cell cycle, and DNA damage response (HRAS, CCND1, ATM, CXCR1, and MIF). Gene-gene interaction and correlation analysis showed positive correlations, including RB1 and CASP8, RB1 and BCL2L1, and CCL20 with CCL18. Notably, six genes exhibited increased expression and showed a strong correlation with enhanced overall survival (OS) and disease-free survival (DFS), indicating their potential utility as prognostic biomarkers. Upregulated genes were positively associated with various immune cells, including B cells, CD8+ and CD4+ T cells, macrophages, neutrophils, and dendritic cells. Functional enrichment analysis revealed involvement in cancer-related processes, inflammatory responses, and cell migration, with key pathways linked to cytokine signaling and chemokine receptor interactions.

Discussion: Through the integration of clinical, experimental, and computational analyses, potential therapeutic targets and prognostic biomarkers were identified that may help improve clinical outcomes. Future studies should focus on the functional assays of identified genes both in vitro and in vivo.

Keywords: CSCC; HPV; TME; cervical squamous cell carcinoma; human papillomavirus; tumor microenvironment.

Graphical abstract

Figure 1

Heatmap Gene-Gene Interactions and Correlation Analysis of Gene Expression. Correlations are color-coded, with red indicating positive and blue indicating negative associations. The intensity of the color corresponds to the magnitude of the correlation coefficient. Each cell in the matrix displays the correlation coefficient between the corresponding gene pairs.

Figure 2

Boxplots of the Most Highly Upregulated Genes, illustrating the expression levels in tumor, normal and metastatic cervical tissues. The Kruskal. Wallis’ p-value was utilized to assess statistical significance. The bars illustrate the proportion of tumor and metastatic tissue samples with elevated gene expression compared to normal tissues, evaluated at multiple quantile thresholds (minimum, 1st quartile, median, 3rd quartile, and maximum).

Figure 3

The prognostic value of highly upregulated genes. (A-J) Kaplan-Meier curves of survival outcomes for each gene by GEPIA2. The p-values were calculated using the Mantel-Cox test to assess statistical significance. Hazard ratios (HR) were calculated from the Cox proportional hazards model. HR >1 indicates that high gene expression is associated with a reduced survival rate, whereas HR < 1 suggests that high gene expression correlates with an increased survival rate.

Figure 4

Gene set variation analysis (GSVA) of the most upregulated gene. The heatmap presents the relationship between GSVA scores and the activity of cancer-associated pathways in CESC. The GSVA score reflects the overall expression level of a gene set and is directly correlated with individual gene expression levels. *P value ≤ 0.05; ^#FDR ≤ 0.05.

Figure 5

Top 6 Enriched Ontology Clusters with their representative enriched terms analyzed by Metascape. (A) Enriched ontology clusters were color-coded based on their cluster IDs, with nodes belonging to the same cluster generally positioned in close proximity. (B) The thickness of the edge represents the similarity score.

Figure 6

Functional enrichment analysis of 21 with similar expression levels in CSCC using g: Profiler. Graphs (A–C) showed the results of the enriched terms with the statistical significance (p-value) according to the (A) gene ontology biological processes (GO BP), (B) gene ontology molecular functions (GO MF), and (C) KEGG.

Figure 7

Topology network analysis illustrating protein-protein interactions (PPI), networks encompassing various biological entities such as chemical, complex, miRNA, phenotype, protein, protein family, small molecule, and stimulus using NetworkAnalyst. The central panel displays the network with distinct modules highlighted in various colors. MicroRNAs are represented by blue squares, while genes are depicted as red circles. The size of each circle corresponds to the degree of the node within the network. This analysis comprises a total of 859 nodes and 1196 edges. Nodes in the network represent proteins and biological entities, while edges represent the interactions between them.

Figure 8

Correlation analysis between (A) GDSC and (B) CTRP drug sensitivity datasets and the mRNA levels of the ten most highly co-upregulated genes across various cancers. Spearman correlation indicates the relationship between gene expression and drug response. Positive correlation suggests that elevated gene expression is associated with increased drug resistance, whereas a negative correlation implies that higher expression corresponds to greater drug sensitivity. The figure summarizes the correlation between the ten most co-upregulated gene expression levels and the sensitivity of the top 50 drugs in pan-cancer.