Metabolic Molecule PLA2G2D Is a Potential Prognostic Biomarker Correlating With Immune Cell Infiltration and the Expression of Immune Checkpoint Genes in Cervical Squamous Cell Carcinoma

Abstract

Cervical squamous cell carcinoma (CSCC) is the major pathological type of cervical cancer (CC), the second most prevalent reproductive system malignant tumor threatening the health of women worldwide. The prognosis of CSCC patients is largely affected by the tumor immune microenvironment (TIME); however, the biomarker landscape related to the immune microenvironment of CSCC and patient prognosis is less characterized. Here, we analyzed RNA-seq data of CSCC patients from The Cancer Genome Atlas (TCGA) database by dividing it into high- and low-immune infiltration groups with the MCP-counter and ESTIMATE R packages. After combining weighted gene co-expression network analysis (WGCNA) and differentially expressed gene (DEG) analysis, we found that PLA2G2D, a metabolism-associated gene, is the top gene positively associated with immune infiltration and patient survival. This finding was validated using data from The Cancer Genome Characterization Initiative (CGCI) database and further confirmed by quantitative reverse transcription-polymerase chain reaction (qRT-PCR). Finally, multiplex immunohistochemistry (mIHC) was performed to confirm the differential infiltration of immune cells between PLA2G2D-high and PLA2G2D-low tumors at the protein level. Our results demonstrated that PLA2G2D expression was significantly correlated with the infiltration of immune cells, especially T cells and macrophages. More importantly, PLA2G2D-high tumors also exhibited higher infiltration of CD8⁺ T cells inside the tumor region than PLA2G2D-low tumors. In addition, PLA2G2D expression was found to be positively correlated with the expression of multiple immune checkpoint genes (ICPs). Moreover, based on other immunotherapy cohort data, PLA2G2D high expression is correlated with increased cytotoxicity and favorable response to immune checkpoint blockade (ICB) therapy. Hence, PLA2G2D could be a novel potential biomarker for immune cell infiltration, patient survival, and the response to ICB therapy in CSCC and may represent a promising target for the treatment of CSCC patients.

Keywords: PLA2G2D; cervical squamous cell carcinoma; immune infiltration; metabolism; multiplex immunohistochemistry; tumor immune microenvironment.

Figure 1

Immune classification and survival analysis of 250 CSCC patients selected from TCGA database. (A) Heatmap of TME composition as defined by the MCP-counter algorithm. (B) Kaplan–Meier survival curve for immune infiltration-high and infiltration-low groups determined by the MCP-counter. (C) Heatmap showing immune and stromal fraction by the ESTIMATE algorithm. (D, E) Kaplan–Meier survival curve for patients with high and low fractions of immune cells and stromal cells classified by the ESTIMATE algorithm.

Figure 2

Gene modules construction by weighted gene co-expression network analysis (WGCNA). (A) Sample dendrogram and trait heatmap. (B) The relationship of soft threshold with scale independence and mean connectivity. (C) Clustering of module eigengenes, height set at 0.5 to merge modules. (D) Cluster dendrogram showing the hierarchical cluster tree for the identified co-expression gene modules with different colors.

Figure 3

Key modules and genes identified based on module–trait relationship. (A) The relationship between co-expression modules and traits; each grid includes the degree of correlation and P-value. (B, C) Scatter plots showing the relationship of module membership (MM) in brown module, with gene significance (GS) for high-immune infiltration determined by the MCP-counter and ESTIMTE algorithms, for which MM >0.8 and GS >0.5 were set as gene filtered standard, respectively. (D) Venn plot showing common filtered genes identified in the high-immune infiltration groups determined by the MCP-counter and ESTIMATE. (E) GO analysis for WGCNA filtered genes in the biological process.

Figure 4

The expression of PLA2G2D is positively correlated with immune infiltration and patient survival. Volcano plot for DEGs identified for high- and low-immune infiltration groups based on the MCP-counter (A) and ESTIMATE (B). The green and red dots represented the significantly downregulated and upregulated genes, respectively, and the gray dots represented the undifferentiated genes; the purple dots showed the WGCNA filtered genes, with |Log₂(FoldChange)| >1 and adjusted P-value <0.01. (C) Cytoscape network plot showed the relationship between PLA2G2D and other co-expressed genes with weight value >0.2. The purple nodes showing the co-expressed genes with weight value >0.3. (D) GO pathways enrichment analysis for genes co-expressed with PLA2G2D. (E) Histogram showing the relationship between the expression levels of PLA2G2D and other five co-expressed genes, ****P < 0.0001. (F) Kaplan–Meier survival curve plot showing the relationship between PLA2G2D expression level and patient survival. (G) Time-dependent ROC curve analysis for minus PLA2G2D TPM value at 2-, 3-, and 5-year cutoffs.

Figure 5

Further validation of the relationship between PLA2G2D expression and immune infiltration. (A) Heatmap for PLA2G2D high- and low-expression groups based on the MCP-counter, ESTIMATE, EPIC, and quanTIseq algorithms using data from the CGCI database. (B) Violin plot showing the differential scores of six immune cell types determined by the MCP-counter algorithm. (C) Violin plot showing the differential StromalScore, ImmuneScore, and ESTIMATEScore determined by the ESTIMATE algorithm. (D) Violin plot showing the differential immune cell types determined by quanTIseq and EPIC algorithms. (E, F) Histograms showing the relationship between the expression levels of PLA2G2D and WGCNA filtered five co-expressed genes using data from the CGCI database (D) and freshly isolated clinical samples (E). *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.

Figure 6

mIHC staining for CD3, CD8, CD68, PCK, and DAPI in cervical squamous cancer specimens with high (A–C) and low expression (D–F) of the PLA2G2D gene. (A, D) Merged mIHC images at low magnification (×4). (B, E) Merged mIHC images at high magnification (×20) indicating the same filed in (A, D). (C, F) Single spectral images indicating the same filed in (B, E).

Figure 7

Statistical analysis for mIHC staining. (A) Pie plots showing the composition of immune cells within different regions in PLA2G2D high- and low-expression groups. (B) Histograms showing the comparison of the percentage of individual immune cell populations between PLA2G2D high- and low- expression groups in different regions. (C) Histograms showing the comparison of the density of individual immune cell populations between PLA2G2D high- and low-expression groups in different regions. (D) Histogram showing the comparison of the CD8/CD68 ratio between PLA2G2D high- and low-expression groups in different regions. *P < 0.05, **P < 0.01, and ***P < 0.001.

Figure 8

Correlation of the expression level of PLA2G2D and ICPs. Histograms showing the comparison of ICP expression level between PLA2G2D-high and PLA2G2D-low groups in TCGA (A) and CGCI database (C). Pearson correlation analysis for the expression level of PLA2G2D and ICPs in TCGA (B) and CGCI (D) database. *P < 0.05, **P < 0.01, and ****P < 0.0001.

Figure 9

Correlation of the expression level of PLA2G2D and response to ICB treatment in melanoma (A–G) and urothelial carcinoma immunotherapy cohorts (H–L). (A, B) Bar and boxplot for the relationship between PLA2G2D expression and response to ICB treatment in melanoma cohort. (C) Heatmap of the relationship of ICPs, cytotoxic genes, and CYT score with different PLA2G2D expression levels and response to treatment in paired pre- and on-treatment patients. Boxplots showing the differential expression of ICPs, cytotoxic genes, and CYT scores between patients with different PLA2G2D expression levels before ICB therapy (D, E) and after ICB therapy (F, G). (H, I) Bar and boxplot for the relationship between PLA2G2D expression and response to ICB treatment in urothelial carcinoma cohort. (J) Heatmap of the relationship of ICPs, cytotoxic genes, and CYT score with different PLA2G2D expression levels. (K, L) Boxplots showing the differential expression of ICPs, cytotoxic genes, and CYT score between patients with different PLA2G2D expression levels. *P < 0.05, **P < 0.01, and ***P < 0.001.