Microenvironment and related genes predict outcomes of patients with cervical cancer: evidence from TCGA and bioinformatic analysis

Cervical cancer (CESC) is one of the most common cancers and affects the female genital tract. Consistent HPV infection status has been determined to be a vital cause of tumorigenesis. HPV infection may induce changes to the immune system and limit the host’s immune response. Immunotherapy is therefore essential to improving the overall survival of both locally advanced and recurrent CESC patients. Using 304 relevant samples from TCGA, we assessed immune cell function in CESC patients to better understand the status of both tumor micro-environment cells and immune cells in CESC. Functional enrichment analysis, pathway enrichment analysis, and PPI network construction were performed to explore the differentially expressed genes (DEGs). The analysis identified 425 DEGs, which included 295 up-regulated genes and 130 down-regulated genes. We established that upregulation of CCL5 was correlated with significantly better survival, meaning that CCL5 expression could serve as a novel prognostic biomarker for CESC patients. We further focused on CCL5 as a hub gene in CESC, as it had significant correlations with increased numbers of several types of immune cells. Cell-type fractions of M1 macrophages were significantly higher in the high-immune-scores group, which was associated with better overall survival. Finally, we concluded that CCL5 is a promising prognostic biomarker for CESC, as well as a novel chemotherapeutic target.


Introduction
Cervical cancer, including cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), is one of the most common cancers. Cervical cancer affects the female genital tract, and its onset is closely related to infection by human papillomavirus (HPV) (He et al., 2016). Currently, the yearly incidence of CESC is dramatically increasing, with a trend towards younger patient age.
As has been well-established, HPV is a vital factor in causing CESC (Sichero et al., 2019). Thanks to the widespread use of HPV vaccinations and extensive cervical screenings, more than 90% of CESC patients are curable when diagnosed during Stages I-II. However, surgical options for CESC patients who have reached an advanced stage by the time of diagnosis may not be available. In advanced cases, radiotherapy and chemotherapy may be used as an optional treatment in addition to surgery. However, some CESC patients with poor sensitivity to conventional radio-chemotherapy had worse overall survival outcomes with this approach.
Immunotherapeutic interventions are essential for advanced-stage patients in order to improve their overall survival (De Felice et al., 2018;Dyer et al., 2019). Immunotherapy is emerging as an available option to treat a variety of malignancies, such as melanoma and renal cancers, and interest in the field has been significantly increasing over the past few decades (Ramanathan et al., 2018). It has been reported that, because it can promote evasion from host T-cell immune-surveillance, HPV infection stimulates the progression of carcinogenesis from low-grade pathological changes to CESC (Chauhan and Bharadwaj, 2018;Whiteside, 2008). HPV infection could induce modifications to the immune system and limit the host's immune response. These changes could include compromised cellular immune response, infiltration of regulatory T cells, an imbalance between Th1 and Th2 cells, macrophage differentiation, and DC activation.
A tumor's growth micro-environment, created by tumor cells and dominated by tumor-induced interactions, is also composed of various immune cells that fail to appropriately exert anti-tumor effects (Mangino et al., 2016). Recently, some clinical studies have shown that immunotherapy can play a significant role in treating both locally advanced and recurrent/metastatic CESC patients (Rotman et al., 2018).
In recent years, the Cancer Genome Atlas (TCGA) has been leveraged to evaluate immune infiltration and improve poor prognoses for advanced-stage CESC patients (De Felice et al., 2018;Wu et al., 2019). Similarly, performed a comprehensive genome-wide study of gene expression data sets downloaded from TCGA in order to identify promising targets for novel anti-tumor therapies for CESC patients.
In our study, we examined immune cell function in CESC patients using 304 relevant samples from TCGA to investigate tumor micro-environment as it relates to CESC carcinogenesis and invasive progression. This analysis identified 425 differentially expressed genes (DEGs) within the "immune" category, of which 295 were up-regulated and 130 were downregulated. GEPIA-sourced data and TIMER were then employed to analyze hub genes and their relationships to infiltration levels. We established that, in CESC patients, CCL5 expression levels showed significant correlations with levels of several immune cells, which could therefore collectively serve as potential prognostic biomarkers to predict survival in CESC patients. Furthermore, CD8 + T cells, M1 macrophages, M0 macrophages, and M2 macrophages may be candidate targets for CESC therapies. Further investigation should be performed to validate the clinical value of these markers.

Stromal and immune scores in CESC
We examined immune cell function in CESC patients using 304 relevant samples from TCGA (https://portal.gdc.cancer. gov/). An R ESTIMATE algorithm (https://www.r-project. org/) (R version 3.5.1), which infers tumor cellularity by exploiting unique transcriptional profile characteristics (Yoshihara et al., 2013), was applied to these samples and used to calculate their CESC immune scores.

Identification of differently expressed genes (DEGs)
Transcriptional data analysis was performed with the limma package in R, using cut-off criteria of |LogFC| > 2 and pvalue < 0.01 to identify DEGs in high-and low-score groups (Ritchie et al., 2015).

Functional enrichment analysis
Enrichment analysis was performed using the clusterProfiler package in R in order to characterize biological processes (BP), cellular components (CC), and molecular functions (MF) associated with the DEGs (2006; Ashburner et al., 2000;Yu et al., 2012). KEGG pathway enrichment analysis was further used to identify relevant signaling pathways. KEGG is an extensive functional and metabolic pathwayascertaining resource; its exhaustive database compiles abundant information on genomes, biological pathways, diseases, chemical substances, and drugs (Ogata et al., 1999). For these analyses, a p-value < 0.05 was considered statistically significant.

Protein-protein interaction (PPI) network construction and module analysis
The Search Tool for the Retrieval of Interacting Genes database (STRING, https://string-db.org/) resource is designed specifically to review PPI information (Szklarczyk et al., 2015). The DEGs identified previously were submitted to this database to evaluate their potential interactions; those interactions with a combined score > 0.4 were deemed significant and extracted for PPI network construction in Cytoscape (Smoot et al., 2011). Molecular Complex Detection (MCODE) was subsequently used to screen significant PPI network modules, with cutoffs set at degree = 2 and node score = 0.1, k-core set at 2 and maximum depth set at 100 (Bandettini et al., 2012). These modules were then analyzed for functional and pathway enrichment.

Selection and validation of hub genes
Module genes were validated via GEPIA (Gene Expression Profiling Interactive Analysis) and TIMER. GEPIA, made up of 8587 normal and 9736 tumor samples from the GTEx (Genotype-Tissue Expression) and TCGA databases (Tang et al., 2017), is an interactive web application tool that was used to confirm expression level reliability of CESC module genes relative to normal samples. We also explored the prognostic significance of these genes on this platform; those exhibiting significant differences were designated as hub genes, as they were closely correlated with patient prognosis.
Correlation between hub genes and immune cells TIMER, as a comprehensive resource that systematically analyzes immune infiltrates across different cancers (https:// cistrome.shinyapps.io/timer/) , applies a statistical deconvolution methodology (Li et al., 2016) to gene expression profiles. In this way, TIMER can infer the abundance of TIICs (tumor-infiltrating immune cells), which include B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells, in a particular sample. In our study, we assessed their correlations with the expression of CESC hub genes via different gene modules.
Estimation of immune cell type fractions CIBERSORT (http://cibersort.stanford.edu/) (Newman et al., 2015), a data set of 547 genes, allows highly sensitive and specific discrimination of 22 human hematopoietic cell phenotypes, including B cells, T cells, natural killer cells, macrophages, dendritic cells and myeloid subsets. This approach has been previously used on gene expression profiles obtained from microarray and was used here to quantify the proportions of immune cell types that were predicted separately for each gene expression series. For each CESC sample, the sum of all cell-type fraction estimates was set equal to 1.
Gene set enrichment analysis GSEA is a computational method used to determine if a set of genes, defined a priori, exhibit statistically significant differences between two biological states (Subramanian et al., 2005). Our study deployed this approach to first generate a systematic list of all genes relative to CCL5 expression, then elucidated the significant survival difference, which was noticed between the high and low expression groups. The recommended 1000 geneset permutations were performed for each analysis, with CCL5 expression levels serving as phenotype labels.

Results
Stromal and immune scores in CESC Relevant data were downloaded from TCGA on April 21, 2019 and used to assess stromal and immune scores in CESC patients. The 304 CESC cases were divided into top and bottom halves based on their immune scores, presented via the Kaplan-Meier survival curves and statue plot in Fig. 1A. We found that a longer median overall survival was displayed by the high-score cases relative to the low-score group, which was found to be statistically significant.

Identification of DEGs in CESC
Next, correlations between global gene expression profiles and immune scores were identified by executing the limma package in R. The volcano plot and heatmap in Fig. 1B, depict 425 identified DEGs, of which 295 were up-regulated and 130 were down-regulated (Fig. 1B).

Enrichment analysis
ClusterProfiler was deployed for functional and pathway enrichment analysis in order to explore the functions of identified DEGs. With reference to biological processes, GO analysis revealed that DEGs were significantly enriched in FIGURE 1. (A) Validation of ImmuneScore in TCGA. Kaplan-Meier survival curves of OS between high-score and low-score patients; distributions of ImmuneScore, survival status associated with CESC. (B) Volcano plot of DEGs between high-score and low-score patients. Red dots represent up-regulated genes; green dots represent down-regulated genes; black dots represent non-differentially expressed genes. (C) Heatmap of DEGs genes associated with lymphocyte-mediated immunity, as well as genes involved in complement activation, classical pathway and humoral immune response mediated by circulating immunoglobulin. Cell component analysis revealed that these DEGs are generally active on the outer surface of the external side of the plasma membrane, in the immunoglobulin complex (circulating) and in the immunoglobulin complex. Likewise, the molecular function of DEGs was significantly enriched in antigen binding, serine-type endopeptidase activity, and serine-type peptidase activity, as shown in Fig. 2A. Furthermore, KEGG pathway analysis indicated that these DEGs were significantly enriched in associations with graft-versus-host disease, cytokine-cytokine receptor interaction, natural killer cell-mediated cytotoxicity, hematopoietic cell lineage, and antigen processing and presentation, as seen in Fig. 2B.
Hub gene selection GEPIA database mining revealed significant differential expression between CESC and normal tissues, further validating the close relationship between module gene expression levels and disease onset. GEPIA-sourced data from 291 CESC patients were divided into high-and lowexpression groups for overall survival (OS) analysis. We established that up-regulation of CCL5 was correlated with significantly better OS in CESC patients, as seen in Fig. 4A. CCL5 expression levels could, therefore, collectively serve as essential prognostic biomarkers to predict survival in CESC patients. For this reason, we further analyzed CCL5 as the designated hub gene.

Correlation between hub genes and immune infiltration
Tumor-infiltrating lymphocytes serve as independent predictors of sentinel lymph node status and survival in cancers (Azimi et al., 2012;Dunn et al., 2007). We therefore assessed potential CCL5 expression across infiltration levels via TIMER. This analysis revealed significant correlations with several immune cells in CESC, as presented in Figs. 4B, 4C Estimation of immune cell type fractions As seen in Figs. 5A, 5B, CIBERSORT revealed that the fractions of CD8 + T cells and M1 macrophages were consistently higher within the high immune scores group relative to the low scores group, while the fraction of M0 macrophages was significantly lower.
GSEA identifies a CCL5-related signaling pathway Gene Set Enrichment Analysis (GSEA) of low-and high-CCL5 data sets facilitated the identification of signaling pathways differentially activated in CESC. This approach revealed significant enrichment of genes involved in the activation of the immune response, cellular response to interferon-gamma, and the WNT signaling pathway, as seen in Fig. 5C.

Discussion
Cervical cancer is the most common tumor in the female genital tract, as well as the leading cause of morbidity and mortality among women worldwide, with approximately more than 270,000 deaths occurring yearly (Sawaya et al., 2019). Outcomes for CESC patients are suboptimal when the disease is diagnosed at an advanced stage. Survival rates for advanced-inoperable and recurrent CESC patients remain low due to locoregional and distant recurrences (Liontos et al., 2019). HPV infection status, widely known as a cause of CIN and cervical CESC, is considered to play an essential role in patients who are locally advanced and have distant recurrences with dismal prognoses . Despite combined approaches such as surgery and radio-chemotherapy, it remains urgent for researchers to identify effective strategies for these currently incurable CESC patients (Cohen et al., 2019). In particular, immunotherapeutic strategies, which play a key role in the treatment of advanced-inoperable and recurrent CESC, are urgently required. Therefore, our study aims to provide insights into novel prognostic markers, ultimately offering important clues into the relationships that local immune infiltration and tumor microenvironment have with cancer etiology and tumor progression.
In the current study, we assessed immune cell function in CESC using 304 relevant samples from TCGA in order to investigate the relationship of micro-environment to CESC carcinogenesis and invasive progression. The analysis showed that 425 DEGs were detected within the assigned immune scores, including 295 up-regulated and 130 downregulated. Functional and pathway enrichment analysis explored the identified DEGs and found them significantly enriched in the categories of lymphocyte-mediated immunity, complement activation, classical pathway, and humoral immune response mediated by circulating immunoglobulin. PPI network construction and module analysis indicated that up-regulated and down-regulated DEGs were significantly enriched in associations with graftversus-host disease, cytokine-cytokine receptor interaction, allograft rejection, Type-I diabetes mellitus, and natural killer cell-mediated cytotoxicity. CCL5, CCR5, CTLA4, CXCL9, CXCR3, FASLG, GZMA, GZMB, IDO1, IFNG, IL2RA, KLRD1, LAG3, PRF1, PTPRC, TBX21, and TLR8 were the designated leading module genes. Using GEPIAsourced data, we established that the up-regulated expression of CCL5 was correlated with significantly better survival in CESC patients. Therefore, CCL5 expression levels could serve as a novel prognostic biomarker for CESC patients. We further focused on CCL5 as the designated hub gene and went on to establish significant correlations with several types of immune cells in CESC samples. Then, we used CIBERSORT to estimate immune cell type fractions, which showed that CD8 + T cells, M1 macrophages, and M0 macrophages may be candidate targets for CESC therapies. According to GSEA, CCL5-related signaling pathways activated in CESC included activation of the immune response, cellular response to interferon-gamma, and the WNT signaling pathway.
It has been reported that macrophage infiltration has a close relationship with poor prognosis and chemotherapy resistance in numerous types of solid cancers (Tanaka et al., 2019). A number of immunotherapeutic approaches have been employed to limit macrophage infiltration and exert potential antitumor functions (Guilbaud et al., 2019). Macrophages, implicated in suppression of antitumor immunity, can promote the initiation and malignant progression of tumor cells by extravasation, migration, invasion, and metastasis in some solid tumors. However, tumor-associated macrophages (TAMs) serve as natural targets for antitumor therapy in humans via the redifferentiation of macrophages from pro-tumoral phenotypes towards M1-type states (Cassetta and Pollard, 2018;. In our current study, cell-type fractions of M1 macrophages were significantly higher in the high immune scores group, which indicated better overall survival. Recent research has revealed that chemokines CCL2 and CCL5 recruit macrophages which strongly transform the phenotype of macrophages in cervical carcinoma mouse models, dependent on a vaccine-induced influx of CD8 T cells (2015). Similarly, in our study, we established that up-regulated CCL5 was correlated with significantly better survival in CESC patients. Furthermore, CCL5-related signaling pathways activated in CESC included activation of the immune response, cellular response to interferon-gamma and the WNT signaling pathway.
However, our study has several limitations. Scientists explored that some immune checkpoints inhibitors, such as anti-cytotoxic T-lymphocyte antigen-4 (CTLA-4), programmed cell death receptor (PD-1), and programmed cell death ligand 1 (PD-L1) inhibitors, are involved in lifting tumor-induced immune response suppression and revealed highly interrelated rates of PD-L1 + macrophages and regulatory T cells in the lymph nodes of metastatic CESC patients. Additional investigation will be required to further characterize the immune checkpoint function of CCL5, including both in vitro and in vivo studies for mechanistic and prospective insights (Rotman et al., 2018).
In conclusion, we provided an overview of associations between immune infiltration and overall survival in CESC patients. CCL5 was identified as a novel prognostic biomarker for CESC, as well as suggested to serve as a possible immune checkpoint and a potential target for new immunotherapies. Further functional studies could be performed to support our conclusions. In addition, several details about the specific mechanism at play here still need to be confirmed and further clarified.