Three Genes Predict Prognosis in Microenvironment of Ovarian Cancer

Ovarian cancer (OC) is the deadliest gynecological cancer in women. Immune cell infiltration has a critical role in regulating carcinogenesis and prognosis in OC. To identify prognostic genes relevant to the tumor microenvironment in OC, we investigated the association between OC and gene expression profiles. Results obtained with the ESTIMATE R tool showed that immune score and stromal score were correlated with lymphatic invasion, and high immune score predicted a favorable prognosis. A total of 342 common differentially expressed genes were identified according to the two scores; these genes were mainly involved in immune response, extracellular region, and serine-type endopeptidase activity. Three immune-related prognostic genes were selected by univariate and multivariate Cox regression analysis. We further established a prognostic model and validated the prognostic value of three hub genes in different databases; our results showed that this model could accurately predict survival and evaluate prognosis independent of clinical characteristics. Three hub genes have prognostic value in OC. TIMER analysis revealed that the three genes were correlated with different immune cells. Low levels of macrophage infiltration and high levels of CD4+ T cell infiltration were associated with favorable survival outcomes. Arm-level gain of GYPC was correlated with neutrophils and dendritic cells. These findings indicate that CXCR4, GYPC, and MMP12 modulate prognosis via effects on the infiltration of immune cells. Thus, these genes represent potential targets for immune therapy in OC.


INTRODUCTION
Ovarian cancer (OC) is the eighth most prevalent cancer and the deadliest gynecological malignancy. Early diagnosis and treatment are associated with favorable prognosis. However, in more than two-thirds of patients, this cancer is diagnosed at a late stage, resulting in a poor prognosis (Guo et al., 2019). Although surgery, chemotherapy, molecular targeted drugs, and immunotherapy have led to rapid advances, the cure rate of patients with advanced OC remains poor, with a 5-year survival rate of approximately 30% Shimizu et al., 2020). Therefore, the effective biomarkers are needed to improve the prognosis of patients with OC.
Accumulating evidence indicates that the tumor microenvironment (TME) has significant clinical value in predicting prognosis (Pu et al., 2019;Hong et al., 2020;Ling et al., 2020;Qu et al., 2020). Previous studies indicated that T cells were positively correlated with good clinical outcomes in OC (Fridman et al., 2012). Infiltration of B cells was reported to be linked to excellent prognosis in epithelial OCs (Giraldo et al., 2014). In glioblastoma, the TME is strongly linked to gene

Database
Gene expression profiles and clinical information for OC were downloaded from TCGA using the UCSC Xena Browser 3 (Pu et al., 2019;Du et al., 2020;Zhang Z. et al., 2020), and TCGA ovarian carcinoma gene expression data were obtained using the AffyU133a array. Scores were computed with the ESTIMATE algorithm 4 .

Identification of DEGs Based on Different Immune Scores and Stromal Scores
Differential expression analysis of the data was performed using the limma package Xie et al., 2020). DEGs were selected using fold change > 1.5 and adj. P < 0.05. The ClustVis web tool 5 was used to generate heat maps and for clustering (Metsalu and Vilo, 2015). The jvenn software 6 was used to display shared upregulated and downregulated genes between the two score groups (Bardou et al., 2014;Liu et al., 2020).

FUNCTIONAL ENRICHMENT ANALYSIS OF 342 INTERSECTING GENES
Gene ontology (GO) analysis of the 342 common DEGs was performed using DAVID (see "footnote 1") (Huang et al., 2009). GO functional enrichment was determined using the criterion of false discovery rate (FDR) less than 0.05 (Zhang Z. et al., 2020).

Selection of Prognostic DEGs
Univariate Cox regression was used for preliminary identification of prognostic DEGs , followed by multivariate analysis to estimate the prognostic value of the identified DEGs . To construct a risk assessment system, the risk score for each patient was calculated according to regression coefficients and mRNA expression levels . Based on the median risk score, patients were allocated to high-and low-risk groups . Kaplan-Meier (KM) survival curves were constructed to evaluate the prognostic value of the DEGs. To estimate the prognostic value of the risk score model, we plotted 3-and 5-year receiver operating characteristic (ROC) curves and calculated the area under the curve (AUC) using the survival and time R packages. Finally, univariate and multivariate independent factor prognostic analyses were used to assess whether the signature could predict patient prognosis independent of clinical characteristics such as age, grade, stage, lymphatic invasion, and subdivision (Hu et al., 2020).

Validation of the Prognostic Value of Three Hub Genes
PROGgenesV2 7 is a web based tool that allows researchers to study the association between genes and overall survival (OS) in multiple cancers based on TCGA and GEO data (Goswami and Nakshatri, 2014). PrognoScan provides a powerful platform for evaluating the relation between gene expression and patient prognosis such as OS and disease-free survival (DFS) across a large collection of publicly available cancer microarray datasets. The database is publicly accessible at http://www.prognoscan.org/ (Mizuno et al., 2009). Online consensus Survival for Ovarian Cancer (OSov) encompasses 22 expression datasets and provides six types of survival terms for 3212 patients of OC, which can be available at http://bioinfo. henu.edu.cn/OV/OVList.jsp (Zheng et al., 2020). In this study, PROGgenesV2, PrognoScan, and OSov database was used to assess the prognostic value of hub genes. P ≤ 0.05 was considered as significant.

Association Between Prognostic Signature and Immune Cell Infiltration
Tumor Immune Estimation Resource contains seven modules: gene, survival, mutation, SCNA, diff exp, correlation, and estimation (Hu et al., 2020). In this study, we used the gene, survival, and SCNA modules to evaluate the association between prognostic biomarkers and immune cell infiltration.

Statistical Analysis
Univariate analysis and multivariate Cox regression analyses were performed using the R survival package. Risk plots and risk survival curves were analyzed and visualized using the survival and pheatmap packages. All these analyses used R version 3.6.3 8 (Luan et al., 2019). 8 https://www.r-project.org/

Immune Scores and Stromal Scores Are Related to Clinical Factors in OC
We estimated immune score and stromal score based on gene expression data using the ESTIMATE algorithm (Xie et al., 2019;Mahal et al., 2020;Pan et al., 2020). Both scores were correlated with lymphatic invasion (P = 0.0005 and P = 0.0008, respectively). There were no associations between the scores and other clinical factors, including anatomic neoplasm subdivision, grade, sample type, and stage ( Figures 1A-J). KM survival analysis demonstrated that patients with low immune scores had poor prognosis (P = 0.0389) ( Figure 1K). The stromal scores showed no association with OS ( Figure 1L).

Identification of DEGs
A total of 593 OC patients were classified into high and low immune/stromal score groups. To determine the association between gene expression and immune and stromal scores, we compared gene expression profiles between groups. As shown in Figure 2A, 531 DEGs were detected in the high immune score group, including 58 upregulated genes and 473 downregulated genes (fold change > 1.5, P < 0.05), and 500 DEGs were obtained in the high stromal score group, including 25 upregulated and 475 downregulated genes (fold change > 1.5, P < 0.05) ( Figure 2B). Volcano plots (Figures 2C,D) were used to display the distribution of the DEGs in the two score groups. According to the Venn analysis, there were 342 common genes, including 324 downregulated DEGs and 18 upregulated DEGs (Figures 2E,F).

Functional Annotation of 342 Intersecting Genes
Forty-four GO terms in the biological process category were significantly enriched. The top five biological process GO terms were immune response, inflammatory response, interferongamma-mediated signaling pathway, regulation of immune response, and extracellular matrix organization. A total of 19 cellular component GO terms were enriched; these were mainly related to the extracellular region. Eleven GO terms in the molecular function category were significantly overrepresented (FDR < 0.05) (Figure 3).

Identification of Prognostic Signature and Evaluation of Prognostic Model
We first conducted univariate Cox regression analysis to identify prognostic genes among the 342 common DEGs. Our result showed that seven DEGs had prognostic value for OC (Table 1). We next performed multivariate Cox regression analysis on these seven DEGs to identify a prognostic signature in OC and obtained three DEGs that constituted a prognostic signature for OC. Survival curves were plotted according to the expression levels of these three DEGs. CXCR4 and MMP12 were identified as favorable prognostic factors, whereas GYPC has poor predictive value (Figures 4A-C). Moreover, we constructed a risk score model for predicting prognosis, which included risk score ranking, survival status, and a heat map of the three DEGs' expression levels (Figures 4D-F). We further plotted survival curves for the high-and low-risk groups. As shown in Figure 4G, the survival rate of patients in the high-risk groups was lower than the survival rate of those in the low-risk groups (P = 0.00021). Besides, the 3-and 5-year AUC values were 0.701 and 0.727, FIGURE 3 | Function enrichment analysis of 342 intersection genes. Top 10 enriched GO terms of 342 common DEGs. Significantly GO terms were identified using DAVID functional annotation tool (https://david.ncifcrf.gov/). FDR smaller than 0.05 was used as the criteria for significantly enriched. Enriched GO terms included three categories: biological process, cellular component, and molecular function.
respectively, indicating that the risk score model can effectively evaluate prognosis (Figure 4H). Univariate and multivariate independent prognostic analyses were used to assess whether the risk score was independent of other clinical characteristics as an evaluation factor. Our results demonstrated that the risk score could accurately predict survival, regardless of the influence of age, grade, stage, lymphatic invasion, and subdivision; thus, it represents an independent prognostic factor for poor survival in OC (Figures 4I,J).

Three Hub Genes Have Prognostic Value in OC
To further confirm the prognostic value of the identified three hub genes, the association between three hub genes and prognosis was analyzed using different databases, including PROGgenesV2, PrognoScan, and OSov. The results indicated that high GYPC expression might be a candidate for poor prognosis in OC patients, and MMP12 expression was positively correlated with good OS and PFS. Based on TCGA dataset, high expression of CXCR4 predicts a favorable prognosis. These results are consistent with our results. However, in GSE8842, GSE13867, and GSE30161, a positive CXCR4 expression indicated a poor prognosis. In summary, these results suggested that three hub genes may be promising biomarkers that predict the prognosis of all OC patients ( Table 2).
Evaluation of the Associations Among the Three Prognostic Genes, Immune Cells, and Survival We assessed the association between the three prognostic biomarkers and immune cell infiltration using TIMER. The   results showed that the expression of CXCR4 was related to CD4+ T cells, neutrophils, and dendritic cells. GYPC was related to four immune cell types: CD8+ T cells, macrophages, neutrophils, and dendritic cells. MMP12 was primarily related to B cells, CD4+ T cells, neutrophils, and dendritic cells (Figures 5A-C). We further used the SCNA module for comparison of tumor infiltration levels among tumors with different SCNAs for the three prognostic biomarkers. Six immune cell types displayed significant downregulation in OC associated with SCNA of GYPC; arm-level gain of GYPC was correlated with neutrophils and dendritic cell. The SCNA level of MMP12 was not significantly associated with immune infiltration (Figures 5D,E). In addition, we used the survival module of TIMER to examine the relationship between clinical characteristics and immune cell infiltration. Our results showed that infiltration levels of CD4+ T cells and macrophages could predict patient survival in OC. High infiltration levels of macrophages and low infiltration levels of CD4+ T cells were associated with poor survival (Figure 5F).

DISCUSSION
Ovarian cancer is among the most prevalent cancers in women . The TME has a major impact on tumor progression, recurrence, and metastasis, resulting in poor prognosis Huang et al., 2019;Zhang F.-P. et al., 2020). Previous studies have shown that the TME influences initiation and progression of OC and affects antitumor treatment . However, the effects of genes associated with the TME on the prognosis of OC patients are still poorly understood. Therefore, identification of prognostic genes related to the TME is critical to enhance the survival of patients with OC. Lymphatic invasion is associated with poor prognosis in OC patients (Matsuo et al., 2012). In the present study, a comprehensive analysis of DEGs in the TME was performed, and associations among DEGs, prognosis, and immune cells were investigated. High immune and stromal scores were associated with lymphatic invasion (P = 0.0005, P = 0.0008), and patients with high immune scores had longer OS (P = 0.0389). Our data suggest that the clinical outcomes of OC patients could be influenced by their TME type (Figures 1A-K).
We next compared gene expression profiles between the different immune/stromal score groups. A total of 342 common genes were identified, including 18 upregulated genes and 324 downregulated genes. GO term analysis resulted in 342 common genes, mainly involved in the TME. The top five biological process GO terms were immune response, inflammatory response, interferon-gamma-mediated signaling pathway, regulation of immune response, and extracellular matrix organization (Figure 3). Our results suggest that DEGs associated with the TME might contribute to the development of OC by affecting the above mentioned biological process.
We further performed univariate and multivariate Cox regression analysis to identify prognostic TME-related genes. Three downregulated DEGs (CXCR4, MMP12, and GYPC) were identified as a prognostic signature for OC, of which CXCR4 and MMP12 were recognized as favorable prognostic factors. GYPC was regarded as a poor prognostic factor (Figures 4A-C). Three databases were carried out to assess the prognostic value of three hub genes. The results indicated high expression of MMP12 and low GYPC expression predicts favorable prognosis in TCGA and GEO datasets. In TCGA dataset, the expression of CXCR4 is positively related to favorable prognosis. These results are in agreement with our results. However, in the GEO dataset, a positive CXCR4 expression was correlated with a poor prognosis. This different result may be related to different statistical methods and sample composition. CXCR4 is a G-proteincoupled chemokine receptor that is upregulated in many cancers, including thyroid, breast, pancreatic, prostate, and kidney cancers, renal cell carcinoma, and OC (Du et al., 2015;Gong et al., 2020). Previous studies showed that upregulation of CXCR4 was strongly correlated with poor prognosis in renal cell carcinoma (Du et al., 2015). In OC, inhibiting the expression of CXCR4 can prevent cell proliferation and promote cell apoptosis, possibly via the MAPK signaling pathway . CXCR4 affects OC progression through enhancing tumor angiogenesis and suppressing immunity/inhibiting immunity (Gil et al., 2014). A previous meta-analysis showed that high CXCR4 expression was linked to poor prognosis in OC (Liu et al., 2014). Our results are inconsistent with this previous study, indicating that further investigation of the underlying mechanism is required. MMP12 is significantly upregulated in gallbladder cancer, in which it has been identified as a novel prognostic biomarker . A meta-analysis and systematic review indicated that MMP12 rs2276109 might have an impact on OC risk (Zhu and Sun, 2017). The G allele of the MMP12 82A/G polymorphism might promote development and progression of epithelial ovarian carcinoma (Li et al., 2009). No association between MMP12 and prognosis in OC has been reported, nor has the role of GYPC in OC been elucidated. In this study, we first identified MMP12 as a favorable prognostic biomarker and GYPC as a negative prognostic biomarker. Taken together, these genes could represent potential prognostic biomarkers for OC.
Increasing evidence indicates that B cells, CAFs, CD4+ T cells, CD8+ T cells, neutrophils, dendritic cells, macrophages, and other cell types have critical roles in OC . In OC, CXCL14 secreted by CAFs stimulates LINC0009 to promote OC growth and invasion (Curtis et al., 2019;Jiang et al., 2020). Moreover, CAFs may affect autophagy in OC cells through a series of pro-inflammatory cytokines, autophagy-derived substrates, and metabolites, leading to OC progression (Thuwajit et al., 2018;Jiang et al., 2020). In high-grade serous OC, high levels of infiltration of CD8+ T cells are usually associated with favorable outcomes (Goode et al., 2017). CD8+ T cells and CD4+ T cells both have prognostic value in OC (Hamanishi et al., 2011). In OC, tumor-associated macrophages are believed to involve in immune inhibition and have a critical role in tumor cell invasion, angiogenesis, metastasis, and early relapse (Yin et al., 2016;Drakes and Stiff, 2018). A previous study showed that dendritic cells have a crucial role in initiating and modulating immune responses and are associated with patient outcomes. Mature dendritic cells are linked to increased immune infiltration and improve the prognosis of patients with OC (Truxova et al., 2018). OC cells can impair dendritic cells' activation, antigen presentation, and differentiation . Formation of neutrophil extracellular traps is associated with improved OS in advanced-grade OC (Muqaku et al., 2020). A meta-analysis suggested that the neutrophilto-lymphocyte ratio could predict OS and progression-free survival of OC patients (Chen et al., 2018;Yin et al., 2019). The DEGs identified in our study were associated with B cells, CAFs, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, dendritic cells, endothelial cells, natural killer cells, and others, consistent with previous studies (Figures 5A-C). Low macrophage infiltration and high CD4+ T cell infiltration are associated with favorable survival outcomes in OC ( Figure 5F). This study confirmed the role of CD4+ T cells in positive outcomes.
We further explored the association between expression of the three prognostic genes and immune cell infiltration. CXCR4 and MMP12 were positively correlated with infiltration of CD4+ T cells, neutrophils, and dendritic cells. GYPC was related to CD8+ T cells, macrophages, neutrophils, and dendritic cells. SCNA of GYPC was downregulated in various immune cells, and arm-level gain of GYPC was associated with neutrophils and dendritic cells (Figures 5D,E). CXCR4 was related to altered immune cells, which can affect the prognosis of gastric cancer . A previous study showed that CXCR4 was upregulated in naive CD4+ and CD8+ T cells and CD4+ central memory T cells, which have a crucial role in modulating immune dysfunction (Ramonell et al., 2017). Macrophages express matrix metalloproteinases, which can promote vascularization (Drakes and Stiff, 2018). These results suggest that CXCR4 might modulate the prognosis of OC via infiltration of CD4+ T cells. GYPC was identified as a poor prognostic factor, which might be linked to macrophages, neutrophils, and dendritic cells. The prognostic value of MMP12 was correlated with various immune cell types. The mechanisms underlying these associations require further investigation.

CONCLUSION
In conclusion, CXCR4, GYPC, and MMP12 appear to affect prognosis via influencing the infiltration of immune cells; therefore, these genes represent potential targets for immune therapy in OC.

DATA AVAILABILITY STATEMENT
The RNAseq data (level 3) and clinical information of ovarian cancer samples can be found in UCSC Xena (http://xena.ucsc.edu/).

AUTHOR CONTRIBUTIONS
YG designed the study and wrote the main text of the manuscript. YW performed the data analysis. WS reviewed the manuscript. PY performed the language editing. JC and HL prepared the figures and the tables.
All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the Shaanxi Provincial Natural Science Foundation (2017JQ8057).