System analysis based on the T cell exhaustion‑related genes identifies CD38 as a novel therapy target for ovarian cancer

Introduction

Ovarian cancer (OV) is one of the deadliest and most aggressive cancers in women, and the incidence has increased in recent years [1]. Since the early symptoms of OV are insidious and the cancer develops rapidly, most patients are not diagnosed until the late stage of the disease [2]. Furthermore, the heterogenous nature of ovarian tumors restricts the predictive accuracy of existing prognostic markers, which also worsens patient prognosis [3]. Currently, surgery and chemoradiotherapy are the most common treatment strategies for OV. However, the considerable adverse effects of these approaches severely reduce the patients’ quality of life [4]. Given that the heterogeneity of ovarian tumors is a key factor affecting cancer progression and overall survival [5], and results in different responses of individual patients to the same treatment [6], it is crucial to use single-cell technologies to identify effective prognostic biomarkers [7,8]. In fact, analyses of tumors at the single-cell level have provided new insights into the molecular mechanisms underlying carcinogenesis, and revealed novel therapeutic possibilities [9–13]. A recent study examined the heterogeneity of OV by single-cell RNA-sequencing (scRNA-seq) [14]. In addition, Hu et al. were able to predict cancer behavior based on the analysis of individual fallopian tube epithelial cells [15]. Therefore, scRNA-seq studies can also improve our understanding of OV.

In recent years, immunotherapy has emerged as one of the most promising modalities of cancer treatment [16–18]. The tumor-infiltrating CD8+ T cells, the main effector cells of anti-cancer immunity [19,20], show impaired viability and proliferation, increased apoptosis, decreased production of effector cytokines, and lower cell killing ability [21]. Studies show that the majority of the anti-tumor CD8+ T cells are in the exhausted state [22], and do not express interleukin-2 (IL-2) or tumor necrosis factor α (TNF-α) [22,23]. CD8+ T cell exhaustion is a significant prognostic factor in liver cancer and colorectal cancer [24,25]. However, the prognostic relevance of CD8+ T cell exhaustion in OV has not been established so far.

In this study, we analyzed the single-cell transcriptional profiles obtained from OV patients, and identified four T cell-associated subtypes. The signature genes of the exhausted CD8+ T cell subset were screened to establish the T-cell related gene score (TRS) by random forest plot method. Based in the TRS, the patients were divided into the low- and high-risk groups, which differed in terms of prognosis, immune cell infiltration and drug sensitivity. Furthermore, the role of CD38, one of the signature genes in the TRS, in the malignant potential of OV cells was verified through in vitro experiments.

This study provides a more comprehensive understanding of the heterogeneity of T cells in ovarian tumors and the clinical significance of T cell exhaustion in OV. The prognostic model based on the TRS can help develop more precise and effective therapies.

Materials and Methods

Data collection

The single-cell sequencing data of 5 OV patients (GSE154600) was retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. The RNA-seq data and accompanying clinical information of 378 OV patients were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). In addition, the GSE140082 dataset including the RNA-seq and clinical data of 379 OV patients was also downloaded from the GEO database as the validation cohort. All available datasets used in this study were based on published reports that had received ethical approval. Dataset inclusion criteria: (1) Pathological diagnosis of ovarian cancer; (2) Complete transcriptome and clinical prognostic data; Exclusion criteria: (1) The pathological diagnosis was unclear; (2) Survival time less than 1 day; (3) The transcriptome information did not contain T cell exhaustion related genes.

Data filtering and correction

The scRNA-seq data was analyzed using the “Seurat”, “SingleR” software packages. The cells with unique feature counts >5000 or <200, or with mitochondrial counts >5% were filtered out. The data of each retained cell was normalized using the default parameters of Seurat’s ‘NormalizeData’ function on feature-expression measurements, and then transferred to a combined Seurat object via the Harmony package. The variable genes were scaled and subjected to principal component analysis (PCA). The significant principal components (PCs) were selected by the “RunUMAP” function (min. dist = 0.2, n. neighbors = 20) and the “FindClusters” function (resolution = 0.5) for umap analysis and cluster analysis.

Cell annotation

The individual cell types were annotated by two modalities. The T-cell subsets were initially clustered using SingleR, an automated annotation method for scRNA-seq data [26], and the identity of each cluster was further determined by manually searching the cell markers (http://biocc.hrbmu.edu.cn/CellMarker/). The differentially expressed genes (DEGs) in the CD8+ exhausted T cell subset were screened using the FindAllMarker function in the R “Seurat” package, with logFC > 0.3 and p < 0.05 as the thresholds.

Enriched pathway analysis

The differences in the biological functions of the T cell-related subsets were determined by gene set enrichment analysis (GSEA), gene ontology (GO) and KEGG pathway analyses using R ‘ABGSEase’ and ‘clusterprofiler’ packages with Hallmark, GO and KEGG reference gene sets. Maximum gene set size (maxGSSize) was set to 100, minimum gene set size (minGSSize) to 50, and p-value < 0.05 was considered statistically significant.

Construction of prognostic model

The DEGs in the CD8+ exhausted T cell subsets in the TCGA and GEO cohorts were corrected to the same level, and the prognostic genes were identified through univariate Cox regression analysis (p < 0.01). A random forest model was then built using the R package “randomForest”. The mean model inaccuracies for all genes were first estimated. The optimal number of trees in the random forest model was set at 500, and the dimensional effect sizes of the model were determined using the Gini coefficient method. Genes with importance indices higher than 0.6 were selected as disease hub genes. TRS was calculated using the following formula:

TRS= Coef1×Geneexpression1+Coef2×Geneexpression2+⋯Coefn×Geneexpressionn

The Cox regression coefficient represents the prognostic value of each gene. Expression values for genes correspond to those of their modeled counterparts. Based on the median TRS, OV patients were classified as high-risk or low-risk. The R “survival” package [27,28] was used to plot Kaplan-Meier survival curves. Based on clinical features, the “rms” package of R software was used to develop a nomogram to predict survival of individual patients [29].

Immune cell infiltration analysis

The infiltration of different immune cell populations was quantified based on the transcriptional profiles of patients using MCPcounter [30] and xCell [31] algorithms. In addition, the absolute content of the four T cell subsets were derived from the reference single-cell dataset. The data was visualized using boxplots, heat maps, and scatter plots.

Drug sensitivity analysis

The sensitivity score of the drugs in the Genomics of Drug Sensitivity in Cancer (GDSC) database was calculated with the R package ‘oncoPredict’ [32]. Wilcox’s test was used to identify the candidate drugs with p < 0.05 as the threshold for statistical significance.

Cell culture

The A2780 OV cell line and IOSE80 normal ovarian cell line were obtained from the American Type Culture Collection (ATCC, Rockville, MD, USA). The cells were cultured in RPMI-1640 medium supplemented with 10% FBS at 37°C and 5% CO2.

Quantitative real-time PCR (qRT-PCR)

RNA was extracted from the cells using TRIzol reagent (Takara, Japan) and reverse transcribed into cDNA using SuperScript II Reverse Transcriptase (Invitrogen, USA), The relative expression level of CD38 was analyzed by qRT-PCR and normalized to GAPDH. The primers are listed in Suppl. Table S1.

Cell transfection

The CD38 siRNA was mixed with the RNAiMAX transfection reagent in Opti-MEM for 20 min at room temperature, and then transfected into the OV cells. After 12 h, the medium was discarded and RPMI-1640 medium was added. The procedure was repeated after 12 h. The CD38-siRNA target sequences were as follows: sense 5′-GCCUGAUCUAUACUCAAAUTT-3′ and antisense 5′-AUUUGAGUAUAGAUCAGGCTT-3′. CD38 knockdown was verified by qRT-PCR.

Apoptosis analysis

The transfected cells were washed with PBS, and harvested with trypsin digestion solution containing no EDTA (Solarbio, China). After centrifuging at 1000 rpm for 5 min, the cells were stained with 7-AAD (BD Biosciences, USA) and annexin-APC (BD Biosciences, USA) for 15 min.

Cell invasion assay

The transfected cells were seeded on matrix gel-coated upper chamber of Transwell inserts in serum-free RPMI-1640 medium, and the lower chambers were filled with complete medium. After culturing for 24 h, the cells remaining on the upper membrane surface were removed with cotton swabs, and those that had adhered to the bottom surface were fixed with 4% paraformaldehyde and then stained with crystal violet. The cells were observed under a light microscope, and counted using ImageJ.

Statistical analysis

R (version 4.1.0) program was used for statistical analyses. Overall survival was evaluated using Kaplan-Meier (KM) and log-rank tests, and differences among subgroups were determined with Wilcox test and Kruskal test.

Results

Cell clustering analysis of OV

The study design is outlined in the flow chart in Fig. 1. The scRNA-seq data of obtained from 5 OV patients was screened, and the cells were annotated in terms of spatial distribution (Fig. 2A, left), as well as characteristic markers. Six cell types were annotated, including fibroblasts, T cells, myeloid cells, epithelial cells, endothelial cells, and B cells (Fig. 2A, right). The expression profiles of the eight classical marker genes are shown in Fig. 2B, and their expression in each cell type is shown in the bubble plots in Fig. 2C. Based on the absolute counts, the T cell subsets were predominant in almost all samples, suggesting a crucial role of these cells in OV progression (Fig. 2D). In addition, the content of the other cell types in each cell subset was significantly different, which confirms the cellular heterogeneity of OV.

Figure 1: Flow chart shows the study design.

Figure 2: Single cell cluster analysis of OV. (A) The tSNE plots show a dimensionality reduction cluster analysis of single cells from five OV samples annotated as a total of six cells. (B) The tSNE plot shows the expression of signature genes in single-cell clusters. (C) Bubble plots showing expression of signature genes. (D) Histogram showing the amount of each cell type in the sample.

Dimensional cluster analysis of T cell subsets

Following dimensionality reduction cluster analysis of T cell subsets, we annotated four subsets, including CD8+ T cells, CD4+ T conv cells, regulatory T cells (Tregs) and CD8+ exhausted T cells (Fig. 3A). The expression of signature genes in each cell subset is shown in the bubble plots in Fig. 3B. The CD8+ T cells were the most abundant subset among the samples (Fig. 3C), whereas the CD8+ exhausted T cells were the least abundant in each sample. Analysis of the differentially expressed genes in the T cell subsets indicated that CXCL13 and HAVCR2 were significantly upregulated in the CD8+ exhausted T cells (Fig. 3D).

Figure 3: Dimensional cluster analysis of T cell subsets. (A) The tSNE plot shows the dimensionality reduction clustering of T cell subsets, annotated as a total of four T cell-associated subsets. (B) Bubble plots showing expression of signature genes. (C) Histogram showing the content of each T-cell related subset in the sample. (D) Violin plot showing expression of two signature genes in CD8+ exhausted T cells.

Enriched pathway analysis

GSEA of tumor-associated pathways was performed to further explore the differences in the biological functions of the T cell subsets. The myogenesis and epithelial mesenchymal transition pathways were significantly activated in CD8+ T cells, and TNFA signaling via NFKB pathway in the CD4+ T conv cells. The Tregs were significantly associated with the oxidative phosphorylation, pancreatic beta cells and peroxidase pathways, whereas the CD8+ exhausted T cells showed significant enrichment of the interferon alpha response, interferon gamma response, MYC target V1, oxidative phosphorylation, E2F target and G2M checkpoint pathways (Figs. 4A and 4B). Subsequently, we analyzed tumor-associated pathways in different cell populations based on PROGENy score. As shown in Fig. 4C, The JAK-STAT and MAPK pathways were activated, and the p53 and PI3K pathways were inhibited in the CD8+ exhausted T cells. Receptor ligand pair analysis between CD8+ exhausted T cells and other T cell subsets revealed high enrichment of the HLA−A−CD8A, HLA−B−CD8A, and HLA−C−CD8A pairs (Fig. 4D). Given the central role of HLA molecules in immunity, our findings indicate that this subset may influence the response to immunotherapy. Furthermore, GO analysis of the DEGs between CD8+ exhausted T cells and controls showed significant enrichment of ATP metabolic process, aerobic respiration, cellular respiration and oxidative phosphorylation in biological processes (BP), redox-driven activity transmembrane transporter activity and NADH dehydrogenase activity in molecular function (MF), and mitochondrial inner membrane, mitochondrial protein complex and mitochondrial inner membrane protein complex in cellular compartment (CC) terms (Fig. 4E). Finally, the KEGG pathways significantly associated with the DEGs included those related to Huntington’s disease, prion disease, thermogenesis, amyotrophic lateral sclerosis, and other neurodegenerative diseases (Fig. 4F).

Figure 4: Enriched pathway analysis. (A, B) Heat map showing the expression of 50 hallmark tumor-associated pathways in each T cell subset. (C) Heat map showing the expression of 14 key pathways in each T cell subset. (D) Bubble plots showing CD8+ T-cell exhaustion vs. expression of ligand pairs for other cell subsets receptors. (E) Histogram of GO enriched pathways for DEGs in the CD8+ exhausted T cells. (F) Histogram of KEGG enriched pathways for the DEGs.

Prognostic model for OV

To identify the prognostic genes associated with CD8+ T cell exhaustion, we performed univariate Cox regression analysis of the DEGs in this subset. We screened 47 prognostic genes that were incorporated into a random forest classifier. Seven genes, including MOB1A, TIMM8B, CD38, HINT2, NAA38, DLEU2 and TXNDC17, were finally identified and prognostic models were constructed using a random forest classifier (Fig. 5A). The TRS was calculated as follows:

Figure 5: Construction of prognostic model of OV. (A) The number of selected trees impacts the error rate. The x-axis represents the number of decision trees, and the y-axis represents the error rate. Random forest classifier results based on Gini coefficient. The x-axis represents genetic variables and the y-axis represents importance indices. Seven genes with prognostic significance were selected from 47 genes based on random forest and prognostic models were constructed. (B) KM survival curves and scatter plot of survival distributions for the TCGA cohort. (C) KM survival curves and scatter plot for the GEO validation cohort. (D) Nomogram of OV prognostic model based on stage, age and TRS. (E) Pie charts show the correlation between each modeled gene. (F) Box plots showing the expression of each modeled gene between high- and low-risk groups.

TRS= 0.25516×MOB1A+(−0.07000)×TIMM8B+(−0.12049)×HINT2+(−0.06368)×NAA38+(−0.15835)×CD38+(−0.04413)×TXNDC17+0.04650×DLEU2

According to the median risk score, the patients in TCGA cohort were divided into the low-risk and high-risk groups. As shown in Fig. 5B, the patients with low-risk score had significantly better OS. Using the median TRS of the training set, we classified patients in the external GSE140082 set as high- or low-risk. As expected, the low-risk group had significantly better OS compared to the high-risk group, thus validating our prognostic model (Fig. 5C). To accurately estimate the survival status of individual OV patients, we developed a predictive nomogram based on TRS, age and tumor stage. As shown in Fig. 5D, the TRS-based nomogram can accurately predict the survival of OV patients.

Differentiation of high-risk and low-risk patients based on immune landscape

Except for CD38, the other genes in the prognostic model showed significant correlation with each other (Fig. 5E). In addition, TIMM8B, CD38, HINT2, NAA38 and TXNDC17 were significantly upregulated in the low-risk patients, whereas higher levels of MOB1A were detected in the high-risk group (Fig. 5F).

Since the immune landscape plays a key role in tumor development and progression, we compared the immune infiltration in the high-risk and low-risk groups using MCPcounter and xCell. According to the MCPcounter algorithm, patients in the low-risk group had significantly higher number of infiltrating CD8+ T cells, cytotoxic lymphocytes, myeloid dendritic cells, and NK cells compared to patients in the high-risk group, whereas the latter had significantly higher abundance of neutrophils (Fig. 6A). We also detected a significant association between CD38 and most immune-related genes (Fig. 6B). In addition, the B lineage cells and cytotoxic lymphocytes showed a significant negative correlation with the TRS (Fig. 6C). The xCell algorithm also indicated significant differences between the immune cell composition of the high- and low-risk groups. As shown in Fig. 6D, the aDC cells, M1 macrophages, M2 macrophages, Th1 cells and Th2 cells were significantly more abundant in the high-risk vs. the low-risk patients. Therefore, the poor prognosis in the high-risk group may be due to the distinct immune landscape. In addition, we found a significant negative correlation between CD4+ memory T cells and TRS, and a significant positive correlation between eosinophils and mast cells and TRS (Fig. 6E).

Figure 6: Immune cell infiltration analysis. (A) Boxplots showing differences in the infiltration of each cell subset in the high- and low-risk groups based on based on MCPcounter algorithm. (B) Heat map showing the correlation between each modeled gene and cell subset. (C) Correlation between the two cell types based on MCPcounter algorithm. (D) Box plots showing differences in the infiltration of each cell subset in the high- and low-risk groups based on the xCell algorithm. (E) Correlation between the three cell types based on the results of immune cell infiltration analysis by xCell algorithm.

CD38 knockdown in OV cells inhibited their malignant potential

CD38 was the only model gene that was significantly associated with prognostic indicators in OV patients. In addition, CD38 expression levels were significantly higher in the OV tissues compared to that in the normal ovarian tissues (Fig. 7A). Consistent with this, the OV cell line A2780 expressed higher levels of CD38 compared to the normal ovarian cells IOSE80 (Fig. 7B). To further evaluate the biological function of CD38 in OV, we silenced the gene in A2780 cells. As shown in Fig. 7C, CD38 was significantly downregulated in the OV cells following gene knockdown. In addition, knocking down CD38 significantly increased the apoptosis rates in the OV cells (Fig. 7D), and inhibited their invasion ability in the Transwell assay (Fig. 7E). Thus, CD38 is a potential therapeutic target for OV.

Figure 7: CD38 inhibited the malignant potential of OV cells. (A) Box-line plot showing that CD38 expression is significantly higher in the OV group. (B) Histogram showing the expression levels of CD38 gene in normal ovarian cells and OV cells. (C) Histogram showing knockdown of CD38 gene expression levels in OV cells. (D) Flow cytometry scatter plot showing apoptosis rates of OV cells after CD38 knockdown. (E) Representative images and histogram of Transwell assay showing invasion ability of OV cells after CD38 knockdown.

Identification of candidate drugs for OV

In addition, we examined the relationship between TRS and the expression of immunomodulators. the TRS was negatively correlated with most immune checkpoints which suggests that TRS may be used to predict the efficacy of immunotherapy (Fig. 8A). To further evaluate the chemosensitivity of OV patients, we assessed the sensitivity scores for each compound. The sensitivity of six drugs was significantly different between the high and low risk groups. The patients in the high-risk group were more sensitive to AZD6738, Ribociclib, AZD3759, LCL161 and Navitoclax, whereas the low-risk group showed greater sensitivity to BI−2536 (Fig. 8B). Furthermore, there was a significant positive correlation between MOB1A and most of these drugs, and CD38 showed correlation with BI−2536, RO−3306 and Tozasertib (Fig. 8C).

Figure 8: Drug sensitivity analysis. (A) Histogram showing the correlation between the TRS and immune checkpoints. (B) Box plots showing drug sensitivity in the high-risk and low-risk groups. (C) Correlation between each modeled gene and potential drugs.

Discussion

OV is one of the deadliest and most aggressive cancers in women, and the incidence has increased in recent years [1]. Given the insidious nature and rapid progression of OV, most patients are not diagnosed until the disease is in the advanced stage [2]. The poor prognosis of OV patients may result from tumor heterogeneity, which restricts prognostic prediction [3]. Currently, surgery and chemoradiotherapy are the most common approaches for treating OV, although the side effects severely reduce the quality of life of the patients [4]. In order to improve treatment outcomes, it is necessary to identify biomarkers for the early diagnosis of OV.

Immunotherapy has emerged as a promising treatment strategy for cancer in recent years, which in turn has spurred the development of immune-related prognostic models. Several studies have reported construction of cancer prognosis models based on genes involved in immune responses, cuproptosis, endoplasmic reticulum stress, angiogenesis, and other processes involved in tumorigenesis, which have shown good prognostic performance, as well as predictive accuracy for immune function [33–37]. Due to considerable intra-tumoral heterogeneity, it is challenging to identify genetic variants through bulk mRNA sequencing. Although the identification of oncogenic drivers by this approach has facilitated the classification of some tumor types, and helped develop new prognostic models, singe cell analysis is crucial to overcome tumor heterogeneity [36,38,39]. ScRNA-seq has been instrumental in identifying the subpopulations within heterogenous tumors, and dissect the molecular mechanisms underlying tumor biology [40,41]. Signature genes screened using single-cell transcriptomic data have been used to establish prognostic models that can not only obviate the limitations of tumor heterogeneity, but also predict prognosis with greater accuracy based on target cell subsets [42–44].

In this study, we identified four distinct T cell sub-populations in ovarian tumors by analyzing the single-cell transcriptional profiles of five OV patients. Pathways related to oxidative phosphorylation and G2M checkpoint were significantly activated in the CD8+ exhausted T cell subset. Oxidative phosphorylation is an integral process in cancer cell growth and development, and a potential therapeutic target [45–47]. In addition, centrosome-associated G2/M checkpoint regulators have also been identified as potential targets for cancer treatment [48,49]. We developed a prognostic model for OV based on seven marker genes for CD8+ T cell exhaustion, including MOB1A, TIMM8B, CD38, HINT2, NAA38, DLEU2 and TXNDC17. The patients were classified as high-risk or low-risk based on the median value of the TRS. Patients in the low-risk group had significantly better prognosis in both the GEO and TCGA cohorts compared to those in the high-risk group. A protein encoded by the MOB1A gene is a component of the Hippo signaling pathway, which promotes apoptosis and controls organ size and tumor growth [50,51]. CD38 encodes a transmembrane glycoprotein that is expressed in several tissues and cells, including those of the immune system. The CD38 protein is a prognostic marker of chronic lymphocytic leukemia and a therapeutic target in multiple myeloma, and has an established role as an immunomodulator in cancer [52]. In addition, CD38 can also predict the prognosis of epithelial OV by enhancing immune infiltration and anti-tumor immunity in the microenvironment [53]. HINT2 encodes a histidine triplet protein that plays a role in various cancers, and HINT2 down-regulation can promote colorectal cancer migration and metastasis [54]. NAA38 has not yet been investigated for its specific role in cancer development, but as part of the NatC complex, it inhibits the apoptotic process and is also a potential target [55]. Apart from modulating cell proliferation and invasion, DLEU2 regulates miRNA levels in pancreatic and non-small cell lung cancers [56,57]. TXNDC17 can induce autophagy and promote paclitaxel resistance in OV cells via the tumor necrosis factor-mediated signaling pathways [58].

Immune cell infiltration analysis by MCPcounter and xCell algorithms indicated significant differences between the immune landscapes of the high-risk and low-risk groups. Interestingly, the aDCs, M1 macrophages, M2 macrophages, Th1 cells and Th2 cells were significantly more abundant in high-risk patients than in low-risk patients. Macrophages in the tumor microenvironment can polarize to the pro-tumorigenic M1 phenotype, or the anti-angiogenic and pro-apoptotic M2 phenotype depending on the stimuli [59,60]. OV cells promote polarization of co-cultured macrophages in vitro, which may regulate the cytokine profile of the tumor microenvironment [61]. In addition, OV patients with high Th2 mRNA levels have a better prognosis than those with high Th1 mRNA levels [62]. Therefore, we speculate that the difference in prognosis between high-risk and low-risk groups may originate from differences in the immune microenvironment.

Finally, we identified six potential agents that may be effective for OV treatment. Patients in the high-risk group showed greater sensitivity to AZD6738, Ribociclib, AZD3759, LCL161 and Navitoclax, whereas BI-2536 was more suitable for patients in the low-risk group. AZD6738 is an ATR inhibitor, and can overcome chemoresistance in OV cells when used in combination with belitecan [63]. Ribociclib has been effective in recurrent ER-positive OV, and is currently in phase II trials for OV treatment [64]. Therefore, the TRS signature may also help in the selection of optimal drugs for individual OV patients.

The role of T cell exhaustion in ovarian cancer was rarely described in the literature before this paper. Zhou’s team [65] constructed a pan-cancer blueprint for the heterogeneity of T cell exhaustion based on the transcriptome information of 32 tumor types to predict tumor prognosis and therapeutic effect. However, the transcriptome information does not contain the specific information of a single cell line. Therefore, the single-cell transcriptome information was added in this paper to further reveal the mechanism of T cell exhaustion. Based on the work of Professor Zhou’s team, T cell depletion genes associated with ovarian cancer were further explored to construct a model for predicting prognosis and therapeutic effect. Meanwhile, there are some limitations in our study. The TRS was developed and validated using data from public databases, and was not further validated in real-world cohorts. Future research will focus on more robust validation of TRS in multicenter, large cohort data.

Conclusion

We developed a T cell exhaustion-based gene signature to predict the prognosis of OV patients, and identified CD38 as a potential therapeutic target.

Acknowledgement: The authors thank all the staff for their contributions.

Funding Statement: This experiment was supported by the following funds: The Shanghai Municipal Key Clinical Specialty (No. shslczdzk06302). National Natural Science Foundation of China (No. 82103029). The Project of The Science and Technology Commission of Shanghai Municipality (No. 21ZR1469500). The Shanghai Jiao Tong University Medicine-Engineering Fund (No. YG2021QN137).

Author Contributions: Tianming Shi was responsible for project design and implementation. Mi Han contributed to the conception of this investigation. Rongrong Yan contributed to data analysis.

Availability of Data and Materials: The data in this study are from public data set TCGA and GEO. In addition, it also contains the experimental data of our research group.

Ethics Approval: No ethical approval was required for this study.

Conflicts of Interest: The authors declare that the study was conducted without any commercial or financial ties.

References

1. Siegel, R. L., Miller, K. D., Jemal, A. (2019). Cancer statistics, 2019. CA: A Cancer Journal for Clinicians, 69(1), 7–34. [Google Scholar] [PubMed]

2. He, C., Huang, F., Zhang, K., Wei, J., Hu, K. et al. (2021). Establishment and validation of an RNA binding protein-associated prognostic model for ovarian cancer. Journal of Ovarian Research, 14(1), 1–11. https://doi.org/10.1186/s13048-021-00777-1 [Google Scholar] [PubMed] [CrossRef]

3. Izar, B., Tirosh, I., Stover, E. H., Wakiro, I., Cuoco, M. S. et al. (2020). A single-cell landscape of high-grade serous ovarian cancer. Nature Medicine, 26(8), 1271–1279. https://doi.org/10.1038/s41591-020-0926-0 [Google Scholar] [PubMed] [CrossRef]

4. Elies, A., Rivière, S., Pouget, N., Becette, V., Dubot, C. et al. (2018). The role of neoadjuvant chemotherapy in ovarian cancer. Expert Review of Anticancer Therapy, 18(6), 555–566. https://doi.org/10.1080/14737140.2018.1458614 [Google Scholar] [PubMed] [CrossRef]

5. Zhao, H., Gao, Y., Miao, J., Chen, S., Li, J. et al. (2021). Single-cell RNA-seq highlights a specific carcinoembryonic cluster in ovarian cancer. Cell Death & Disease, 12(11), 1–11. https://doi.org/10.1038/s41419-021-04358-4 [Google Scholar] [PubMed] [CrossRef]

6. Nguyen, H. M., Vessella, R. L., Morrissey, C., Brown, L. G., Coleman, I. M. et al. (2017). LuCaP prostate cancer patient-derived xenografts reflect the molecular heterogeneity of advanced disease and serve as models for evaluating cancer therapeutics. The Prostate, 77(6), 654–671. https://doi.org/10.1002/pros.23313 [Google Scholar] [PubMed] [CrossRef]

7. Lawson, D. A., Kessenbrock, K., Davis, R. T., Pervolarakis, N., Werb, Z. (2018). Tumour heterogeneity and metastasis at single-cell resolution. Nature Cell Biology, 20(12), 1349–1360. https://doi.org/10.1038/s41556-018-0236-7 [Google Scholar] [PubMed] [CrossRef]

8. Yuan, X., Wang, J., Huang, Y., Shangguan, D., Zhang, P. (2021). Single-cell profiling to explore immunological heterogeneity of tumor microenvironment in breast cancer. Frontiers in Immunology, 12, 643692. https://doi.org/10.3389/fimmu.2021.643692 [Google Scholar] [PubMed] [CrossRef]

9. Navin, N., Hicks, J. (2011). Future medical applications of single-cell sequencing in cancer. Genome Medicine, 3(5), 1–12. https://doi.org/10.1186/gm247 [Google Scholar] [PubMed] [CrossRef]

10. Guo, X., Zhang, Y., Zheng, L., Zheng, C., Song, J. et al. (2018). Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nature Medicine, 24(7), 978–985. https://doi.org/10.1038/s41591-018-0045-3 [Google Scholar] [PubMed] [CrossRef]

11. Savas, P., Virassamy, B., Ye, C., Salim, A., Mintoff, C. P. et al. (2018). Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis. Nature Medicine, 24(7), 986–993. https://doi.org/10.1038/s41591-018-0078-7 [Google Scholar] [PubMed] [CrossRef]

12. Ren, X., Kang, B., Zhang, Z. (2018). Understanding tumor ecosystems by single-cell sequencing: Promises and limitations. Genome Biology, 19(1), 1–14. https://doi.org/10.1186/s13059-018-1593-z [Google Scholar] [PubMed] [CrossRef]

13. Bartoschek, M., Oskolkov, N., Bocci, M., Lövrot, J., Larsson, C. et al. (2018). Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing. Nature Communications, 9(1), 1–13. https://doi.org/10.1038/s41467-018-07582-3 [Google Scholar] [PubMed] [CrossRef]

14. Shih, A. J., Menzin, A., Whyte, J., Lovecchio, J., Liew, A. et al. (2018). Identification of grade and origin specific cell populations in serous epithelial ovarian cancer by single cell RNA-seq. PLoS One, 13(11), e0206785. https://doi.org/10.1371/journal.pone.0206785 [Google Scholar] [PubMed] [CrossRef]

15. Liu, Y., Wu, J., Huang, W., Weng, S., Wang, B. et al. (2021). Development and validation of a hypoxia-immune-based microenvironment gene signature for risk stratification in gastric cancer. Journal of Translational Medicine, 18(1), 201. https://doi.org/10.1186/s12967-020-02366-0 [Google Scholar] [PubMed] [CrossRef]

16. Schuster, M., Nechansky, A., Kircheis, R. (2006). Cancer immunotherapy. Biotechnology Journal, 1(2), 138–147. https://doi.org/10.1002/(ISSN)1860-7314 [Google Scholar] [CrossRef]

17. Couzin-Frankel, J. (2013). Breakthrough of the year 2013. Cancer immunotherapy. Science, 342(6165), 1432–1433. https://doi.org/10.1126/science.342.6165.1432 [Google Scholar] [PubMed] [CrossRef]

18. Farkona, S., Diamandis, E. P., Blasutig, I. M. (2016). Cancer immunotherapy: The beginning of the end of cancer? BMC Medicine, 14(1), 1–18. https://doi.org/10.1186/s12916-016-0623-5 [Google Scholar] [PubMed] [CrossRef]

19. Dolina, J. S., Van Braeckel-Budimir, N., Thomas, G. D., Salek-Ardakani, S. (2021). CD8+ T cell exhaustion in cancer. Frontiers in Immunology, 12, 715234. https://doi.org/10.3389/fimmu.2021.715234 [Google Scholar] [PubMed] [CrossRef]

20. Ma, X., Bi, E., Lu, Y., Su, P., Huang, C. et al. (2019). Cholesterol induces CD8+ T cell exhaustion in the tumor microenvironment. Cell Metabolism, 30(1), 143–156.e5. https://doi.org/10.1016/j.cmet.2019.04.002 [Google Scholar] [PubMed] [CrossRef]

21. Hashimoto, M., Kamphorst, A. O., Im, S. J., Kissick, H. T., Pillai, R. N. et al. (2018). CD8 T cell exhaustion in chronic infection and cancer: Opportunities for interventions. Annual Review of Medicine, 69(1), 301–318. https://doi.org/10.1146/annurev-med-012017-043208 [Google Scholar] [PubMed] [CrossRef]

22. He, Q. F., Xu, Y., Li, J., Huang, Z. M., Li, X. H. et al. (2019). CD8+ T-cell exhaustion in cancer: Mechanisms and new area for cancer immunotherapy. Briefings in Functional Genomics, 18(2), 99–106. https://doi.org/10.1093/bfgp/ely006 [Google Scholar] [PubMed] [CrossRef]

23. McLane, L. M., Abdel-Hakeem, M. S., Wherry, E. J. (2019). CD8 T cell exhaustion during chronic viral infection and cancer. Annual Review of Immunology, 37(1), 457–495. https://doi.org/10.1146/annurev-immunol-041015-055318 [Google Scholar] [PubMed] [CrossRef]

24. Guan, Y., Kraus, S. G., Quaney, M. J., Daniels, M. A., Mitchem, J. B. et al. (2020). FOLFOX chemotherapy ameliorates CD8 T lymphocyte exhaustion and enhances checkpoint blockade efficacy in colorectal cancer. Frontiers in Oncology, 10, 586. https://doi.org/10.3389/fonc.2020.00586 [Google Scholar] [PubMed] [CrossRef]

25. Wang, X., He, Q., Shen, H., Xia, A., Tian, W. et al. (2019). TOX promotes the exhaustion of antitumor CD8+ T cells by preventing PD1 degradation in hepatocellular carcinoma. Journal of Hepatology, 71(4), 731–741. https://doi.org/10.1016/j.jhep.2019.05.015 [Google Scholar] [PubMed] [CrossRef]

26. Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V. et al. (2019). Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nature Immunology, 20(2), 163–172. https://doi.org/10.1038/s41590-018-0276-y [Google Scholar] [PubMed] [CrossRef]

27. Therneau, T. (2013). R survival package. R Core Team. [Google Scholar]

28. Moore, D. F. (2016). Applied survival analysis using R, vol. 473. Germany: Springer. [Google Scholar]

29. Harrell Jr, F. E., Harrell, M. F. E., Hmisc, D. (2017). Package ‘rms’. USA: Vanderbilt University, USA. [Google Scholar]

30. Becht, E., Giraldo, N. A., Lacroix, L., Buttard, B., Elarouci, N. et al. (2016). Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biology, 17(1), 218. https://doi.org/10.1186/s13059-016-1070-5 [Google Scholar] [PubMed] [CrossRef]

31. Aran, D., Hu, Z., Butte, A. J. (2017). xCell: Digitally portraying the tissue cellular heterogeneity landscape. Genome Biology, 18(1), 220. https://doi.org/10.1186/s13059-017-1349-1 [Google Scholar] [PubMed] [CrossRef]

32. Maeser, D., Gruener, R. F., Huang, R. S. (2021). oncoPredict: An R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Briefings in Bioinformatics, 22(6), bbab260. https://doi.org/10.1093/bib/bbab260 [Google Scholar] [PubMed] [CrossRef]

33. Wang, Z., Zhu, J., Liu, Y., Liu, C., Wang, W. et al. (2020). Development and validation of a novel immune-related prognostic model in hepatocellular carcinoma. Journal of Translational Medicine, 18(1), 1–16. https://doi.org/10.1186/s12967-020-02255-6 [Google Scholar] [PubMed] [CrossRef]

34. Sun, Y. L., Zhang, Y., Guo, Y. C., Yang, Z. H., Xu, Y. C. (2020). A prognostic model based on the immune-related genes in colon adenocarcinoma. International Journal of Medical Sciences, 17(13), 1879–1896. https://doi.org/10.7150/ijms.45813 [Google Scholar] [PubMed] [CrossRef]

35. Li, R., Liu, X., Zhou, X. J., Chen, X., Li, J. P. et al. (2020). Identification of a prognostic model based on immune-related genes of lung squamous cell carcinoma. Frontiers in Oncology, 10, 1588. https://doi.org/10.3389/fonc.2020.01588 [Google Scholar] [PubMed] [CrossRef]

36. Yan, C., Niu, Y., Ma, L., Tian, L., Ma, J. (2022). System analysis based on the cuproptosis-related genes identifies LIPT1 as a novel therapy target for liver hepatocellular carcinoma. Journal of Translational Medicine, 20(1), 1–18. https://doi.org/10.1186/s12967-022-03630-1 [Google Scholar] [PubMed] [CrossRef]

37. Sun, W., Xu, Y., Zhao, B., Zhao, M., Chen, J. et al. (2022). The prognostic value and immunological role of angiogenesis-related patterns in colon adenocarcinoma. Frontiers in Oncology, 12, 1003440. https://doi.org/10.3389/fonc.2022.1003440 [Google Scholar] [PubMed] [CrossRef]

38. Ye, F., Hu, Y., Gao, J., Liang, Y., Liu, Y. et al. (2021). Radiogenomics map reveals the landscape of m6A methylation modification pattern in bladder cancer. Frontiers in Immunology, 12, 722642. [Google Scholar] [PubMed]

39. Ye, F., Liang, Y., Hu, J., Hu, Y., Liu, Y. et al. (2021). DNA methylation modification map to predict tumor molecular subtypes and efficacy of immunotherapy in bladder cancer. Frontiers in Cell and Developmental Biology, 9, 760369. [Google Scholar] [PubMed]

40. Papalexi, E., Satija, R. (2018). Single-cell RNA sequencing to explore immune cell heterogeneity. Nature Reviews Immunology, 18(1), 35–45. https://doi.org/10.1038/nri.2017.76 [Google Scholar] [PubMed] [CrossRef]

41. Der, E., Ranabothu, S., Suryawanshi, H., Akat, K. M., Clancy, R. et al. (2017). Single cell RNA sequencing to dissect the molecular heterogeneity in lupus nephritis. JCI Insight, 2(9), e93009. [Google Scholar] [PubMed]

42. Lu, J., Chen, Y., Zhang, X., Guo, J., Xu, K. et al. (2022). A novel prognostic model based on single-cell RNA sequencing data for hepatocellular carcinoma. Cancer Cell International, 22(1), 1–12. https://doi.org/10.1186/s12935-022-02469-2 [Google Scholar] [PubMed] [CrossRef]

43. Fan, F., Zhang, H., Dai, Z., Zhang, Y., Xia, Z. et al. (2021). A comprehensive prognostic signature for glioblastoma patients based on transcriptomics and single cell sequencing. Cellular Oncology, 44(4), 917–935. https://doi.org/10.1007/s13402-021-00612-1 [Google Scholar] [PubMed] [CrossRef]

44. Zhang, X., Qiu, J., Huang, F., Han, P., Shan, K. et al. (2022). Construction and verification of a hypoxia-related nine-gene prognostic model in uveal melanoma based on integrated single-cell and bulk RNA sequencing analyses. Experimental Eye Research, 223(8), 109214. https://doi.org/10.1016/j.exer.2022.109214 [Google Scholar] [PubMed] [CrossRef]

45. Sica, V., Bravo-San Pedro, J. M., Stoll, G., Kroemer, G. (2020). Oxidative phosphorylation as a potential therapeutic target for cancer therapy. International Journal of Cancer, 146(1), 10–17. https://doi.org/10.1002/ijc.32616 [Google Scholar] [PubMed] [CrossRef]

46. Solaini, G., Sgarbi, G., Baracca, A. (2011). Oxidative phosphorylation in cancer cells. Biochimica et Biophysica Acta (BBA)-Bioenergetics, 1807(6), 534–542. https://doi.org/10.1016/j.bbabio.2010.09.003 [Google Scholar] [PubMed] [CrossRef]

47. Ashton, T. M., McKenna, W. G., Kunz-Schughart, L. A., Higgins, G. S. (2018). Oxidative phosphorylation as an emerging target in cancer therapy. Clinical Cancer Research, 24(11), 2482–2490. https://doi.org/10.1158/1078-0432.CCR-17-3070 [Google Scholar] [PubMed] [CrossRef]

48. Löbrich, M., Jeggo, P. A. (2007). The impact of a negligent G2/M checkpoint on genomic instability and cancer induction. Nature Reviews Cancer, 7(11), 861–869. https://doi.org/10.1038/nrc2248 [Google Scholar] [PubMed] [CrossRef]

49. Wang, Y., Ji, P., Liu, J., Broaddus, R. R., Xue, F. et al. (2009). Centrosome-associated regulators of the G2/M checkpoint as targets for cancer therapy. Molecular Cancer, 8(1), 1–13. https://doi.org/10.1186/1476-4598-8-8 [Google Scholar] [PubMed] [CrossRef]

50. Mohamed, A. D., Tremblay, A. M., Murray, G. I., Wackerhage, H. (2015). The Hippo signal transduction pathway in soft tissue sarcomas. Biochimica et Biophysica Acta (BBA)-Reviews on Cancer, 1856(1), 121–129. https://doi.org/10.1016/j.bbcan.2015.05.006 [Google Scholar] [PubMed] [CrossRef]

51. Nishio, M., Hamada, K., Kawahara, K., Sasaki, M., Noguchi, F. et al. (2012). Cancer susceptibility and embryonic lethality in Mob1a/1b double-mutant mice. The Journal of Clinical Investigation, 122(12), 4505–4518. https://doi.org/10.1172/JCI63735 [Google Scholar] [PubMed] [CrossRef]

52. Li, Y., Yang, R., Chen, L., Wu, S. (2020). CD38 as an immunomodulator in cancer. Future Oncology, 16(34), 2853–2861. https://doi.org/10.2217/fon-2020-0401 [Google Scholar] [PubMed] [CrossRef]

53. Zhu, Y., Zhang, Z., Jiang, Z., Liu, Y., Zhou, J. (2020). CD38 predicts favorable prognosis by enhancing immune infiltration and antitumor immunity in the epithelial ovarian cancer microenvironment. Frontiers in Genetics, 11, 369. https://doi.org/10.3389/fgene.2020.00369 [Google Scholar] [PubMed] [CrossRef]

54. Li, W., Cai, S., Wang, L., Yang, C., Zhou, B. et al. (2017). HINT2 downregulation promotes colorectal carcinoma migration and metastasis. Oncotarget, 8(8), 13521–13531. https://doi.org/10.18632/oncotarget.14587 [Google Scholar] [PubMed] [CrossRef]

55. Deng, S., Gottlieb, L., Pan, B., Supplee, J., Wei, X. et al. (2021). Molecular mechanism of N-terminal acetylation by the ternary NatC complex. Structure, 29(10), 1094–1104.e4. https://doi.org/10.1016/j.str.2021.05.003 [Google Scholar] [PubMed] [CrossRef]

56. Xu, B., Gong, X., Zi, L., Li, G., Dong, S. et al. (2019). Silencing of DLEU 2 suppresses pancreatic cancer cell proliferation and invasion by upregulating microRNA-455. Cancer Science, 110(5), 1676–1685. https://doi.org/10.1111/cas.13987 [Google Scholar] [PubMed] [CrossRef]

57. Zhou, Y., Shi, H., Du, Y., Zhao, G., Wang, X. et al. (2019). lncRNA DLEU2 modulates cell proliferation and invasion of non-small cell lung cancer by regulating miR-30c-5p/SOX9 axis. Aging, 11(18), 7386–7401. https://doi.org/10.18632/aging.102226 [Google Scholar] [PubMed] [CrossRef]

58. Zhang, S. F., Wang, X. Y., Fu, Z. Q., Peng, Q. H., Zhang, J. Y. et al. (2015). TXNDC17 promotes paclitaxel resistance via inducing autophagy in ovarian cancer. Autophagy, 11(2), 225–238. https://doi.org/10.1080/15548627.2014.998931 [Google Scholar] [PubMed] [CrossRef]

59. Yuan, A., Hsiao, Y. J., Chen, H. Y., Chen, H. W., Ho, C. C. et al. (2015). Opposite effects of M1 and M2 macrophage subtypes on lung cancer progression. Scientific Reports, 5(1), 1–12. https://doi.org/10.1038/srep14273 [Google Scholar] [PubMed] [CrossRef]

60. Najafi, M., Hashemi Goradel, N., Farhood, B., Salehi, E., Nashtaei, M. S. et al. (2019). Macrophage polarity in cancer: A review. Journal of Cellular Biochemistry, 120(3), 2756–2765. https://doi.org/10.1002/jcb.27646 [Google Scholar] [PubMed] [CrossRef]

61. Hagemann, T., Wilson, J., Burke, F., Kulbe, H., Li, N. F. et al. (2006). Ovarian cancer cells polarize macrophages toward a tumor-associated phenotype. The Journal of Immunology, 176(8), 5023–5032. https://doi.org/10.4049/jimmunol.176.8.5023 [Google Scholar] [PubMed] [CrossRef]

62. Kusuda, T., Shigemasa, K., Arihiro, K., Fujii, T., Nagai, N. et al. (2005). Relative expression levels of Th1 and Th2 cytokine mRNA are independent prognostic factors in patients with ovarian cancer. Oncology Reports, 13(6), 1153–1158. https://doi.org/10.3892/or.13.6.1153 [Google Scholar] [CrossRef]

63. Hur, J., Ghosh, M., Kim, T. H., Park, N., Pandey, K. et al. (2021). Synergism of AZD6738, an ATR inhibitor, in combination with belotecan, a camptothecin analogue, in chemotherapy-resistant ovarian cancer. International Journal of Molecular Sciences, 22(3), 1223. https://doi.org/10.3390/ijms22031223 [Google Scholar] [PubMed] [CrossRef]

64. Colon-Otero, G., Zanfagnin, V., Hou, X., Foster, N. R., Asmus, E. J. et al. (2020). Phase II trial of ribociclib and letrozole in patients with relapsed oestrogen receptor-positive ovarian or endometrial cancers. ESMO Open, 5(5), e000926. https://doi.org/10.1136/esmoopen-2020-000926 [Google Scholar] [PubMed] [CrossRef]

65. Zhang, Z., Chen, L., Chen, H., Zhao, J., Li, K. et al. (2022). Pan-cancer landscape of T-cell exhaustion heterogeneity within the tumor microenvironment revealed a progressive roadmap of hierarchical dysfunction associated with prognosis and therapeutic efficacy. eBioMedicine, 83(2), 104207. https://doi.org/10.1016/j.ebiom.2022.104207 [Google Scholar] [PubMed] [CrossRef]

Supplementary Materials