Identification of an Metabolic Related Risk Signature Predicts Prognosis in Cervical Cancer and Correlates With Immune Infiltration

Introduction

Cervical cancer (CC) is the]aqvskip-2pc top two cause of cancer death affecting women worldwide (Siegel et al., 2020). Due to limited therapeutic options, the 5-year overall survival (OS) of CC patients with recurrent, persistent and metastatic is only about 15% (Liu et al., 2019). As an increasing understanding of tumor microenvironment (TME), immunotherapy has been a promising approach to advanced/metastatic CC (Attademo et al., 2020). How to screen the patients who can benefit from immunotherapy and increase the sensitivity of the chemotherapy, is a hot topic and needs to be further explored.

The reprogramming of energy metabolism which an emerging hallmarks of cancer, is well-established prominent hallmarks of cancer (Hanahan and Weinberg, 2011). Interestingly, growing evidence indicates that the fatty acid (FA) metabolism confers a selective advantage for tumor metastatic process (Nath and Chan, 2016; Pascual et al., 2016). Previous reports have identified that the metabolism of tumor cells could affect the metastasis to improve patient survival (Shang et al., 2018). Moreover, tumor cells and immune cells experience metabolically completion, modulating antitumor immune responses (Li et al., 2019). Lipid droplets (LDs) are regarded as cytoplasmic organelles within tumor cells and immune cells to store lipids, which play an important role in modulating cancer metasatsis (Shang et al., 2020). Thus, Targeting metabolic pathways in terms of anticancer immunotherapy are important directions for the treatment of tumors.

In this study, a novel independent prognostic model based on metabolic-related gene (MRG) signature was identified and constructed by microarray analysis and The Cancer Genome Atlas (TCGA) database to reflect the tumor immune status in CC. Additionally, nomograms associated with MRGs were established. This prognostic model could delineate the metabolic features of CC and further monitor the status of immune infiltration. Surprisingly, P4HA2 which was included in this prognostic risk model, could promote the LDs accumulation in CC cells and exhibited a close relationship with the abundance of different cell types in the TME and immune checkpoint blockade. Our findings manifested that these prognositic MRGs system play crucial roles in patient management and provide potential drug targets to develop effective immunotherapies for CC.

Materials and Methods

Patients and Sample Collection

This study was conducted in accordance with the Declaration of Helsinki. Each participant signed an informed consent form, and the Ethical Committee of the Peking University Third Hospital approved this study. Five CC tissues and six normal cervical tissues were collected for mRNA microarray. Another 15 normal or 15 CC tissues were selected for validation. All freshly frozen CC tissues were collected from the Peking University Third Hospital (Beijing, China). All of the cases were histologically confirmed by pathologists and no patients had received chemotherapy or radiotherapy prior to surgery. The normal controls were collected from women who underwent hysterectomy for non-malignant conditions. All tissue specimens were immediately frozen at −80°C.

Cell Culture

The human CC cell line SiHa and HeLa was purchased from the American Type Culture Collection (Manassas, VA, United States) and cultured in Dulbecco’s Modified Eagle’s medium (Gibco; Thermo Fisher Scientific, Inc.) in a humid atmosphere with 5% CO2 at 37°C. Authenticity of SiHa and HeLa cell lines was verified by using short tandem repeat (STR) genotyping.

Data Analysis of Metabolic-Related mRNA Expression Profiles

The differentially mRNA expression profiles were analyzed by Agilent Feature Extraction software (version 11.0.1.1). A comprehensive MRGs set was extracted from the Metabolic Atlas in the Human Protein Atlas (HPA) database. Differentially expressed metabolic-related mRNAs were screened with p < 0.05 and —fold change| > 1.5.

RNA Extraction and Quantitative Real-Time PCR

Total RNA was extracted using TRIzol reagent (TaKaRa, Japan). Quantitative real-time PCR (qRT-PCR) was performed with SYBR Premix Ex Taq (TaKaRa, Dalian, China) using the 7500 Fast Real-Time PCR system (Applied Biosystems, United States). GAPDH was utilized as an internal standard control. The gene-specific primers were as follows: P4HA2, forward 5′-AGTACCAGGCAATGCTGAGT-3′, reverse 5′-CCTCTTCTGTCTACGGGGTG-3′; PD-L1 forward, 5′-TGGCATTTGCTGAACGCATTT-3′, reverse, 5′-TGCAGCCAGGTCTAATTGTTTT-3′; GAPDH, forward 5′-CCTGTTCGACAGTCAGCCGCAT-3′, reverse 5′-GACTCC GACCTTCACCTTCCCC-3′. Relative RNA expression levels were calculated by the relative quantification method (2^–ΔΔCT).

Cell Transfection

SiRNAs targeting P4HA2 and a negative control were provided by GeneCreate (Wuhan, China). The following oligonucleotide against genes were used: siRNA against P4HA2 (5′-CGAGATACTTTCAAGCATTTA-3′). Transfection was conducted using Lipofectamine RNAiMAX Reagent (Invitrogen, United States) according the manufacturer’s recommended protocol.

Quantification of Lipid Droplets Accumulation

For cell LD quantification, the lipophilic fluorescence dye BODIPY 493/503 (Invitrogen) was employed to monitor LD accumulation in SiHa cells. Cells grown on coverslips were fixed with 4% paraformaldehyde and stained with BODIPY 493/503 (1 μg/ml) for 45 min. Nuclear recognition was based on DAPI staining (1 μg/mL). Fluorescence was analyzed by Olympus FV-1000 confocal microscope (Olympus, Tokyo, Japan).

Transwell Assays

2 × 10⁵ cells in 100 μl serum-free medium were loaded in the upper chamber pre-coated with 50 μl Matrigel (BD Biosciences). Medium containing 10% bovine serum albumin were added to the lower chambers. After incubation, cells which adhered to the lower surface were fixed and stained. Invasiveness was determined by counting cells in five randomly selected visual fields.

Data Collection

After excluding patients with missing expression information, the RNA expression data containing corresponding clinical information of CC, were downloaded from the TCGA database, including 303 CC tissues (Supplementary Table 1)¹.

Identification and Construction of the Metabolic Gene Signature

Genes with P values less than 0.0001 in maximally separated Kaplan–Meier analysis were defined as prognostic genes. Then, a least absolute shrinkage and selection operator (LASSO) regression was performed to picks the optimal number of potential metabolic genes to build a prognostic MRGs model. Finally, the Risk Score of each sample was calculated based on the LASSO-Cox regression co-efficiency through the following formula: $\text{Risk\hspace{0.17em}score\hspace{0.17em}}={\displaystyle {\sum }_{i=1}^{n}Coe{f}_i}{}^{\ast }{X}_i$, Coef_i is the risk coefficient of prognostic MRGs and X_i is the expression value. Then the optimal cutoff value was determined by the “sur_cutpoint” function of the “survminer” R package. It automatically calculated the segmentation point with the minimum P value. Patients were divided into Low- and High-risk groups according to the best cutoff values or median values.

Development and Evaluation of the Nomogram

Plotted nomogram calibration curves based on the proper MRGs determined by univariate and multivariate cox regression analysis through “rms” R package. The forest was used to show the P value, HR and 95% confidence interval (CI) of each variable through “forestplot” R package.

Correlation Analysis Between MRGs Model and Tumor Microenvironment Scores

The ESTIMATE algorithm was employed to calculate the microenvironment scores of tumors (including stromal, immune, and TME scores) based on expression data from TCGA database² in the CC tissues. ESTIMATE algorithm is based on single-sample gene set enrichment analysis. For each sample from TCGA database, the gene expression values were normalized and sorted by rank. The empirical cumulative distribution functions of each gene and the other genes in the signature were calculated based on absolute expression rather than differential expression. Then, the microenvironment score is automatically output (Yoshihara et al., 2013). We compared the differences between the three kinds of scores in the high- and low-MRG risk model groups, and explore their impact on the prognosis of CC considering the best cutoff values possessed by three kinds of scores.

Analysis of the Relationship Between Immune Infiltration and MRGs

CIBERSORT is a versatile computational method for quantifying cell fractions from bulk tissue gene expression profiles (Newman et al., 2015). We estimated the fraction of the 22 subtype immune cells in CC patients from TCGA database through the CIBERSORT algorithm to extend the utility of this MRGs signature³. The differences in immune cell infiltration were also evaluated by Tumor Immune Estimation Resource (TIMER⁴) including CD4+T cells, B cells, CD8+T cells, neutrophils, macrophages, and dendritic cells. The correlations between hub MRGs expression and the abundances of immune infiltrates, immune checkpoint molecules were confirmed by TIMER database.

Bioinformatic Analysis

A part of bioinformatic analysis including heatmap, volcano plot, gene set enrichment analysis (GSEA) and correlation coefficient matrix was performed using the OmicStudio tools at https://www.omicstudio.cn/tool. The pan-cancer analysis of immune infiltrates by xCell (Aran et al., 2017) and the different expression of hub MRGs in CC tissues and normal tissues was drawn at www.aclbi.com. All data of normal tissue samples were obtained from GTEx V8 release version⁵. Genetic alterations analysis, mainly mutations and copy number alterations, were performed on the cBioportal platform⁶. The protein expression of MRGs, as reflected by antibody staining, was investigated using the Pathology Atlas portal in the HPA database (Uhlen et al., 2015, 2017; Thul et al., 2017).

Statistical Analysis

All statistical analyses were performed using the SPSS software package (version 13.0) and Prism 5.0 software (GraphPad Prism, Inc., La Jolla, CA, United States). The Kaplan–Meier (KM) method was used for disease-free survival (DFS) and OS analysis. Two-tailed Student’s t-test was used for the comparison of two independent groups. Multiple statistical packages were used in R software (version 4.0.2, R Core Team, Foundation for Statistical Computing, Vienna, Austria) to download TCGA data and generate forest plots, survival curves and the receiver operating characteristic (ROC) curve. P-value < 0.05 was considered significant.

Results

Identification of Differentially Expressed MRGs

A flowchart was graphed to describe our study more visually in Figure 1. To identify genes involved in metabolism that affect the prognosis of CC, the MRGs sets including 3,616 MRGs and 137 metabolic pathways were firstly analyzed using KM survival analysis. A total of 189 MRGs include in 93 metabolic pathways were defined as prognostic MRGs. Then, we identified the profiles of differentially expressed MRGs between CC and normal cervical tissues (Figure 2A). There were 23 differentially expressed MRGs, containing 17 downregulated and six upregulated genes, were confirmed by consistency analysis between gene expression and prognostic value (Figure 2B). The results of GSEA analysis confirmed that the differential MRGs were related to carbohydrate derivative metabolic process, nucleobase containing small molecule metabolic process, organophosphate metabolic process and lipid metabolic process (Figures 2C–F).

FIGURE 1

Figure 1. Schematic overview of the whole study.

FIGURE 2

Figure 2. Acquisition of differentially expressed MRGs and construction of the metabolic prognostic signature. (A) Clustered heatmaps to explain the different expressions between CC and normal tissues. (B) The volcano plot of the 54 differentially expressed metabolic mRNAs. Red dots represent upregulated mRNAs, and blue dots represent downregulated mRNAs. (C–F) Gene set enrichment analysis (GSEA) results of significant metabolic-associated biological processes. (G) The LASSO was utilized to validate the parameter selection adjustment. (H) The distribution of LASSO coefficient profiles for metabolic model as prognostic factors. Forest plots of hazard ratios of survival-associated MRGs obtained using univariate (I) and multivariate (J) Cox regression analysis.

Identification of the Prognostic MRGs Model

The differentially expressed prognostic MRGs were analyzed in a LASSO-Cox regression model (Figures 2G,H). Finally, five candidate prognostic MRGs signature were identified, namely, Prolyl 4-Hydroxylase Subunit Alpha 2 (P4HA2), Prolyl 4-Hydroxylase Subunit Alpha 1 (P4HA1), ABL Proto-Oncogene 2 (ABL2), Cytochrome P450 Family 4 Subfamily F Member 12 (CYP4F12), and glycolipid transfer protein (GLTP) (Figures 2I,J). Through KM survival analysis, we further identified that high expression of P4HA2, P4HA1 and ABL2 with best cutoff value or median value were unfavorable in CC; however, CYP4F12 and GLTP were favorable prognostic factors for survival (Figure 3A and Supplementary Figure 1A).

FIGURE 3

Figure 3. The survival analysis for 5-MRGs and the nomograms for the prediction of OS. (A) KM curves for OS (top), and DFS (bottom) in patients from TCGA database. (B,C) KM curves of OS and DFS for patients with stage I-IV divided by the best cutoff values. (D,E) KM curves of OS and DFS for patients with stage IIB-IV divided by the best cutoff values. (F,G) KM curves of OS and DFS for stage I-IV or stage IIB-IV patients divided by the median values. (H,I) ROC curve shows the sensitivity and specificity of MRGs model for predicting OS and DFS. (J) The forest plots of the correlation of clinical features and MRGs model with OS. (K) Nomogram for predicting the 1-, 3- or 5-year OS.

Then, we calculated a prognostic risk score with the following formula: Risk Score = [P4HA2 × (0.485)] + [P4HA1 × (0.211)] + [ABL2 × (0.179)] + [CYP4F12 × (−0.167)] + [GLTP × (−0.347)] (Table 1). The patients were divided into high-risk and low-risk groups according to the best cutoff values. The KM analysis showed that patients (stage I–IV) in the high-risk group had worse OS and DFS than patients in the low-risk group (Figures 3B,C and Supplementary Figure 1B). All patients with advanced-stage (IIB-IV) were further evaluated by the risk score. Similar significant predictive power of the 5-MRGs signature was found in both OS (n = 109) and DFS (n = 87) (Figures 3D,E and Supplementary Figure 1C). At the same time, OS was also significantly worse in the high-risk group with the median cutoff value (Figures 3F,G). Significantly, the overall 5-year survival rate was 47.548% /19.017 in high risk group vs. 81.400%/78.739% in low risk group in all stage patients; for advanced-stage (IIB-IV) patients, the overall 5-year survival rate was 35.997%/7.856% in high risk group vs. 76.977%/67.389% in low risk group (Table 1). Additionally, the area under the ROC curve (AUC) value of 0.783/0.737 indicated great specificity and sensitivity of the MRGs signature in predicting OS for patients with CC (Figures 3H,I and Supplementary Figure 1D). As shown in Table 1, the negative predictive values of MRGs model for DFS was 88% and the negative predictive values for OS was 59%. Taken together, these data may provide a powerful candidate prognostic MRGs signature model to serve as a promising biomarker for predicting outcomes in CC patients.

TABLE 1

Table 1. Clinical significance of the 5-metabolic gene signature-based risk score in predicting prognosis of patients in cervical cancer.

Independent Prognostic Validation of MRGs Model and Construction of Prognostic Nomograms

To explore whether MRGs model as an independent OS predictor, offers a measurable method to estimate the prognosis of CC patients, we constructed a nomogram that contained the MRGs model and other clinical characteristics, including age, nodal metastasis, grade, and stage. The multivariate Cox regression showed high MRGs risk score (HR 2.50; 95% CI, 1.49–4.2) and nodal metastasis (HR 3.02; 95% CI, 1.49–6.1) were independent prognostic factors (Figure 3J). A nomogram including MRGs model and nodal metastasis could provide a clinically applicable method for predicting the prognosis of CC patients. The C-index was 0.71 (95% CI: 0.617–0.803) for the TCGA cohorts (Figure 3K). Thus, our model had an approximately moderate accuracy in the prognostic prediction of CC to some extent.

Association Between the MRGs Model and TME

To further understand the relationship between metabolic status and TME in CC, immune score and stromal score were enrolled into the analysis. Based on the ESTIMATE method, we obtained stromal scores (range: −2586.99∼778.01) and immune scores (range: −712.91∼141.82) for all these CC patients. Both immune score and TME score of CC cases were significantly lower in high-risk group than those of low-risk cases (P = 0.0003, 0.0174, respectively), though no significant differences between high- and low-risk group were found for stromal scores (Figures 4A–C). The MRGs model were significantly negatively correlated with the immune scores (Figure 4D).

FIGURE 4

Figure 4. The clinical value of ESTIMATE score. (A–C) The relationship between immune, stromal, TME score and MRGs model. (D) The MRGs model were significantly negatively correlated with the immune scores. (E–G) The predictive value of stromal, immune and TME score in evaluating OS. (H–J) The predictive value of stromal, immune and TME score in evaluating DFS. (K–M) The immune landscape in clinicopathologic features: metastasis, lymph node metastasis, and margins involvement.

The predictive value of stromal/immune score in evaluating OS and DFS was explored based on TCGA database. As shown in Figures 4E–G and Supplementary Figure 1E, the OS of patients with high immune score or TME score was longer than those in low score group; however, the patients with low stromal score only showed a trend toward better OS. Consistently, patients with low immune or TME score and high stromal score was found to have significantly worse DFS (Figures 4H–J and Supplementary Figure 1F).

The Immune Landscape in Clinicopathologic Features, Including MRGs Model

The tumor infiltrating immune cells were variously distributed in different clinicopathologic characteristics. We found that the infiltrating levels of memory activated CD4+T cells were significantly lower and regulatory Tregs cells were higher in primary tumor microenvironment (PTME) with metastasis (Figure 4K). It’s worth noting that the infiltrating levels of resting mast cells were significantly higher and CD8+T cells, memory activated CD4+T cells, resting NK cells were lower in PTME with lymph node metastasis (Figure 4L). The infiltrating levels of neutrophils cells displayed downward trend and macrophage (M0) cells, activated mast cells exhibited rising trend in the PTME with margins involvement (Figure 4M, p > 0.05).

Furthermore, we calculated the relationship between MRGs model and immune infiltration through CIBERSORT system. A total of six subtypes of immune cells [memory B cells, follicular helper T cell, regulatory T cells, macrophage (M0) cells, resting mast cells and activated mast cells] had an obvious negative correlation with MRGs model (Figure 5A). Meanwhile, there was a significant positive correlation between naïve B cells and plasma cells (Figure 5B).

FIGURE 5

Figure 5. The correlation between MRGs model and immune cells. (A) The relationships between MRGs model and infiltration levels of 22 subtypes of immune cells calculated by CIBERSORT system. (B) Correlation heatmap analysis of immune cells.

Validation of MRGs Expression

All of the five hub prognostic MRGs were validated in protein and mRNA levels using TCGA and The HPA database. We found that P4HA2, P4HA1, and ABL2 had the higher expression levels in CC tissues than that in normal cervical tissues, while CYP4F12 and GLTP had the lower expression levels by immunohistochemistry (IHC) analysis (Figure 6A). Furthermore, the RNA-seq results showed that P4HA2, P4HA1, and CYP4F12 had the consistent expression trend, while the mRNA expression of ABL2 and GLTP in clinical sample tissues is just the opposite (Figure 6B). However, we did not observe significant expression differences of MRGs in patients with pathological staging, grade, and nodal metastasis status (Supplementary Figure 1G).

FIGURE 6

Figure 6. Validation of expression and pan-analysis of immune infiltration. (A) The representative protein expression of the 5-MRGs in CC tissues and normal tissues. (B) The mRNA expression levels of 5-MRGs. (C) Genetic alterations of the 5-MRGs in CC. (D) The differential expression between tumor and normal tissues for P4HA2 across all TCGA tumors. (E) Correlation heatmap between P4HA2 and immune infiltration across all TCGA tumors. (F) Correlation heatmap between P4HA2 and immune checkpoint genes across all TCGA tumors. *P < 0.05, **P < 0.01, ***P < 0.001. CNS, central nervous system; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; CESC, cervical and endocervical cancer; DLBC, diffuse large B-cell lymphoma; GBM, glioblastoma multiforme; HNSC-HPVpos, head and neck cancer-HPV positive; HNSC-HPVneg, head and neck cancer-HPV negative; KICH, kidney chromophobe; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; TGCT, testicular germ cell tumors; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma.

Genetic alteration are considered to be a key factor in immune tolerance breakdown (Small et al., 2018). Thus, the genetic alteration in the MRGs model was analyzed with cBioPortal software. Among them, ABL2, CYP4F12, and P4HA2 were altered in six (2%) and four (1.3%) from the 308 patients, respectively (Figure 6C). The amplification and deep deletion were the most common genetic alteration forms (Figure 6C).

Pan-Cancer Analysis of P4HA2 Expression

P4HA2 showed consistent expression patterns at mRNA and protein levels and was an independent risk factor for evaluating patient prognosis based on the foregoing results (Figures 2J, 6B). Thus, we selected P4HA2 to compare its RNA sequencing data with that in corresponding normal tissues across cancers in TCGA using TIMER (Figure 6D). Consistent with the expression trend in CC, P4HA2 expression was upregulated in BRCA (breast invasive carcinoma), CHOL (cholangiocarcinoma), ESCA (esophageal carcinoma), HNSC (head and neck cancer), KIRC (kidney renal clear cell carcinoma), LIHC (liver hepatocellular carcinoma), PCPG (pheochromocytoma and paraganglioma), STAD (stomach adenocarcinoma), and THCA (thyroid carcinoma). The opposite trend was only observed in PRAD (prostate adenocarcinoma). Furthermore, the higher expression of P4HA2 was found in HNSC patients with HPV negative than those with HPV positive. Collectively, P4HA2 acts as a tumor promoter in most cancer types.

Association Between P4HA2 Expression and Immune Infiltrates or Immune Checkpoint Across Cancers

It has long been recognized that immune cells are intimately associated with tumor cells. We first investigated the correlation between infiltration of immune cell subtypes and P4HA2 expression is based on the xCell across 33 cancers (Figure 6E). The results indicated that P4HA2 was notably correlated with CD8+T cells in 11 cancer types, B cells and memory B cells in 12 cancer types, plasmacytoid dendritic cells and naïve CD4+ T cells in 13 cancer types, and CD4+Th1 cells in 17 cancer types (P < 0.0001). Three top cancer types, BRCA, LGG (lower grade glioma) and THCA, exhibited most strongly correlation between P4HA2 with immune cells (P < 0.0001). Moreover, P4HA2 expression was significantly correlated with the immune score in 12 cancer types, mircroenvironment score in 14 cancer types, and stromal score in four cancer types (P < 0.0001). Consistent with the result in Figures 4B,D, the P4HA2 expression are significantly negatively correlated with the immune score and other immune cell subtypes in CC.

To clarify the potential mechanisms underlying the involvement of P4HA2 in the enhancement of immune cell infiltration, we calculated the correlations of P4HA2 with immune checkpoints (CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, SIGLEC15, and TIGIT) (Figure 6F). All of the eight immune checkpoints were significantly correlated with the P4HA2 expression in BRCA, COAD (colon adenocarcinoma), LGG, and PRAD. Furthermore, P4HA2 expression was positively correlated with immune checkpoints in the majority of cancers. Therefore, P4HA2 was not only correlated with the extent of immune infiltration, but potentially played a key role in regulating immune evasion.

Correlation of P4HA2 Expression With Immune Infiltration

To further explore the potential value of combined immunotherapy and targeting P4HA2 for the synergistic treatment of CC, we obtained the coefficient of P4HA2 and immune cells using TIMER database. No correlations was found between P4HA2 and tumor purity (Figure 7A). The results also revealed that P4HA2 was negatively correlated with CD8+T cells, naïve CD8+T cells, central memory CD8+T cells, effector memory CD8+T cells, naïve CD4+T cells, CD4+Th1 cells, B cells, memory B cells, macrophage (M0, M1, and M2), (activated) myeloid dendritic cells, plasmacytoid dendritic cells, activated mast cells, and follicular helper T cells (Figures 7B–D,F–I). On the contrary, P4HA2 expression was positively correlated with CD4+T cells, memory resting CD4+T cells, neutrophil cells, resting mast cells, MDSC cells, and NK cells (Figures 7C,E,H,J,K). Further information was available in the Supplementary Figures 2A–G. Moreover, P4HA2 positively correlated with endothelia cells (r = 0.299, P = 3.87e-07) and cancer associated fibroblast (CAF) cells (r = 0.35, P = 2.10e-09) (Figures 7L,M).

FIGURE 7

Figure 7. P4HA2 were noticeably correlated with immune infiltration. (A–K) The relationships between P4HA2 and immune cell infiltration levels calculated by TIMER database. (L,M) The relationships between endothelial cells, cancer associated fibroblast cells and P4HA2 expression. (N,O) The relationships between TMB, MSI, and P4HA2. (P) The differential expression between POLE/STK11 wild-type and POLE/STK11 mutation-type group for P4HA2.

P4HA2 Is Associated With the Microsatellite Instability and Immune Checkpoint

Microsatellite instability (MSI), tumor mutational burden (TMB) and PDL-1 are the most important biomarkers for immune therapy in clinical practice (Mazloom et al., 2020). According to correlation analysis, P4HA2 was significantly positively correlated with MSI but not associated with TMB (Figures 7N,O). In addition, patients with POLE/STK11 mutations are more sensitive to immune checkpoint inhibitors (Skoulidis et al., 2018; Wang et al., 2019). P4HA2 was upregulated in POLE/STK11 mutation groups of CC patients (Figure 7P). Next, we analyzed the correlation between P4HA2 and immune checkpoint genes. P4HA2 was positively correlated with CD274 (PD-L1) and negatively correlated with CTLA4, LAG3, and PDCD1 (Figures 8A–D). However, P4HA2 was not significantly correlated with HAVCR2, PDCD1LG2, TIGIT, and SIGLEC15 (Supplementary Figure 2H).

FIGURE 8

Figure 8. P4HA2 might promote CC metastasis through regulating the LDs accumulation and immune escape. (A–D) The relationships between P4HA2 and immune checkpoint genes. (E) The qRT-PCR results showed that P4HA2 was more highly expressed at the CC tissues than normal uterine cervical tissues. (F) LDs were detected by double staining with BODIPY 493/503 dye and Hoechst in the indicated cells. Scale bars, 25 um. (G) P4HA2 inhibited SiHa cells invasion by transwell assay. (H) The qRT-PCR analysis of PD-L1 expression in SiHa cells transfected with siP4HA2. **P < 0.01, ***P < 0.001.

Knockdown of P4HA2 Might Inhibit Tumor Metastasis Through Regulating the Lipid Droplet Accumulation and Immune Escape

To further confirm the effects of P4HA2 in regulating immune microenvironment of CC, we firstly detected the expression levels of P4HA2 in CC tissues and normal tissues. Our data showed that P4HA2 was significantly higher in CC tissues compared with normal tissues (Figure 8E). Next, the transfection efficiencies of interference reagents for P4HA2 were confirmed by qRT-PCR (Supplementary Figure 2I). Knockdown of P4HA2 strongly decreased the number of LDs in CC cells, suggesting that P4HA2 modulated LD accumulation in vitro (Figure 8F). Transwell assays indicated that P4HA2 knockdown significantly inhibited the invasion ability of SiHa and HeLa cells (Figure 8G). Importantly, significant downregulation of PD-L1 was observed in CC cells after P4HA2 knockdown (Figure 8H). Thus, we reasoned that P4HA2 might inhibit CC metastasis through regulating tumor immune escape.

Discussion

The appropriate strategy for treating CC depends on tumor stage at diagnosis, but the prognosis of patients with locally advanced or metastatic disease remains dismal (Pfaendler and Tewari, 2016). Immunotherapy as the new frontier of the anticancer treatment, is being an available option in treating this part of CC patients (Ventriglia et al., 2017). Here, for the first time, we identify the metabolic-related prognostic signature to provide an effective prognostic model of CC patients and potential predictive biomarkers of benefit from immunotherapy.

There is growing evidence that metabolism reprogramming has considered as a hallmark of cancer (Ward and Thompson, 2012). It is worth nothing that the Warburg effect may contribute to the early stages of tumor development. Instead, the lipid metabolism confers a selective advantage for tumor metastasis (Nath and Chan, 2016; Pascual et al., 2016). Our previous studies had been also confirmed that non-coding RNA mediated FA metabolism reprogramming could promote the process of lymph node metastasis (Shang et al., 2018). To establish an effective and simple scheme to explore the metabolic status and evaluate clinical outcomes of CC patients, we identify a prognostic MRGs model and patients with high risk scores exhibit a poor prognosis. Moreover, a nomograph composed of MRGs model and other clinical factors was constructed to better predict the OS. Our results illustrated that MRGs model could be applied as independent prognostic factors and indexes of immune infiltration status. Other study have found that immune score can also be used to evaluate the prognosis of CC patients (Chen et al., 2021). However, a larger sample size is needed to evaluate whether MRGs scores can substitute immune scores.

Oncogene-driven metabolism reprogramming can profoundly influence the TME to limit immune responses and put up more barriers to cancer therapy (Bader et al., 2020). And that’s why we focus on the potential of targeting metabolic pathways to improve the TME for favoring cancer immunotherapy. LDs as cytoplasmic lipid-rich organelles, have been verified to play an indispensable role in tumor pre-metastatic microenvironment (Shang et al., 2020). In our research, the MRGs were mainly enriched on the lipid metabolism pathway. Based on the MRGs model, the risk scores were inversely related to the immune score and TME scores. Meanwhile, memory B cells, follicular helper T cells, regulatory Tregs cells and resting mast cells were all negatively associated with the risk scores. Previous studies also showed that inhibition of glucose metabolism results in a drastic reduction of the frequency and number of follicular helper T cells (Choi et al., 2018). Thus, interventions targeting the metabolic pathways might reshape the immune state of the TME.

Collagen prolyl 4-hydroxylases (P4Hs) are located within the lumen of the endoplasmic reticulum and catalyzes post-translational formation of 4-hydroxyproline in -Xaa-Pro-Gly-sequences in collagens (Kivirikko and Pihlajaniemi, 1998). P4HA2 as one of the subtypes of collagen prolyl-4-hydroxylases α isoforms, was activated by HIF-1 to inducing extracellular matrix remodeling under hypoxic conditions and promoting cancer metastasis (Gilkes et al., 2013; Xiong et al., 2014). Our research indicated that P4HA2 was a MRG with independent prognostic value. Significant missense mutation in P4HA2 was observed, and different protein or mRNA expression levels were verified between the CC tissues and normal tissues. Consistent with the study of Cao et al. (2020), P4HA2 was markedly upregulated in CC tissues and functions as an oncogene in promoting cell metastasis by inducing epithelial-mesenchymal transition.

With the help of co-expression analysis, P4HA2 was predicted to be highly connected with multiple immune infiltration cells, including CD8+T cells, CD4+ T cells, B cells, macrophage cells (M0, M1, and M2), dendritic cells, mast cells, neutrophil cells, and NK cells. Moreover, CD8+T cells, memory activated CD4+T cells, resting NK cells and resting mast cells exhibited a significantly difference in evaluating lymph node metastasis. Our results were in accordance with previous research showing that infiltration with CD8+T cells was negatively associated with pelvic lymph node metastasis and predicted poor survival outcomes in CC patients (Ohno et al., 2020). CD8+T cells are preferred immune cells for targeting cancer. PD-1 and CTLA-4 can be targeted for relieving CD8+T cells exhaustion and thereby eliminating antigen-expressing cancer cells (Farhood et al., 2019). Approximately 96% patients with locally advanced CC expressed PD-L1 and tended to have a worse progression-free survival (Enwere et al., 2017). In this study, P4HA2 was negatively associated with CD8+T cells and positively correlated with CD274 (PD-L1) and CAF cells. As previous research illustrate, CAFs could specifically excluding CD8+ T cells from tumors to inhibit immune response and promoted immune evasion (Ford et al., 2020). Thus, P4HA2 might be a highly reliable predictive biomarkers to facilitate patient selection for immune checkpoint inhibitor-based therapies.

There were still several limitations to our research. First, the verification with another independent database was lacked. Second, although the LDs accumulation modulated by P4HA2 had been validate though BODIPY 493/503 staining, the further experimental exploration was still needed to confirm the potential mechanism and clinical utility. At the last, our work only reflect certain aspects of the immune infiltration from the perspective of metabolism reprogramming.

Conclusion

In conclusion, we identified 5-MRGs model and a prognostic nomogram to predict survival of CC patients, which could reflect the immune infiltration status. Our nomogram including MRGs model may provide a reference tool for clinicians to guide follow-ups for CC patients. The hub gene-P4HA2 could modulate the LDs accumulation in CC cells and was closely associated with tumor infiltrating cells lymphocytes and immune checkpoint genes. This study provides a new and feasible method for evaluating the prognosis of CC and explores the potential value to navigate metabolic pathways to enhance anti-tumor immunity and immunotherapy.

Data Availability Statement

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 below: GEO accession GSE173097, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173097.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethical Committee of the Peking University Third Hospital approved this study. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

HG worked on design and conception of this study. JH processed the microarray data and revised the manuscript. CS collected the data, performed the data analysis, and drafted the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by The National Natural Science Foundation of China (81902648) and The China Postdoctoral Science Foundation (2021T140028 and 2018 M641115).

Conflict of Interest

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.677831/full#supplementary-material

Abbreviations

CC, cervical cancer; OS, overall survival; FA, fatty acid; TCGA, The Cancer Genome Atlas; MRG, metabolic-related gene; STR, short tandem repeat; qRT-PCR, quantitative real-time PCR; LASSO, least absolute shrinkage and selection operator; CI, confidence interval; TME, tumor microenvironment; TIMER, Tumor Immune Estimation Resource; HPA, Human Protein Atlas; KM, Kaplan–Meier analysis; GSEA, gene set enrichment analysis; DFS, disease-free survival; ROC, receiver operating characteristic; P4HA2, prolyl 4-hydroxylase subunit alpha 2; P4HA1, prolyl 4-hydroxylase subunit alpha 1; ABL2, ABL proto-oncogene 2; CYP4F12, cytochrome P450 family 4 subfamily F member 12; GLTP, glycolipid transfer protein; AUC, area under the curve; PTME, primary tumor microenvironment; IHC, immunohistochemistry; CAF, cancer associated fibroblast cells; MSI, microsatellite instability; TMB, tumor mutational burden; P4Has, collagen prolyl 4-hydroxylases; LDs, lipid droplets.

Footnotes

1. ^ https://portal.gdc.cancer.gov/
2. ^ https://bioinformatics.mdanderson.org/estimate/index.html
3. ^ https://gdc.cancer.gov/about-data/publications/panimmune
4. ^ https://cistrome.shinyapps.io/timer/
5. ^ https://gtexportal.org/home/datasets
6. ^ www.cbioportal.org

References

Aran, D., Hu, Z., and Butte, A. J. (2017). XCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Attademo, L., Tuninetti, V., Pisano, C., Cecere, S. C., Di Napoli, M., Tambaro, R., et al. (2020). Immunotherapy in cervix cancer. Cancer Treat. Rev. 90:102088. doi: 10.1016/j.ctrv.2020.102088

PubMed Abstract | CrossRef Full Text | Google Scholar

Bader, J. E., Voss, K., and Rathmell, J. C. (2020). Targeting metabolism to improve the tumor microenvironment for cancer immunotherapy. Mol. Cell 78, 1019–1033. doi: 10.1016/j.molcel.2020.05.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Han, Q., Li, J., Jia, Y., Zhang, R., and Shi, H. (2020). P4HA2 contributes to cervical cancer progression via inducing epithelial-mesenchymal transition. J. Cancer 11, 2788–2799. doi: 10.7150/jca.38401

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Lin, H., Pi, Y., Chen, X., Zhou, H., Tian, Y., et al. (2021). Development and validation of a prognostic signature based on immune genes in cervical cancer. Front. Oncol. 11:616530. doi: 10.3389/fonc.2021.616530

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, S. C., Titov, A. A., Abboud, G., Seay, H. R., Brusko, T. M., Roopenian, D. C., et al. (2018). Inhibition of glucose metabolism selectively targets autoreactive follicular helper T cells. Nat. Commun. 9:4369. doi: 10.1038/s41467-018-06686-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Enwere, E. K., Kornaga, E. N., Dean, M., Koulis, T. A., Phan, T., Kalantarian, M., et al. (2017). Expression of PD-L1 and presence of CD8-positive T cells in pre-treatment specimens of locally advanced cervical cancer. Mod. Pathol. 30, 577–586. doi: 10.1038/modpathol.2016.221

PubMed Abstract | CrossRef Full Text | Google Scholar

Farhood, B., Najafi, M., and Mortezaee, K. (2019). CD8+ cytotoxic T lymphocytes in cancer immunotherapy: a review. J. Cell. Physiol. 234, 8509–8521. doi: 10.1002/jcp.27782

PubMed Abstract | CrossRef Full Text | Google Scholar

Ford, K., Hanley, C. J., Mellone, M., Szyndralewiez, C., Heitz, F., Wiesel, P., et al. (2020). NOX4 inhibition potentiates immunotherapy by overcoming Cancer-Associated Fibroblast-Mediated CD8 t-cell exclusion from tumors. Cancer Res. 80, 1846–1860. doi: 10.1158/0008-5472.CAN-19-3158

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilkes, D. M., Bajpai, S., Chaturvedi, P., Wirtz, D., and Semenza, G. L. (2013). Hypoxia-inducible factor 1 (HIF-1) promotes extracellular matrix remodeling under hypoxic conditions by inducing P4HA1, P4HA2, and PLOD2 expression in fibroblasts. J. Biol. Chem. 288, 10819–10829. doi: 10.1074/jbc.M112.442939

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of cancer: the next generation. Cell 144, 646–674. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Kivirikko, K. I., and Pihlajaniemi, T. (1998). Collagen hydroxylases and the protein disulfide isomerase subunit of prolyl 4-hydroxylases. Adv. Enzymol. Relat. Areas Mol. Biol. 72, 325–398. doi: 10.1002/9780470123188.ch9

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Wenes, M., Romero, P., Huang, S. C., Fendt, S., and Ho, P. (2019). Navigating metabolic pathways to enhance antitumour immunity and immunotherapy. Nat. Rev. Clin. Oncol. 16, 425–441. doi: 10.1038/s41571-019-0203-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Wu, L., Tong, R., Yang, F., Yin, L., Li, M., et al. (2019). PD-1/PD-L1 inhibitors in cervical cancer. Front. Pharmacol. 10:65. doi: 10.3389/fphar.2019.00065

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazloom, A., Ghalehsari, N., Gazivoda, V., Nimkar, N., Paul, S., Gregos, P., et al. (2020). Role of immune checkpoint inhibitors in gastrointestinal malignancies. J. Clin. Med. 9:2533. doi: 10.3390/jcm9082533

PubMed Abstract | CrossRef Full Text | Google Scholar

Nath, A., and Chan, C. (2016). Genetic alterations in fatty acid transport and metabolism genes are associated with metastatic progression and poor prognosis of human cancers. Sci. Rep. U. K. 6:18669. doi: 10.1038/srep18669

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohno, A., Iwata, T., Katoh, Y., Taniguchi, S., Tanaka, K., Nishio, H., et al. (2020). Tumor-infiltrating lymphocytes predict survival outcomes in patients with cervical cancer treated with concurrent chemoradiotherapy. Gynecol. Oncol. 159, 329–334. doi: 10.1016/j.ygyno.2020.07.106

PubMed Abstract | CrossRef Full Text | Google Scholar

Pascual, G., Avgustinova, A., Mejetta, S., Martín, M., Castellanos, A., Attolini, C. S., et al. (2016). Targeting metastasis-initiating cells through the fatty acid receptor CD36. Nature 541, 41–45. doi: 10.1038/nature20791

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfaendler, K. S., and Tewari, K. S. (2016). Changing paradigms in the systemic treatment of advanced cervical cancer. Am. J. Obstet. Gynecol. 214, 22–30. doi: 10.1016/j.ajog.2015.07.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer statistics, 2020. CA Cancer J. Clin. 70, 7–30. doi: 10.3322/caac.21590

PubMed Abstract | CrossRef Full Text | Google Scholar

Shang, C., Qiao, J., and Guo, H. (2020). The dynamic behavior of lipid droplets in the pre-metastatic niche. Cell Death Dis. 11:990. doi: 10.1038/s41419-020-03207-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Shang, C., Wang, W., Liao, Y., Chen, Y., Liu, T., Du, Q., et al. (2018). LNMICC promotes nodal metastasis of cervical cancer by reprogramming fatty acid metabolism. Cancer Res. 78, 877–890. doi: 10.1158/0008-5472.CAN-17-2356

PubMed Abstract | CrossRef Full Text | Google Scholar

Skoulidis, F., Goldberg, M. E., Greenawalt, D. M., Hellmann, M. D., Awad, M. M., Gainor, J. F., et al. (2018). STK11/LKB1 mutations and PD-1 inhibitor resistance in KRAS-Mutant lung adenocarcinoma. Cancer Discov. 8, 822–835. doi: 10.1158/2159-8290.CD-18-0099

PubMed Abstract | CrossRef Full Text | Google Scholar

Small, M., Treilleux, I., Couillault, C., Pissaloux, D., Picard, G., Paindavoine, S., et al. (2018). Genetic alterations and tumor immune attack in Yo paraneoplastic cerebellar degeneration. Acta Neuropathol. 135, 569–579. doi: 10.1007/s00401-017-1802-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Thul, P. J., Akesson, L., Wiking, M., Mahdessian, D., Geladaki, A., Ait, B. H., et al. (2017). A subcellular map of the human proteome. Science. 356:eaal3321. doi: 10.1126/science.aal3321

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlen, M., Fagerberg, L., Hallstrom, B. M., Lindskog, C., Oksvold, P., Mardinoglu, A., et al. (2015). Proteomics. Tissue-based map of the human proteome. Science 347:1260419. doi: 10.1126/science.1260419

PubMed Abstract | CrossRef Full Text | Google Scholar

Uhlen, M., Zhang, C., Lee, S., Sjostedt, E., Fagerberg, L., Bidkhori, G., et al. (2017). A pathology atlas of the human cancer transcriptome. Science 357:eaan2507. doi: 10.1126/science.aan2507

PubMed Abstract | CrossRef Full Text | Google Scholar

Ventriglia, J., Paciolla, I., Pisano, C., Cecere, S. C., Di Napoli, M., Tambaro, R., et al. (2017). Immunotherapy in ovarian, endometrial and cervical cancer: state of the art and future perspectives. Cancer Treat. Rev. 59, 109–116. doi: 10.1016/j.ctrv.2017.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Zhao, Q., Wang, Y. N., Jin, Y., He, M. M., Liu, Z. X., et al. (2019). Evaluation of POLE and POLD1 mutations as biomarkers for immunotherapy outcomes across multiple cancer types. JAMA Oncol. 5, 1504–1506. doi: 10.1001/jamaoncol.2019.2963

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, P. S., and Thompson, C. B. (2012). Metabolic reprogramming: a cancer hallmark even warburg did not anticipate. Cancer Cell 21, 297–308. doi: 10.1016/j.ccr.2012.02.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, G., Deng, L., Zhu, J., Rychahou, P. G., and Xu, R. (2014). Prolyl-4-hydroxylase alpha subunit 2 promotes breast cancer progression and metastasis by regulating collagen deposition. BMC Cancer 14:1. doi: 10.1186/1471-2407-14-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar