Investigation of Immune ‐ related Makers Associated with Tumor-inltrating T Cells in Ovarian Cancer

Background: Tumor-inltrating T cells in the tumor microenvironment are the biological basis of immunotherapy and promising predictors of cancer prognosis. Aim: The aim of this study was to investigate immune ‐ related RNAs associated with tumor-inltrating T cells in ovarian cancer (OV). Methods: The gene expression data of patients with OV were downloaded from The Cancer Genome Atlas (TCGA) database. The immune and stromal scores were calculated and the differentially expressed mRNAs (DEGs) were screened. The abundance of six types of inltrating immune cells was investigated, and the immune-related DEGs associated with tumor-inltrating CD4 + and CD8 + T cells were explored by correlation analyses. Subsequently, multiple analyses, i.e., protein-protein interaction (PPI) network analysis, competing endogenous RNA (ceRNA; lncRNA-miRNA-target) network analysis, and small-molecule target network analysis, were performed. Results: In total, 37 and 49 immune-related DEGs of CD4 + and CD8 + T cells were screened, respectively. PPI network results showed that granzyme B (GZMB) was a hub node in the two PPI networks constructed by immune-related DEGs of CD4 + and CD8 + T cells. Moreover, the ceRNA chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E was obtained from the two constructed ceRNA networks related to CD4 + and CD8 + T cells. Survival analysis revealed that key immune-related DEGs of CD4 + and CD8 + T cells, such as GZMB and CD3E, were positively correlated with patient prognosis. Conclusion: GZMB and ceRNAs, such as chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E, may mediate the role of tumor-inltrating T cells in OV. GZMB and CD3E may be used as promising T cell-related biomarkers with prognostic value in OV.

The prognosis of patients with OV is poor due to lack of early symptoms and the high proportion of patients with advanced stage OV (> 60%), which have a ve-year survival rate of just 40% (4,5). Therefore, the exploration of key mechanisms associated with cancer prognosis will facilitate the design of effective therapeutics to improve the clinical outcomes of OV.
Cancer immunotherapy has recently achieved remarkable success in the treatment of late-stage tumors (6,7) and in elucidating the potential role of the tumor immune microenvironment in tumor progression (8). Accumulating evidence has revealed that OV is an immunogenic disease recognized by the host immune system (9). Tumor-in ltrating lymphocytes (TILs), which indicate an ongoing host immune response against cancer cells, exist in the tumor microenvironment and impact the prognostic value of patient-reported outcomes in various cancer entities (10)(11)(12)(13). TILs (e.g., T-cells, macrophages, B-cells, or natural killer cells) are types of white blood cells capable to in ltrate solid tumors. Cytotoxic CD8 + T cells are shown to play a key role in the direct recognition and elimination of tumor cells, and can be recruited and activated by non-cytotoxic CD4 + T cells to promote anti-tumor response. Tumor-in ltrating T cells, especially CD4 + and CD8 + T cells, are considered as powerful prognostic indicators of clinical outcomes in patients with OV (14,15).
In addition, gene expression in tumor tissues is found to be signi cantly affected by the tumor microenvironment, and subsequently in uences the prognosis of patients (16). Immune and stromal cells are important components of the tumor microenvironment that have exhibited diagnostic and prognostic values in tumors (17,18). Increasing studies have revealed tumor microenvironment-associated biomarkers with prognostic value in various cancers on the basis of gene expression data of immune and/or stromal cells (19,20). However, the immune-related markers associated with tumor-in ltrating T cells in OV have not been fully explored.
In the present study, we downloaded the gene expression data of patients with OV from The Cancer Genome Atlas (TCGA) database, and then analyzed the differentially expressed mRNAs (DEGs) on the basis of the immune and stromal scores, which were calculated by the Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) algorithm. The ESTIMATE algorithm could provide scores for estimating the levels of stromal cells and in ltrating immune cells in tumor tissues as well as tumor purity (21). Subsequently, the immune-related DEGs associated with tumor-in ltrating T cells were screened by correlation analysis, followed by protein-protein interaction (PPI) network analysis, competing endogenous RNA (ceRNA; lncRNA-miRNA-target) network analysis, small-molecule target network analysis, and survival analysis. The ow chart of the comprehensive bioinformatic analysis is shown in Supplemental Fig. 1. We aimed to analyze the immune-related RNAs associated with tumor-in ltrating T cells in OV. The ndings of our study will provide theoretical guidance for studying the molecular mechanism of immune escape in OV and exploring the promising T cell immune-related biomarkers for OV therapy and prognosis.

Data source
This study used data available in the public domain. The gene expression data and clinical data of patients with OV (version07-20-2019) were downloaded from the TCGA database (https://xenabrowser.net/). A total of 373 tumor samples were included.

Calculation of immune and stromal scores
In accordance with the annotation information of the platform, the empty probes were removed, and the mean value of multiple probes mapped to the same gene was considered as the mRNA expression level. Immune and stromal scores were then calculated by means of the ESTIMATE algorithm (version 1.0.13) (21) in R package.

DEGs screening
Based on the median value, all samples were assigned as high and low immune score groups as well as high and low stromal score groups. The DEGs in high vs low immune score as well as in high vs low stromal score groups were respectively screened using the limma package (version 3.10.3, http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) (22) in R package. Using the method of Benjamini-Hochberg (23), the adjusted p-value was obtained. The cutoff criteria for DEGs were adjusted to p < 0.05 and |log fold change (FC)| 0.585. The volcano plots for DEGs were visualized using ggplot2 (version: 3.2.1). The intersection DEGs were then obtained by a VENN analysis.

Analysis of immune-related DEGs in T cells
By applying CIBERSORT (26), the abundance of six types of in ltrating immune cells, including CD4 + T cells, CD8 + T cells, B cells, macrophages, neutrophils, and dendritic cells, was calculated. Twenty-two functionally de ned human immune subsets (LM22) were used as a reference expression signature. The histogram for displaying the abundance of immune cell in ltration in the samples was then generated using the ggplot2 package (version: 3.2.1) in R package. The Pearson correlation coe cient between DEGs and the abundance of CD4 + or CD8 + T cells was analyzed, followed by the analysis of their immune-related DEGs with a threshold value of |r| > 0.2.

PPI network construction
The interactions between immune-related DEGs of T cells were retrieved based on the information on the STRING (version: 10.0, http://string-db.org/) (27) database with a combined score of 0.4. The PPI network was visualized using Cytoscape (version: 3.2.0, http://www.cytoscape.org/) (28).
Construction of thechemical small molecule-target network Using Ovarian Neoplasms (Synonym: Ovarian Cancer) as the keyword, disease-related molecules were analyzed using the Comparative Toxicogenomics Database (CTD, http://ctd.mdibl.org/) (31). By comparing the results with those obtained from the ceRNA network, the chemical small molecule-target interactions were obtained, and its related network was constructed using Cytoscape.

Survival analysis of immune-related DEGs in T cells
The prognosis-related clinical data in TCGA were retrieved, including overall survival (OS) and OS status.
Based on the median expression value, the immune-related DEGs in T cells were divided into high and low expression groups. Subsequently, Kaplan-Meier (K-M) survival curves were plotted together with a logrank statistical test. p < 0.05 was set as the signi cant threshold value to select prognosis-related genes.

Results
DEGs were screened based on immune score and stromal score groups In total, 723 DEGs (63 upregulated and 610 downregulated) and 751 DEGs (113 upregulated and 688 downregulated) were screened in high immune score and stromal score groups, in comparison to those in low immune and stromal score groups, respectively.

Function analysis for DEGs
The results of the functional enrichment analysis revealed that the upregulated intersection genes were signi cantly enriched in 46 GO BP terms but were not enriched in any KEGG pathway. Analysis of immune-related DEGs associated with tumorin ltrating T cells The abundance of six types of in ltrating immune cells, including CD4 + and CD8 + T cells, was calculated by CIBERSORT, and the histogram of the abundance of immune cell in ltration in the samples are shown in Fig. 1. By Pearson correlation analyses, 37 and 49 immune-related DEGs of CD4 + and CD8 + T cells, respectively, were obtained. Notably, these immune-related genes were all downregulated.

Interaction analyses between immune-related DEGs in T cells based on PPI network
The PPI network constructed by immune-related DEGs in CD4 + T cells consisted of 36 nodes and 219 interaction relation pairs (Fig. 2A). The PPI network constructed by immune-related DEGs in CD8 + T cells contained 46 nodes and 331 interaction relation pairs (Fig. 2B). Interestingly, granzyme B (GZMB) was identi ed as the hub node in two PPI networks with the highest degree.
Prediction of immune-related miRNAs based on miRNAtarget regulatory network In total, 18 miRNAs that could target six immune-related DEGs of CD4 + T cells were predicted, and 18 miRNA-target pairs were obtained. The miRNA-target regulatory network of CD4 + T cells is shown in Fig. 3A. The apolipoprotein L6 (APOL6) was targeted by the greatest number of miRNAs in this network, such as has-miR-6867-5p. In addition, 12 miRNAs that could target ve immune-related DEGs of CD8 + T cells were predicted, and 12 miRNA-target pairs were obtained. The miRNA-target regulatory network of CD8 + T cells is shown in Fig. 3B. CD3E was targeted by the greatest number of miRNAs, such as has-miR-1249-5p.

Analysis of the immune-related ceRNA network in T cells
Using the DIANA-LncBase v2, 45 lncRNAs that could interact with 16 miRNAs in the miRNA-target regulatory network of CD4 + T cells were obtained. From these, chr22-38_28785274-29006793.1 could interact with the greatest number of miRNAs. By the integration of lncRNA-miRNA and miRNA-target interactions, the ceRNA network in CD4 + T cells (which included 45 lncRNA, 16 miRNAs, 6 immune-related DEGs in CD4 + T cells, and 79 ceRNA interactions) was constructed (Fig. 3C). Moreover, 9 lncRNAs that could interact with 11 miRNAs in the miRNA-target regulatory network of CD8 + T cells were predicted, and chr22-38_28785274-29006793.1 could also interact with the greatest number of miRNAs. By the integration of lncRNA-miRNA and miRNA-target interactions, the ceRNA network of CD8 + T cells (which included 9 lncRNA, 11 miRNAs, 5 immune-related DEGs in CD8 + T cells, and 31 ceRNA interactions) was constructed (Fig. 3D). Intriguingly, chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E was obtained from both ceRNA networks.
Screening chemical small molecules associated with tumor-in ltrating T cells Based on the information of the CTD, 90 chemical small molecule-target pairs in CD4 + T cells were obtained. There were three lncRNAs, ve immune-related DEGs, and 43 chemical small molecules in the small molecule-target network of CD4 + T cells (Fig. 4A). In addition, 86 chemical small molecule-target pairs in CD8 + T cells were obtained. There was one lncRNA, four immune-related DEGs, and 42 chemical small molecules in the small-molecule target network of CD8 + T cells (Fig. 4B). Notably, CD3E was found to be targeted by 14 chemical small molecules, including 1,3-butadiene, arsenic, benzene, benzo(a)pyrene, bisphenol A, cyclophosphamide, dexamethasone, diethylstilbestrol, estradiol, (+)-JQ1 compound, progesterone, resveratrol, tetrachlorodibenzodioxin, and topotecan.

GZMB and CD3E were correlated with the survival of OV patients
Survival analyses were performed to analyze prognosis-related genes. The results showed that 19 and 19 immune-related DEGs of CD4 + and CD8 + T cells, which all included GZMB and CD3E, had a signi cant positive correlation with OS of OV patients (Fig. 5). GZMB and CD3E were common immune-related DEGs associated with both CD4 + and CD8 + T cells.

Discussion
OV is one of the most common oncogenic malignancies among women. Despite multiple progress in OV treatment, 70-80% of patients with OV who initially respond to treatment eventually relapse and die (32). Immunotherapy as a cancer treatment has become a growing eld, but reports of success in treating epithelial OV are limited (33). Tumor-in ltrating T cells in the tumor microenvironment are reported to be the biological basis of immunotherapy and promising predictors of cancer prognosis (14). Therefore, investigation of immune-related RNAs associated with tumor-in ltrating T cells in OV patients could be helpful in nding prognostic markers for OV.
In the present study, we calculated the immune and stromal scores on the basis of gene expression data of patients with OV, and the immune-related DEGs associated with tumor-in ltrating T cells were analyzed. The results showed that 37 and 49 immune-related DEGs of CD4 + and CD8 + T cells, respectively, were screened. Subsequently, the PPI network showed that GZMB was a hub node in both PPI networks constructed by immune-related DEGs of CD4 + and CD8 + T cells. Moreover, the ceRNA chr22-38_28785274-29006793.1/ has-miR-1249-5p/CD3E was obtained from the constructed ceRNA networks associated with CD4 + and CD8 + T cells. Furthermore, survival analysis results revealed that key immunerelated DEGs of CD4 + and CD8 + T cells, such as GZMB and CD3E, were positively correlated with patient prognosis. These results merit further discussion.
GZMB is a serine protease that is primarily expressed by multiple immune cells, including macrophages, lymphocytes, and mast cells (34). It is reported that GZMB is used by activated cytotoxic T lymphocytes to induce apoptosis of target cells (35). In tumor-in ltrating cytotoxic T lymphocytes, type I interferon is found to inhibit tumor growth by inducing STAT3 activation to promote the expression of GZMB (36). Moreover, Prizment et al. has demonstrated that GZMB expression is correlated with the improved survival of patients with colorectal cancer and might be used as a useful prognostic factor (37). In our study, GZMB was a hub node in the two PPI networks constructed by immune-related DEGs of CD4 + and CD8 + T cells. Moreover, GZMB was found to be positively correlated with patient prognosis. Thus, we speculate that GZMB could be used as a potential prognostic marker of tumor-in ltrating T cells for OV.
CD3E is one of the four invariant chains physically associated with the T-cell receptor (TCR), and the TCR-CD3 complex has been reported to exert a signi cant function in the immune system (38). The combined immunoscore of CD103 and CD3 may be used to identify patients with early relapse or longterm survival in serous OV (39). We found that CD3E was targeted by the greatest number of miRNAs in the miRNA-target regulatory network, and CD3E showed signi cant positive correlation with patient OS. Therefore, CD3E may also be a promising prognostic marker of tumor-in ltrating T cells for OV. In addition, CD3E was found to be targeted by several chemical small molecules, such as 1,3-butadiene and benzene. It is reported that 1,3-butadiene and benzene exposure can increase lifetime risk of cancer (40).
Thus, we speculate that monitoring the expression of CD3E could be used to predict lifetime risk of cancer caused by exposure to chemical small molecules.
Furthermore, miRNAs able to target CD3 subunit molecules and affect the immune system are considered as promising biomarkers for early cancer detection (41). In our study, we found that CD3E was targeted by multiple miRNAs, such as has-miR-1249-5p. miR-1249 is found to act as either an oncogene or tumor suppressor in different tumor types. For instance, miR-1249 is upregulated in patients with small cell carcinoma of the esophagus, and its enhanced expression is associated with tumor relapse (42). Nevertheless, miR-1249 is shown to be markedly downregulated in colorectal cancer tissues, and it may serve as a novel therapeutic target (43). Our data predicted the target relationship between miR-1249-5p and CD3E. Although little is known about miR-1249 in OV, our result prompts us to speculate that miR-1249-5p may be a potential regulator mediating the key role of tumor-in ltrating T lymphocytes in OV via the targeting of CD3E. In addition, existing evidence has con rmed that lncRNAs function as ceRNAs to competitively bind to the miRNA binding sites, consequently reducing their regulatory effect on target mRNAs (44,45). Our results showed that the ceRNA chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E was obtained from the two constructed ceRNA networks related to CD4 + and CD8 + T cells. Therefore, we speculate that ceRNAs, such as chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E, could be a target of tumor-in ltrating T cells and are involved in the immune escape of tumors in OV. It would be interesting to further investigate the clinical signi cance of ceRNA chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E in OV.
In conclusion, our data reveal that tumor-in ltrating T cells may contribute to the immune escape of tumor cells in OV by regulating GZMB and ceRNAs, such as chr22-38_28785274-29006793.1/has-miR-1249-5p/CD3E. GZMB and CD3E may serve as promising T cell immune-related biomarkers with prognostic value in OV.