key: cord-0022754-hypc8nd4
authors: Liu, Hong; Xu, Ruiyi; Gao, Chun; Zhu, Tong; Liu, Liting; Yang, Yifan; Zeng, Haihong; Huang, Yafei; Wang, Hui
title: 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
date: 2021-10-18
journal: Front Oncol
DOI: 10.3389/fonc.2021.755668
sha: c056fea07813290cdf1cf9a2de87a3d7b62e8178
doc_id: 22754
cord_uid: hypc8nd4

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.

Cervical cancer (CC) is the fourth common malignant disease for female with an estimated 604,000 new cases and 342,000 deaths worldwide in 2020 (1) . Cervical squamous cell carcinoma (CSCC) and adenocarcinoma are the most common pathological types accounting for approximately 70% and 25% of all CC, respectively (2) . Among them, CSCC is mainly related to human papillomavirus (HPV) 16 subtype infection and viral gene integration, whereas adenocarcinoma is often complicated with HPV18 infection (3) . Despite the substantial efforts being made in promoting HPV vaccination and early screening, the incidence of CC remains alarming in developing countries (2) . In fact, CC is the most frequently diagnosed female cancer in 23 countries according to the latest report (2) . Current therapies for CC including surgical treatment, chemotherapy, and radiotherapy have greatly improved the clinical outcome of CC; however, the therapeutic efficacy remains limited for patients with advanced and distant conditions, which estimated a median overall survival of 17 months and 5-year survival of 17% (4, 5) . Therefore, there is an urgent need for developing novel therapeutic strategies that can effectively treat these patients (6, 7) .

Immunotherapy is one of such strategies that has become a rapidly developing field for cancer treatment including cervical cancer. Through replenishing a sufficient number of expanded autologous T cells that can specifically kill cancer cells, adoptive cell transfer (ACT) has shown great promise in treating CC patients with advance diseases. Two major strategies of ACT, tumor infiltrating lymphocyte (TIL) and T-cell receptorengineered T cell (TCR-T), were reported to have the objective response rate (ORR) of 44.4% and 50% (8, 9) , respectively. Immune checkpoint inhibitors (ICBs), another type of immunotherapy that target programmed death-1(PD-1) receptor or ligand and cytotoxic T-lymphocyte-associated protein 4 (CTLA4), have achieved great success in various kinds of tumors. However, in several clinical trials for advanced cervical cancer patients, the response rate to ICBs was relatively low (10) (11) (12) (13) . Previous studies have shown that multiple factors are correlated with the efficacy of ICB therapy. Among them, more immune cell infiltration, especially the density and localization of CD8 + T cells in the tumor immune microenvironment (TIME) has been demonstrated to be correlated with a favorable response to ICBs in various cancer types (14, 15) . However, the prognostic value and underlying molecular mechanism of immune infiltration in CC with or without immunotherapy remain less characterized.

Here, we utilized the RNA-seq data of CSCC patients from The Cancer Genome Atlas (TCGA) database for the sake of finding biomarkers related to prognosis and immune cell infiltration. To this end, we used two immune scoring algorithms to divide them into two groups with high-and lowimmune infiltration. By applying weighted gene co-expression network analysis (WGCNA) and differentially expressed gene (DEG) analysis, phospholipase A 2 Group IID (PLA2G2D), an immune-and metabolism-associated molecule, was identified to be the biomarker which is predictive for patient prognosis and immune cell infiltration of CC. Furthermore, five immunerelated genes (i.e., SLAMF6, SLAMF1, SH2D1A, TRAT1, and ZNF831) were found to be co-expressed with PLA2G2D. In addition, PLA2G2D expression was also found to be positively correlated with the expression of multiple ICP genes. Finally, through a series of bioinformatics analysis and experimental verification approaches at both the transcriptional and protein levels, we proved that PLA2G2D could be a novel biomarker correlating with immune infiltration, especially CD8 + T cells, in CSCC.

The Cancer Genome Atlas (TCGA) CSCC RNA-seq gene expression matrix based on Illumina platform and phenotype data were downloaded from UCSC Xena (https://xenabrowser. net/datapages/). Two hundred fifty CSCC patients were selected for subsequent analysis according to pathological diagnosis of squamous cell carcinoma. Two types of data including raw counts and fragments per kilo base per million mapped (FPKM) reads were applied. FPKM data were converted to transcripts per kilobase of per million mapped (TPM) data. In addition, the validation cohort dataset also based on Illumina high-throughput platform was downloaded from The Cancer Genome Characterization Initiative (CGCI) database (16) . RNAseq data from two immunotherapy cohorts for melanoma (17) and urothelial cancer (18) were downloaded from a public database. The gene names of all expression matrix were deannotated through "ClusterProfiler" and "org.Hs.eg.db" package based on R language (V4.1.1) (19) .

Estimation of Tumor Immune Infiltration and Cytolytic Activity Score cell infiltration in tumor, based on which samples were classified into high-and low-immune infiltration groups using the hierarchical clustering method (20) . Alternatively, the level of immune and stromal fraction was scored by Estimation of Stromal Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) based on log 2 -transformed TPM data, and samples were further divided into high-and lowimmune infiltration groups equally with the cutoff score set at median (21) . In addition, EPIC (22) and quanTIseq (23) algorithm were also utilized to calculate immune cell proportion through TPM expression matrix. CYT score was calculated by the log-transformed geometric mean of GZMA and PRF1 TPM value (24) .

The log 2 -transformed TPM value expression matrix was put into WGCNA R package to select immune-associated genes (25) . Firstly, 12,417 genes with coefficient of variation (CV) >0.5 were selected for further analysis by WGCNA R package. Of note, five samples were excluded for the abnormal height value in sample dendrogram analysis. The power of b value was set at 3 to construct the topological overlap matrix (TOM). Next, we set the minimum module gene size at 30 and generated 63 gene modules with different colors based on hierarchical clustering method, then merged multiple similar gene modules into one. Finally, the correlation between gene modules and traits was calculated to determine the most relevant module and WGCNA filtered genes were screened by setting the standard of module membership (MM) >0.8 and gene significance (GS) >0.5. Genes network was constructed by Cytoscape (V3.8.0) based on the genes with weight value >0.20 (26) .

In order to determine the DEGs between the high-and lowimmune infiltration groups, the raw read count matrix of 250 CSCC patients from TCGA database was brought into DESeq2 (27) . Genes with |Log 2 (FoldChange)| >1 and adjusted P-value <0.01 were considered as DEGs. ClusterProfiler R package was used to perform Gene Ontology (GO) term enrichment analysis for biological pathway. CGCI data matrix was imported into GSEA software (V4.1.0) for gene set enrichment analysis using the Hallmark gene sets database. 

All FFPE slides were deparaffinized by dipping in xylene for 1 h and then rehydrated using the gradient ethanol method. After deparaffinization and rehydration, all slides were put into boiled AR6 retrieval solution for heat-induced epitope retrieval in a microwave for 15 min. Endogenous peroxidase was eliminated with 3% H 2 O 2 for 15 min. Slides were cooled to room temperature followed by washing with 1× Tris-buffered saline-Tween 20 (TBST) buffer. Then, Opal 7-color manual IHC kit (Akoya Biosciences, NEL811001KT) was used to stain several markers in a single FFPE slide. To block non-specific protein binding sites, slides were incubated with blocking buffer (Antgene, ANT041) at room temperature for 10 min. Subsequently, FFPE tissue slides were incubated with primary antibody at 37°C for 30 min or 4°C overnight, HRP-labeled secondary antibody for 10 min, and Opal fluorescein for 10 min, successively. To stain four markers in a single slide, four rounds of the above staining procedure were performed with indicated primary antibody and matched Opal fluorescein pairs (Supplementary Table 2 ). The slides were always washed three times between each step with 1× TBST buffer. After four rounds, samples were stained for cellular DNA with 4′,6-diamidino-2phenylindole (DAPI) (1:10, Servicebio, G1012-10ML) for 10 min followed by mounting with Fluoromount-G (SouthernBiotech, 0100-01) and preserved in the dark at 4°C. All slides were scanned by Verctra V3.0 System, and ×4 and ×20 objective lens were used to acquire low-power and high-power images, respectively. Images were analyzed with Phenochart V1.0 software and inForm V2.4 software for tissue component segmentation of tumor, stroma, and glass regions, respectively. To identify cell phenotypes, a threshold of fluorescein intensity was set for each marker. The ratio of the targeted cell counts and all cell counts in a certain region was used to calculate the percentage of immune cells in each region. The ratio of cell counts and area (cell/mm 2 ) was utilized to reveal the density of a certain cell population.

Survival information of TCGA CSCC patients was downloaded from the UCSC Xena database. Patients who survived less than 30 days were excluded from this analysis. Kaplan-Meier survival analysis was applied to educe the correlation between 5-year overall survival and the expression level of indicated genes. Based on the relationship with survival time, survival status, and minus PLA2G2D TPM value, ROC curve analysis was applied by survivalROC R package. Statistical analysis was performed with R statistical package. Wilcoxon rank-sum test was employed for comparison between two groups.

P-values <0.05 were considered significant: *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.

The RNA-seq data and clinical information of a CSCC cohort including 250 patients were downloaded from TCGA database for bioinformatics analysis. Firstly, the MCP-counter, a widely used algorithm, was used to evaluate the level of immune infiltration for each sample, according to which CSCC patients were divided into high-and low-immune infiltration clusters based on the unsupervised stratification method ( Figure 1A ). Subsequently, the 5-year overall survival of the high-and lowimmune infiltration groups was compared. The results demonstrated that the high-immune infiltration cluster determined by the MCP-counter had significantly higher survival than the low-immune infiltration cluster in CSCC ( Figure 1B) . Next, ESTIMATE, another algorithm which can similarly stratify CSCC patients into high-and low-immune fraction clusters by a median cutoff ( Figure 1C) , was chosen to verify this finding. Again, the high-immune infiltration cluster classified by ESTIMATE was found to have higher survival than the low-immune infiltration cluster in CSCC ( Figure 1D ). Of note, ESTIMATE can also classify CSCC patients into high-and low-stromal fraction clusters ( Figure 1C ). However, no difference was found between the two stromal fractionstratified clusters in terms of patient survival ( Figure 1E , P = 0.22). Together, our results demonstrated that immune infiltration, but not stromal fraction, was correlated with patient survival in CSCC. Moreover, the other two algorithms, EPIC and quanTIseq, were used to calculate several immune cell proportions compared with the MCP-counter and ESTIMATE. Moreover, the immune cell infiltration results of EPIC and quanTIseq were consistent with the MCP-counter and ESTIMATE ( Figure S1 ).

CSCC RNA-seq data from TCGA was imported into WGCNA R package to identify gene modules containing similarly expressed genes, especially immune-related genes. Sample clustering and trait distribution are shown in Figure 2A . Power value was set at b = 3 to build a scale-free network ( Figure 2B ). Gene modules were calculated and the cutoff height was set at 0.5 to merge similar gene modules into one ( Figure 2C ). This analysis resulted in 63 merged dynamic gene modules with each module containing dozens to more than a thousand genes ( Figure 2D , Supplementary Table 3 ).

In order to identify immune-related modules and genes, we constructed the relationship between gene modules and immune infiltration traits ( Figure 3A ). We focused on the brown module which was most relevant to immune infiltration or immune fraction determined by the MCP-counter and ESTIMATE, respectively. To further identify genes in the brown module that are most correlated with the MCP-counter immune infiltration and ESTIMATE immune fraction, the screening criteria of MM >0.8 and GS >0.5 were applied, through which 77 genes (brown module vs. MCP-counter immune infiltration) and 73 genes (brown module vs. ESTIMATE immune fraction) were filtered (Figures 3B, C). Next, Venn plot was used to reveal the overlapped filtered genes of these two groups. As shown in Figure 3D , all 73 brown module genes contained in the ESTIMATE high-immune fraction group were found to be presented in the MCP-counter high-immune infiltration group as well. Finally, GO pathways enrichment analysis for the shared 73 WGCNA genes showed that these genes were mostly involved in immune-associated pathways, such as T-cell activation, lymphocyte differentiation, and regulation of T-cell activation ( Figure 3E ).

To further explore hub genes from these 73 genes, DEG analysis was performed using DESeq2 R package. DEGs were determined between the high-and low-immune infiltration groups determined by the MCP-counter, through which 968 upregulated genes and 323 downregulated genes were found in the high-and low-immune infiltration groups, respectively ( Figure 4A) . Similarly, 1,060 upregulated genes and 650 downregulated genes were found in the high-and low-immune fraction groups determined by ESTIMATE, respectively ( Figure 4B ). Interestingly, all 73 WGCNA filtered genes were presented in these two upregulated gene sets. Notably, PLA2G2D, a metabolism-related gene, was identified as the most differentially expressed gene among the 73 genes, as determined by the log 2 (FoldChange) value of these genes ( Table 1 and Figures 4A, B) .

Next, we sought to identify the genes that were co-expressed with PLA2G2D. The co-expression network between PLA2G2D and other genes was constructed by using the criteria of weight >0.20 ( Figure 4C and Table 2 ). To gain insight into the function of these co-expressed genes, GO analysis of the biological process was conducted and the most co-expressed genes were enriched in immune-related pathways ( Figure 4D ). Of note, several genes with weight >0.3 including SLAMF6, SLAMF1, SH2D1A, TRAT1, and ZNF831 were most co-expressed with PLA2G2D ( Figure 4E) . Importantly, Kaplan-Meier survival curve analysis showed that the expression level of PLA2G2D was positively correlated with patient survival ( Figure 4F , P = 0.00017). Finally, ROC curve was constructed to demonstrate the prognostic ability of PLA2G2D in TCGA-CSCC cohort. The AUCs of PLA2G2D TPM minus value for 2, 3, and 5 years were 0.773, 0.683, and 0.639, respectively ( Figure 4G ).

To validate the correlation between overexpressed PLA2G2D and more immune cell infiltration demonstrated in Figure 1 , we analyzed RNA-seq data of CSCC from the CGCI database that was also sequenced on an Illumina high-throughput sequence platform. The immune characteristics of patients in the PLA2G2D high-and low-expression clusters are illustrated in Figure 5A by the analysis with the four algorithms.

Overexpressed PLA2G2D was bound up with more immune cell infiltration such as T cells, B cells, NK cells, and myeloid dendritic cells ( Figure 5B) . Similarly, ImmuneScore, StromalScore, and ESTIMATEScore were higher in the PLA2G2D high-expression cluster than in the PLA2G2D lowexpression cluster ( Figure 5C ). CD8 + T cells, macrophages, and B cells were also higher in the PLA2G2D-high group based on quanTIseq and EPIC algorithms ( Figure 5D ). GSEA analysis showed that the PLA2G2D high-expression cluster was enriched in several immune-associated pathways such as allograft rejection, interferon gamma response, interferon alpha response, complement, and oxidative phosphorylation pathways ( Figure S2A ), whereas PLA2G2D low-expression cluster was enriched in hypoxia, apical junction, glycolysis, mitotic spindle, and angiogenesis pathways ( Figure S2B ). We also validated whether PLA2G2D was co-expressed with other genes as indicated by the aforementioned WGCNA analysis.

Our analysis demonstrated that the expression level of PLA2G2D was also positively associated with SLAMF6, SLAMF1, SH2D1A, TRAT1, and ZNF831 when the CGCI database was used ( Figure 5E ). To further validate these results using fresh-isolated clinical samples, 18 cervical squamous carcinoma samples were collected and the corresponding clinical data are shown in Table 3 . Total RNA per sample was isolated and qRT-PCR was performed to examine the gene expression levels of PLA2G2D and other five genes ( Table 4 ). Our cohort was classified into PLA2G2D high-and low-expression groups according to the median expression level of PLA2G2D. Similarly, the expression levels of four immune-related genes were significantly higher in the PLA2G2D high-expression group, except for ZNF831, which could hardly be detected by qRT-PCR in both groups ( Figure 5F and Table 4 ).

Next, the multiplex immunohistochemistry (mIHC) method was performed to validate whether immune cell infiltration was correlated with PLA2G2D expression in our cohort of 18 CSCC clinical samples. As our aforementioned bioinformatic analysis indicated that T cells and macrophages were the most prevalent cells among the infiltrated immune cells positively correlated with PLA2G2D expression, and previous studies have shown that the number and percentage of T cells and macrophages were more than the other kind of immune cells in cervical cancer (28, 29) , we chose these two cells along with tumor cells to determine the percentage and density of individual immune cell populations in each slide. For this purpose, the mIHC panel (Supplementary Table 2 ) was designed to simultaneously stain several markers including CD3 (for all T cells), CD8 (for cytotoxic T cells), CD68 (for macrophages), PCK (for tumor cells), and DAPI (for nucleic identification). Of note, CD3 + CD8 − TILs were defined as CD4 + TILs. Our results indicated that CD3 + T cells, including both CD4 + and CD8 + TILs, and CD68 + macrophages were more abundant in the PLA2G2D-high group than in the PLA2G2D-low group, as shown by representative merged images at low magnification ( Figures 6A, D) and high magnification ( Figures 6B, E) and by single spectrum images for each marker (Figures 6C, F) . After comparing PLA2G2D-high and PLA2G2D-low groups for the compositions of individual cell populations in different regions ( Figure 7A ) and the percentages of individual cell populations in each sample ( Figure S3 ), we found that the percentage values of CD3 + TILs, CD8 + TILs, and macrophages in the stromal region and all regions of the PLA2G2D high-expression group were significantly higher compared with that of the PLA2G2D low-expression group ( Figure 7B) . Importantly, in the tumor region, CD8 + TILs were the only cell population that were more frequently observed in the PLA2G2D high-expression group ( Figure 7B ). Next, cell density (cell counts/mm 2 area) was analyzed in the tumor and stromal region, respectively. Similar to the comparison of cell percentage, cell densities of CD3 + TILs, CD8 + TILs, and macrophages were significantly higher in the stromal region of the PLA2G2D high-expression group than those of the PLA2G2D low-expression group. Moreover, in the tumor region, only CD8 + TIL density in the PLA2G2D highexpression group was higher compared with that in the PLA2G2D low-expression group ( Figure 7C) . 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. Tumor-resident macrophages were previously considered as a major cell population in suppressing the function of cytotoxic T cells and were reported to be associated with poor prognosis of cancer (30) . Therefore, we assessed the CD8/CD68 ratio between these two groups. The ratio was higher in all areas in the PLA2G2D high-expression group ( Figure 7D ).

The expression of immune checkpoint genes [e.g., PDCD1 (PD1), CD274 (PDL1), CTLA4, LAG3, and HAVCR2 (TIM3)] has been utilized in predicting the response of patients to ICB therapy in a variety of cancers including CC. To evaluate the possible relationship between the expression of PLA2G2D and these ICPs, correlation analysis was performed using transcriptomic data from TCGA database. Our analysis indicated that the expression of ICP genes was significantly higher in the PLA2G2D high-expression group compared with that in the PLA2G2D low-expression group. Through Pearson correlation analysis, we found that these common ICP genes were positively correlated with PLA2G2D expression (Figures 8A, C) . Similar findings were noted when using the CGCI validation cohort except for CD274, the expression of which was comparable in the PLA2G2D low-expression and high-expression groups and was not correlated with the expression of PLA2G2D. This discrepancy may be attributable to the small sample size in this cohort (Figures 8B, D) . Positive Correlation Between the Expression of PLA2G2D and the Response to ICB Therapy in Immunotherapy Cohorts ICB therapy has been performed in patients with advanced cervical cancer, but transcriptional data for these cohorts are not accessible. To validate the relationship between PLA2G2D and the response to ICB treatment, we analyzed the transcriptome data from immunotherapy cohorts of melanoma and urothelial carcinoma. In the melanoma anti-PD1 therapy cohort, we selected patients who have not previously received immunotherapy for enrollment analysis and divided them based on PLA2G2D expression at pretreatment stage. The percentage of compete response or partial response (CR/PR) was higher in PLA2G2D-high patients ( Figure 9A ). PLA2G2D expression level of CR/PR patients also showed a higher tendency compared with stable disease (SD) and progressive disease (PD) patients; however, these differences failed to reach statistical significance ( Figure 9B ). In this cohort, 9 and 10 pairs of pre-and ontreatment patients were assigned to the PLA2G2D-high and PLA2G2D-low cluster, respectively. Before ICB therapy, PLA2G2D-high patients show higher ICP expression and CYT score ( Figures 9C-E) . Interestingly, after ICB therapy, PLA2G2D-high patients also have higher ICP expression and CYT score (Figures 9C, F, G) . Importantly, compared with pretreatment patients, on-treatment patients have stronger cytotoxicity and higher ICP expression after anti-PD-1 therapy regardless of PLA2G2D-high or PLA2G2D-low groups ( Figure 9C) . Furthermore, in the urothelial carcinoma anti-PDL1 therapy cohort, the percentage of CR/PR patients was higher in PLA2G2D-high patients, which is similar to that in melanoma immunotherapy cohort ( Figure 9H) . CR/PR patients showed a higher tendency of PLA2G2D expression level than SD and PD patients ( Figure 9I) . Unfortunately, only pretreatment transcriptome data are available for this cohort, thus preventing us from performing similar analysis for on-treatment data. Nevertheless, PLA2G2D-high expression was also accompanied by high expression of ICPs despite no statistical significance ( Figures 9J, K) , while expression level of cytotoxicity markers (GZMA and PRF1) and the CYT score were significantly higher in the PLA2G2D-high group (Figures 9J-L) .

The tumor microenvironment (TME) is a complex milieu that comprises diverse cell populations including malignant cells, immune cells, stromal cells, and other cell types (31) . Each cell population can impact others through releasing soluble molecules and/or by direct cell-cell interaction. As a crucial part of the tumor microenvironment, TIME plays a pivotal role in tumorigenesis and the clinical outcomes of cancer patients (32) . Accumulating lines of evidence have suggested that a certain type of TIME characterized by more immune cell infiltration is indicative of a better prognosis and a favorable response to ICBs (15, 33) . Furthermore, several factors from tumor cells themselves, immune cells, and stromal cells were reported to be implicated in shaping the landscape of TIME (34) (35) (36) . However, for CSCC, the connection between TIME (especially immune infiltration) and cancer prognosis and the possible underpinning molecular mechanism are far from elucidated. This limitation greatly prevents the identification of new therapeutic targets that can be utilized to improve patient prognosis through transforming the TIME with poor-immune infiltration to high-immune infiltration. Therefore, this study was designed to address this issue. Through bioinformatics analyses, we identified a key immune-and metabolismassociated molecule PLA2G2D that was significantly related to better prognosis of CSCC patients and more immune cell infiltration. Metabolism-associated molecules are critical components of TIME (37) . Metabolic reprograming of tumor cells may promote tumor growth and metastasis directly or indirectly through influencing the functions of a variety of immune cells, especially T cells (38) (39) (40) , and the shape of TIME. Phospholipase A2 (PLA 2 ) proteins are a group of lipid metabolism-associated molecules that can catalyze the hydrolysis of the sn-2 position of membrane glycerophospholipids to release unsaturated fatty acid and lysophospholipid, thereby playing a pivotal role in multiple biochemical processes including inflammatory response. While both cytosolic and secreted forms of PLA 2 have been identified, more than one-third of the PLA 2 enzymes belong to the secreted PLA 2 (sPLA 2 ) family, which have both pro-and antiinflammatory functions as a result of the production of diverse lipid mediators (41, 42) . As a member of sPLA 2 , group IID PLA 2 (encoded by PLA2G2D), is broadly expressed in human body, such as the spleen, lymph nodes, squamous epithelium, and colorectal cancer tissue (43) (44) (45) . In the current investigation, PLA2G2D was also found to be expressed in CSCC specimens at the transcriptional level ( Figure 5E ). Previous studies have also shown that PLA2G2D was preferentially expressed in lymphoid tissue-resident dendritic cells and macrophages and implicated in anti-inflammation response in an array of inflammation-related conditions including contact hypersensitivity (43) , viral infectionassociated inflammation (46) , experimental encephalomyelitis and colitis (47) , and contact dermatitis and psoriasis (48) . The role of Pla2g2d in cancer has been investigated as well. Using both Pla2g2d-deficient and Pla2g2d-transgenic mouse models, Miki et al. found that Pla2g2d could facilitate the development of chemical-induced skin cancer accompanied by macrophage polarization toward M2 phenotype and decreased number of cytotoxic T cells (48) . The cancer-promoting effect of PLA2G2D has also been suggested in human. Recently, PLA2G2D was identified as one of the high-risk genes for colorectal and rectum adenocarcinoma that were negatively correlated with patient survival. However, PLA2G2D expression was also found to be positively correlated with immune infiltration and better prognosis in head and neck squamous cell carcinoma and breast cancer in human (49, 50) , which is consistent with our current findings in CSCC. The opposite effects of PLA2G2D in different cancer types may result from the distinct microenvironments and different downstream lipid mediators presented in each cancer type, which in turn result in either an elevated or dampened inflammatory responses, thereby leading to the distinct clinical outcomes reported in these investigations (48) . In this study, the relationship between PLA2G2D overexpression and more immune cell infiltration in CSCC was validated through two approaches. Firstly, another dataset from the CGCI database was used for validation at the transcriptional level. Secondly, mIHC, a multispectral microscopy technique that can reveal the TIME profile on FFPE slides, was applied to determine the major immune cell populations infiltrated in tumor including CD4 + T cells, CD8 + T cells, and macrophages, which were subsequently used to confirm their relationship with PLA2G2D expression. Our results demonstrated that PLA2G2D high-expressed samples had more immune cells infiltrated in the stroma region. More importantly, CD8 + TILs were the only cell population that were more frequently found inside the tumor region in PLA2G2D high-expressed samples, suggesting that PLA2G2D could participate in the recruitment of CD8 + TILs into the cancer nests, either directly or indirectly. The presence of tumor-associated macrophages is often associated with tumorigenesis, immune suppression, and poor prognosis of patients (35, 51) . Indeed, macrophages can directly suppress CD8 + T-cell proliferation in vitro (30, (52) (53) (54) . Therefore, we compared the CD8/CD68 ratio between PLA2G2D-high and PLA2G2D-low samples. Again, a significantly higher ratio of CD8/CD68 was found in PLA2G2D-high samples compared with that in PLA2G2D-low samples. Taken together, these results not only confirmed the aforementioned relationship between immune infiltration and PLA2G2D expression, but also revealed the spatial characteristics of immune infiltration in PLA2G2D-high samples.

The mechanism by which PLA2G2D affects immune infiltration in CSCC is currently unclear, and lipid mediators downstream of PLA2G2D might contribute. For example, the enzyme activity of PLA2G2D could result in the production of leukotriene B4, a potent chemotactic molecule with the capacity to recruit neutrophils, monocytes/macrophages, CD8 + cytotoxic T lymphocyte, and Th17 cells (55) . However, molecules other than these lipid mediators may be involved as well. To explore the possible candidates, genes that co-expressed with PLA2G2D were determined by WGCNA analysis, which resulted in the identification of the top 5 genes, namely, SLAMF6, SLAMF1, SH2D1A, TRAT1, and ZNF831. These genes were successfully validated using another dataset from the CGCI database and confirmed by qRT-PCR using fresh-isolated CSCC samples, with the exception of ZNF831, the expression of which was extremely low and no statistical significance was found between PLA2G2Dhigh and PLA2G2D-low samples. It is noteworthy that these five genes are all associated with immunity. SLAMF1 (CD150) and SLAMF6 belong to the signaling lymphocytic activation molecule (SLAM) family and express on several immune cells including T cells, B cells, and NK cells (56) . Previous studies have shown that the expressions of SLAMF1 in cervical cancer and SLAMF6 in liver cancer were closely related to tumor infiltration of immune cells and patient prognosis (57, 58) . The SH2D1A gene encodes SH2 domain-containing protein 1A, which is often known as SLAM-associated protein. SLAM can transduce the tyrosine signaling pathway to promote interferon-g production in the presence of SH2D1A (56, 59) . TRAT1 encodes the tripartite motif (TRIM) protein which is essential for T-cell activation and positively correlated with the survival of patients with metastatic melanoma (59) . ZNF831 is a transcription factor gene and was also significantly correlated with immune cell infiltration in triple negative breast cancer (60) . Together, in light of the role of the aforementioned co-expressed genes in various immune-related pathways, PLA2G2D may interact with these co-expressed genes indirectly or directly, thereby influencing immune cell infiltration.

Immune cell infiltration, along with PDL1 expression and tumor mutation burden (TMB), has been frequently used as independent biomarkers in predicting the response of patients to ICB therapy in multiple cancer types (61) . Therefore, the positive correlation between immune cell infiltration and the expression of PLA2G2D in CSCC identified in this study indicates that PLA2G2D expression may also represent another predictive biomarker for ICB in CSCC. In support of this notion, we found that the levels of ICP genes, another independent biomarker for ICB therapy, were markedly higher in the PLA2G2D high-expression group than those in the PLA2G2D low-expression group at the transcriptional level. Furthermore, a positive correlation between the expression of PLA2G2D and ICP genes was noted in both TCGA and CGCI cohorts. Finally, we directly examined the predictive value of PLA2G2D expression for patient response to ICB therapy in melanoma and urothelial carcinoma immunotherapy cohorts, given the lack of transcriptional data in CSCC immunotherapy cohorts. Before ICB therapy, PLA2G2D-high patients have higher ICP expression and CYT score. After ICB therapy, the expression level of ICPs and CYT score were further increased compared with the pretreatment stage, as indicated by paired analysis. Hence, PLA2G2D-high patients have higher cytotoxicity and favorable response to ICB treatment. Taken together, PLA2G2D expression may represent a novel biomarker with a better predictive power for ICB therapy.

In conclusion, through integrating bioinformatics and experimental verification, we demonstrated that the metabolic molecule PLA2G2D was positively correlated with immune infiltration and patient prognosis in CSCC, suggesting that PLA2G2D could be a novel prognosis biomarker for CSCC patients. Furthermore, PLA2G2D might be a promising biomarker for the evaluation of immune infiltration situation across different patients and represent another independent predictive biomarker for ICB therapy in CSCC. However, the underpinning mechanism regarding how PLA2G2D works in TME is unclear and warrants further investigations.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

The studies involving human participants were reviewed and approved by the Ethical Committee of Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology (approval number: TJ-IRB2021207). The patients/ participants provided their written informed consent to participate in this study.

HW, YH, HL, and RX conceived and designed the study. HL did bioinformatics analysis and wrote the manuscript. HL and RX did the PCR and mIHC experiments and analyzed the data. CG designed the mIHC panel. HL, CG, TZ, LL, YY, and HZ handled the clinical samples and collected clinical data. HW and YH had supervised this project and contributed to writing and revision of the manuscript. All authors contributed to the article and approved the submitted version. 

Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries

Cervical Cancer

Genome-Wide Profiling of HPV Integration in Cervical Cancer Identifies Clustered Genomic Hot Spots and a Potential Microhomology-Mediated Integration Mechanism

The Need for Improvement of the Treatment of Advanced and Metastatic Cervical Cancer, the Rationale for Combined Chemo-Immunotherapy. Anti-Cancer Agents in Medicinal Chemistry

Improved Survival With Bevacizumab in Advanced Cervical Cancer

New Strategies in Advanced Cervical Cancer: From Angiogenesis Blockade to Immunotherapy

TCR-Engineered T Cells Targeting E7 for Patients With Metastatic HPV-Associated Epithelial Cancers

Safety and Efficacy of Adoptive Cell Transfer Using Autologous Tumor Infiltrating Lymphocytes (LN-145) for Treatment of Recurrent, Metastatic, or Persistent Cervical Carcinoma

Association of Ipilimumab With Safety and Antitumor Activity in Women With Metastatic or Recurrent Human Papillomavirus-Related Cervical Carcinoma

The Tumor Microenvironment in the Response to Immune Checkpoint Blockade Therapies

Efficacy and Safety of Pembrolizumab in Previously Treated Advanced Cervical Cancer: Results From the Phase II KEYNOTE-158 Study

Safety and Efficacy of Pembrolizumab in Advanced, Programmed Death Ligand 1-Positive Cervical Cancer: Results From the Phase Ib

PD-1 Blockade Induces Responses by Inhibiting Adaptive Immune Resistance

Predictive Correlates of Response to the Anti-PD-L1 Antibody MPDL3280A in Cancer Patients

Analysis of Ugandan Cervical Carcinomas Identifies Human Papillomavirus Clade-Specific Epigenome and Transcriptome Landscapes

Tumor and Microenvironment Evolution During Immunotherapy With Nivolumab

Contribution of Systemic and Somatic Factors to Clinical Response and Resistance to PD-L1 Blockade in Urothelial Cancer: An Exploratory Multi-Omic Analysis

Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters

Estimating the Population Abundance of Tissue-Infiltrating Immune and Stromal Cell Populations Using Gene Expression

Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data

Simultaneous Enumeration of Cancer and Immune Cell Types From Bulk Tumor Gene Expression Data

Molecular and Pharmacological Modulators of the Tumor Immune Contexture Revealed by Deconvolution of RNA-Seq Data

Molecular and Genetic Properties of Tumors Associated With Local Immune Cytolytic Activity

WGCNA: An R Package for Weighted Correlation Network Analysis

Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks

Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2

Characteristics of Infiltrating Immune Cells and a Predictive Immune Model for Cervical Cancer

The Prognostic Landscape of Tumor-Infiltrating Immune Cells in Cervical Cancer

Macrophages and Therapeutic Resistance in Cancer

Immune Cells Within the Tumor Microenvironment: Biological Functions and Roles in Cancer Immunotherapy

Tumor Immune Microenvironment Characteristics and Their Prognostic Value in Non-Small-Cell Lung Cancer

Understanding the Tumor Immune Microenvironment (TIME) for Effective Therapy

Tumor Cell-Intrinsic Factors Underlie Heterogeneity of Immune Cell Infiltration and Response to

Macrophages as Regulators of Tumour Immunity and Immunotherapy

0 0 E ) P r o m o t e s S t r o m a l C e l l -M ed i a t e d Immunosuppression via Induction of Interleukin-1 in Melanoma

Tumor Microenvironment: A Metabolic Player That Shapes the Immune Response

Metabolic Modulation of the Tumor Microenvironment Leads to Multiple Checkpoint Inhibition and Immune Cell Infiltration

Glutamine Blockade Induces Divergent Metabolic Programs to Overcome Tumor Immune Evasion

Metabolic Regulation of T Cell Longevity and Function in Tumor Immunotherapy

Increased Gene Expression of Novel Cytosolic and Secretory Phospholipase A(2) Types in Human Airway Epithelial Cells Induced by Tumor Necrosis Factor-Alpha and IFN-Gamma

Receptors for a Growing Family of Secreted Phospholipases A2

Lymphoid Tissue Phospholipase A2 Group IID Resolves Contact Hypersensitivity by Driving Antiinflammatory Lipid Mediators

Distinct Expression Pattern of the Full Set of Secreted Phospholipases A2 in Human Colorectal Adenocarcinomas: Spla2-III as a Biomarker Candidate

Characterization and Differentiation-Dependent Regulation of Secreted Liu et al. Prognostic Biomarker for CSCC Frontiers in Oncology | www

Phospholipases A2 in Human Keratinocytes and in Healthy and Psoriatic Human Skin

Critical Role of Phospholipase A2 Group IID in Age-Related Susceptibility to Severe Acute Respiratory Syndrome-CoV Infection

Secretory Phospholipase A2-IID Is an Effector Molecule of CD4+CD25+ Regulatory T Cells

Dual Roles of Group IID Phospholipase A2 in Inflammation and Cancer

A Novel Risk Model Based on Lipid Metabolism-Associated Genes Predicts Prognosis and Indicates Immune Microenvironment in Breast Cancer

Prognostic Value of Lipid Metabolism-Related Genes in Head and Neck Squamous Cell Carcinoma

Tumour-Associated Macrophages as Treatment Targets in Oncology

Macrophage Expression of Hypoxia-Inducible Factor-1 Alpha Suppresses T-Cell Function and Promotes Tumor Progression

Leukocyte Complexity Predicts Breast Cancer Survival and Functionally Regulates Response to Chemotherapy

Macrophage IL-10 Blocks CD8+ T Cell-Dependent Responses to Chemotherapy by Suppressing IL-12 Expression in Intratumoral Dendritic Cells

Lymphatic Dysfunction, Leukotrienes, and Lymphedema

The SLAM Family of Immune-Cell Receptors

Identification of a Tumor Microenvironment-Related Gene Signature to Improve the Prediction of Cervical Cancer Prognosis

A Five-Gene Signature Based on Stromal/ Immune Scores in the Tumor Microenvironment and Its Clinical Implications for Liver Cancer

Immune Profile and Mitotic Index of Metastatic Melanoma Lesions Enhance Clinical Staging in Predicting Patient Survival

Classification of Triple-Negative Breast Cancers Based on Immunogenomic Profiling

Predictive Biomarkers for Checkpoint Inhibitor-Based Immunotherapy

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.