T cell exhaustion-related exosome genes for predicting survival and immunotherapy efficacy in colorectal cancer

doi:10.21203/rs.3.rs-4933597/v1

Download PDF

Research Article

T cell exhaustion-related exosome genes for predicting survival and immunotherapy efficacy in colorectal cancer

https://doi.org/10.21203/rs.3.rs-4933597/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Treatment options for colorectal cancer are limited. T cell exhaustion is one of the barriers to tumor immunotherapy. No comprehensive analysis of T cell exhaustion-related exosome prognostic models for colorectal cancer (CRC) has been conducted.

Method

Samples were collected from the Cancer Genome Atlas (TCGA) database, exoRBase database and Gene Expression Omnibus (GEO) database. The single sample gene set enrichment analysis (ssGSEA) algorithm screened out T cell exhaustion-related exosome differential expression genes, signature genes were screened by univariate Cox regression and Lasso regression, and risk score models were constructed and validated. A nomogram containing risk scores and clinical parameters was established and evaluated. In addition, single cell analysis and tumor immune microenvironment assessment were also performed.

Results

Sixteen signature genes were identified, based on which the risk score model was constructed and validated. This model can predict the overall survival (OS) of TCGA and GEO queues well. Scores were identified as independent risk factors for OS and correlated with certain clinicopathological features. A nomogram was developed that integrated clinical parameters and risk scores and showed higher predictive accuracy. Finally, significant differences in immune microenvironment were found between the high- and low-risk groups. Thus, scores can also be used to predict the response to immunotherapy.

Conclusions

In general, we screened out T cell exhaustion-related exosome genes of CRC, constructed a risk score model which could predict survival and immunotherapy efficacy, and found correlations between risk scores and clinicopathologic features and immune microenvironment.

T cell exhaustion

exosome

colorectal cancer

tumor immune microenvironment

bioinformatics

Colorectal cancer (CRC) has the third highest incidence and mortality in the world nowadays, and is characterized by high heterogeneity and invasiveness, posing a serious threat to human health [1]. New prevention and treatment strategies are urgently needed. Clinical decisions on chemotherapy for CRC depend primarily on clinicopathological staging, regardless of molecular biological characteristics, and this inadequate decision-making approach may lead to potential over- or under-treatment. In the era of individualized therapy, it is imperative to identify reliable biomarkers to optimize the prognosis of colorectal cancer pharmacotherapy.

Exosomes are tiny goblet vesicles secreted by cells. A growing number of studies have clarified the correlation between exosome production and tumorigenesis [2–4]. Tumor-derived exosomes participate in the exchange of genetic information between tumor cells and basal cells, affecting angiogenesis and promoting tumor growth, metastasis and invasion [5]. Numerous pieces of evidence strongly indicate that exosomes play a pivotal regulatory role in the Tumor Microenvironment (TME) of CRC, influencing factors such as flora, hypoxia, inflammation, and the immunological microenvironment [6]. As tumor markers for diagnosis and treatment, exosomes have become one of the focuses of current research. Exosome circSATB2 is associated with non-small cell lung cancer progression and therefore can potentially be used as a diagnostic marker for this cancer type [7]. Exosome circ0048117 regulates esophageal squamous cell carcinoma progression, and higher serum exosome circ0048117 is significantly and positively correlated with TNM stage [8]. Furthermore, circulating exosome miR-203 levels correlate with metastasis, and low miR-203 expression in tumor tissues is a poor prognostic factor in colorectal cancer [9].

In chronic infections or cancer, T cells gradually lose effector function and memory T cells begin to be depleted due to prolonged exposure to antigens and inflammation, this process is known as T Cell exhaustion [10], which is one of the major obstacles to anti-cancer immunotherapy. T cell exhaustion represents a unique state of T cell differentiation [11], has considerable clinical relevance, but our understanding of it is still incomplete.

To date, no comprehensive analysis of T cell exhaustion-related exosome prognostic models for colorectal cancer has been conducted. The biological process of T cell exhaustion-related exosomes in colorectal cancer and the relationship between them and clinical treatment and immunotherapy of patients are still unclear. Therefore, we systematically screened out the signature genes of T cell exhaustion-related exosomes in colorectal cancer based on the existing database data, and analyzed the correlation between the signature genes and clinical treatment and immunotherapy of patients.

2.1 Data acquisition and preprocessing

Clinical information and expression data for colonic adenocarcinoma (COAD), rectum adenocarcinoma (READ) were obtained from The Cancer Genome Atlas of America database (TCGA, https://cancergenome.nih.gov/), and were combined as a training set after removing the batch effects. Colorectal cancer exosome and normal exosome expression profiles were downloaded from the exoRBase database (http://www.exorbase.org/). The set GSE41258 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE41258) was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) as a validation set, and the set GSE132465 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132465) as a single cell quality control data set. T cell exhaustion related genes were extracted from the study of Zhang et al [12] and shown in Supplementary table 1.

2.2 Screening of differentially expressed genes (DEGs)

DEGs were filtered through the “limma ” [13] package. The DEGs threshold point was an adjusted p value, corrected by the BH method, less than 0.05 and a |log₂ fold change| greater than 0.585.

2.3 Single sample gene set enrichment analysis (ssGSEA)

Enrichment analysis for the expression profile of TCGA-COADREAD tumor samples was performed by R package “GSVA” [14] based on 40 T cell exhaustion-related genes obtained.

2.4 Screening of TEX-exo genes

The intersection of the two groups of DEGs was named “exo-genes”. Spearman correlation between exo-genes expression and T cell exhaustion-related enrichment score was calculated in the training set, and the significantly related genes were retained as TEX-exo genes.

2.5 Functional enrichment analysis

Gene ontology (GO), Kyoto Encyclopaedia of Genes and Genomes (KEGG) analyses and top10 visual display of the TEX-exo genes were performed using the R package “clusterProfiler” [15], the parameters were pvalueCutoff = 0.05, pAdjustMethod = “BH”, and qvalueCutoff = 0.5.

2.6 Construction of the Signature model

The hazard ratio (HR) and prognostic significance of TEX-exo genes were determined by univariate Cox regression analysis. The genes with p < 0.05 were prognostic related genes. The LASSO algorithm was implemented through the R package “glmnet” [16] for further selection and contraction of prognostic related genes. The selected genes were the signature genes. Then calculated the risk scores:

Score = \(\:{\sum\:}_{i=0}^{n}\beta\:i\text{*}\chi\:i\)

n for the number of signature genes, 𝛽\(\:i\) for the weight coefficient of each gene, 𝜒\(\:i\) for the expression of each gene.

Patients were divided into high - and low-risk groups based on the median Score, the “survminer” package was used to analyze the OS and the "timeROC" package was used to plot a time-dependent receiver operator characteristic (ROC) curve. Finally, univariate and multivariate Cox analyses were performed to explore the independent prognostic value of risk Scores and to validate them in the validation set.

2.7 Nomogram construction

A nomogram including clinical features and risk Scores was constructed using the R package “rms”, and a calibration curve was then drawn to assess the accuracy of the nomogram in predicting survival rate.

2.8 Single-cell analysis

R package "Seurat" [17] was used for quality control of colorectal cancer samples to exclude low-quality cells and low-expression genes, and the Score of each cell was calculated according to the formula above.

2.9 Evaluation of tumor immune microenvironment (TIM)

CIBERSORTx algorithm in R package “IOBR” [18] was used to calculate the Scores of 22 immune cells in tumor immune microenvironment in the training set.

2.10 Statistical analysis

All analyses in this paper were carried out by R version 4.1.2. The Wilcoxon rank sum test was used to analyze the significance of two groups of numerical variables, and Kruskal − Wallis test for more than two groups. In the plots display, ns means p > 0.05, * means p < 0.05, ** means p < 0.01, *** means p < 0.001, **** means p < 0.0001. Survival curves were generated by Kaplan-Meier method, and the significance of the differences was determined by log-rank test.

3.1 Screening of TEX-exo genes

The DEGs between normal exosome and colorectal cancer exosome samples was calculated by the “limma” algorithm, and 1016 DEGs were obtained (Fig. 1a). The expression difference of normal and tumor samples in TCGA-COADREAD was calculated by the same method, and 4352 DEGs were obtained (Fig. 1b). Taking the intersection of the DEGs of the above two resulted in 276 differential exo-genes.

The T cell exhaustion-related enrichment scores for each sample in TCGA-COADREAD are generated based on 40 T cell exhaustion-related genes with the "GSVA" algorithm, as detailed in Supplementary Table 2. Spearman correlation was calculated between the expression of 276 exo-genes and T cell exhaust-related enrichment score in TCGA-COADREAD, and 179 genes with significant correlation (p < 0.05) were retained as the TEX-exo genes for further analysis. Figure 1c shows 8 genes with the highest correlation among the genes with significant positive correlation, including DOCK2, LSP1, HCLS1, NCKAP1L, GIMAP1, ARHGEF6, GYPC and ARHGAP25.

3.2 Functional enrichment analysis

To understand the biological mechanisms related to T cell exhaustion in exosomes, KEGG and GO enrichment of TEX-exo genes were analyzed, and the results are shown in Fig. 2a. The biological process of GO enrichment showed that these TEX-exo genes were mainly distributed in the lymphocyte proliferation, mononuclear cell proliferation, leukocyte proliferation, and so forth. The celluar component of GO enrichment suggested that these TEX-exo genes were involved in the cell cortex, secrebory granule lumen, spliceosomal snRhP complex and so on. The molecular function of GO enrichment indicated that these TEX-exo genes enriched in the carbohydrate binding, MHC protein complex binding and other functions. In addition, the KEGG results demonstrated that these TEX-exo genes enriched in Hematopoietic cell lineage, Intestinal immune network for IgA production and other pathways.

To further confirm the function of proteins encoded by TEX-exo genes, we download COAD and READ level 4 protein expression data from the TCGA database (https://www.tcpaportal.org/tcpa/). Among the 179 TEX-exo genes, only the MS4A1 gene was able to correspond to its matching encoded protein, CD20. Therefore, We further compared the differences of CD20 in six different clinical feature groups (age, sex, stage, colon polyps, lymphatic invasion and venous invasion) of COAD and READ. There were no significant differences in all six groups. The results can be seen in Supplementary Fig. 1.

To gain further insight into the genetic variation of the TEX-exo genes, we also analyzed the incidence of copy number variation (CNV) and somatic mutations in the 179 TEX-exo genes. First of all, Fig. 2b showed the genes with single nucleotide variation (SNV) frequency above 3%. Among them, DOCK2 exhibited the highest mutation frequency, followed by AKAP12. Then, investigation of the frequency of CNV revealed that CNV was widespread in the 179 TEX-exo genes, for example, there were extensive amplification in IGF2 and MMP9, and extensive deletion in TCEA3 and STMN1, as shown in Fig. 2c. This systematically revealed genetic variation in TEX-exo genes, illustrating the potential pathogenic landscape of these genes in CRC.

3.3 Building the signature model

Univariate Cox regression analysis was used to evaluate the effect of TEX-exo genes on OS. Nineteen genes were found to be significant prognostic factors (p < 0.05) (Fig. 3a). Of these genes, seven (HR > 1; RBPMS2, AKAP12, SCARA3, TIMP1, HSPA1A, RBMX2 and C1orf35) were risk factors for OS and 12 (HR < 1; PPP1R16B, MS4A1, POU2AF1, BIRC3, FKBP5, DNASE1L3, CTNND1, FAM177B, SULT1B1, CEP70, UQCRFS1 and PSRC1) were protective factors for OS. Survival analysis was performed on these 19 prognostic factors in tumor samples, and the median expression level was used as the cut-off value to divide the high and low expression groups to compare the influence of expression level on prognosis. Figure 3b showed the survival curves of 8 genes with the most significant p value (p < 0.05). There are 4 genes with low expression that have a better prognosis (HR < 1; HSAP1A, RBMX2, C1orf35 and SCARA3). On the contrary, 4 genes with high expression have a superior prognosis (HR > 1; UQCRFS1, MS4A1, FAM177B and CEP70).

Given the prognostic value of the 19 genes, we would like to build a more reliable and useful predictive model from these genes to guide clinical practice in colorectal cancer patients. The expression profiles of the 19 genes above were analyzed by LASSO-Cox regression analysis, and the prognostic model was constructed. 16 signature genes were identified based on the optimal value of λ 0.0087, as shown in Fig. 3c-e. Furthermore, the regression coefficients of each signature gene were shown in Table 1. The results of univariate COX analysis of these 16 signature genes in OS are displayed in the forest plot of Fig. 3f, and all of their expressions were strongly associated with prognosis.

Table.1 Regression coefficients for 16 signature gene.

Symbol	Coef
C1orf35	0.619111702
TIMP1	0.281871688
CTNND1	0.215980486
HSPA1A	0.170078448
RBMX2	0.1222027
CEP70	0.088315003
MS4A1	0.03913604
POU2AF1	0.032869455
RBPMS2	0.004424399
AKAP12	-0.010104709
DNASE1L3	-0.0677533
SULT1B1	-0.08324181
FAM177B	-0.089640959
UQCRFS1	-0.151520918
FKBP5	-0.206617298
PSRC1	-0.401296281

Based on the coefficients calculated for each gene by the LASSO-COX algorithm and the expression values of the 16 signature genes for each CRC patient, we were able to compute the Score for each patient according to the formula in the Methods section. Patients were divided into two groups with high Score and low Score according to the median value (Fig. 4a). Survival times for each high Score patient and low Score patient were displayed in Fig. 4b. Overall, patients with high Score had the shorter survival time and those with low Score had the longer survival time. Moreover, the expression of the 16 signature genes in patients with high and low Score was exhibited in Fig. 4c. Kaplan-Meier curve showed that the survival probability of patients with higher Scores was significantly lower than that of patients with lower Scores (HR = 2.97; 95% CI = 1.98–4.47; P < 0.001; Fig. 4d). The area under ROC curve (AUC) was 0.712 in year 1, 0.676 in year 3, and 0.721 in year 5, which means that this risk scoring model has a high value for prognostic guidance (Fig. 4e).

In order to verify the stability of the model, the same algorithm was used to analyze the verification set (GSE41258). Similarly, the validation set patients were categorized into high and low Score based on the median Score (Fig. 4f), whereas the survival time of the high Score patients remained longer overall, while the survival time of the low Score patients continued to be shorter on the whole (Fig. 4g). In addition, the expression of 16 signature genes in high and low Score patients was demonstrated in Fig. 4h. Consistent with results obtained in the training set TCGA-COADREAD cohort, patients in the higher Score group had a shorter survival time compared to those in the lower Score group (HR = 1.75; 95% CI = 1.18–2.61; P = 0.005; Fig. 4i). In addition, in terms of survival effect prediction, the 1-year AUC was 0.724, 3-year AUC was 0.640, and 5-year AUC was 0.603, which reaffirmed the excellent efficacy of this risk score model in predicting prognosis (Fig. 4j).

In order to exclude the influence of other factors, univariate and multivariate Cox analyses were used to determine whether Score was an independent prognostic factor for OS. In univariate Cox analysis, Score obtained by the training set TCGA-COADREAD cohort was significantly correlated with OS (High vs. Low; HR = 2.97, 95% CI = 1.98–4.47, P < 0.001; Fig. 5a). After adjusting for other confounding factors, multivariate Cox analysis showed that Score was still an independent predictor of OS (High vs. Low; HR = 3.47, 95% CI = 1.74–6.93, P < 0.001), as shown in Fig. 5b.

Similarly, univariate and multivariate Cox analysis was used in the verification set. In univariate Cox analysis, Score obtained by the validation set cohort was significantly correlated with OS (High vs. Low; HR = 1.75, 95% CI = 1.18–2.61, P = 0.006; Fig. 5c). After adjusting for other confounding factors, multivariate Cox analysis showed that Score was still an independent predictor of OS (High vs. Low; HR = 1.52, 95% CI = 1.01–2.29, P = 0.043), as shown in Fig. 5d.

3.4 Nomogram construction

We combined Score with five different clinical features (age, sex, stage, lymphatic invasion and venous invasion) to construct a nomogram to predict the probability of 1-year, 3-year, and 5-year OS. Each factor was assigned in proportion to its risk contribution to survival (Fig. 6a), and calibration curves showed that the combined model (nomogram) showed high accuracy over 1 -, 3 -, and 5-year OS (Fig. 6b). Compared with a single prognostic factor, the nomogram constructed using a combination model may be a better predictor of patients' OS (Fig. 6c). We analyzed the robustness of Score in different clinical features, and the results were shown in Fig. 6d. In most clinical groups, we observed that the prognosis of the patients with low Scores was better than that of the patients with high Scores (HR < 1, P < 0.05). We also compared the differences in Scores among six different clinical feature groups, as shown in Fig. 6e. Scores were higher in stage III/IV patients than in stage I/II patients, in patients with lymphatic invasion than in patients without lymphatic invasion, and in patients with venous invasion than in patients without venous invasion. Additionally, there were no significant differences in Scores among different age, sex and colon polyps.

3.5 Single-cell analysis

The role and value of the Score was further validated using the single-cell analysis of colorectal cancer. The single-cell analysis of 23 colorectal cancer tumor samples showed that 40,000 cells were retained from the original 47,285 cells after quality control. Figure 7a and Fig. 7b showed the landscape of cell distribution before batch effect and after debatching effect, respectively. The distribution of the six cell types (B cells, Epithelial cells, Mast cells, Myeloids, Stromal cells and T cells) can be derived from the gene expression values of each cell (Fig. 7c). The same formula was used to calculate the Score of each cell, and the results were displayed by using UMAP, as shown in Fig. 7d. Among the six cell types, Stromal cells had the highest Scores, while T cells had the lowest (Fig. 7e). The differences in Scores between the other five cells and the Epithelial cells were shown in Fig. 7f. Among them, compared with Epithelial cells, the Scores of Mast cells, Myeloids and Stromal cells were relatively high (log2FC > 0), while the Scores of B cells and T cells were relatively low (log2FC < 0).

3.6 Evaluation of TME

We further explored the tumor immune microenvironment between the high and low Scores groups with the aim of identifying the underlying immune mechanisms. Different immune cell subsets were quantified using CIBERSORTx, and rank sum tests were used to compare the significance of infiltration degree between groups. We found that the infiltration degree of B cells memory, Dendritic cells activated Macrophages M0 cells, Plasma cells, T cells CD4 memory resting, and T cells regulatory were significantly different in higer and lower Score groups. The results are shown in Fig. 8a.

We also investigated the differences in the expression of immune checkpoints in the higher and lower Score groups. Expression of all active, inhibit and two-side immune checkpoint genes differed significantly between high and low Scores (Fig. 8b and Supplementary Fig. 2).

We also explored whether prognostic signature could predict a patient's response to immune checkpoint blocking therapy. Immunotherapy data for a metastatic melanoma was used [19], it was found that patients with high Scores have a better prognosis than those with low Scores (P = 0.038; Fig. 8c). Furthermore, patients who responded to immunotherapy drugs had higher Scores than those who did not respond to the drugs (P = 0.026; Fig. 8d). Among the high Scores patients, 28.57% responded to immunotherapy, while 71.43% did not. In addition, all patients in the low Scores patients did not respond to immunotherapy (Fig. 8e). Finally, the high Scores group had higher TIDE scores than the low Scores group in TCGA-COADREAD (P = 3.9e-11; Fig. 8f). In conclusion, the immunotherapy effect of the high Scores group was better than that of the low Scores group. Therefore, Scores can be used as a novel marker to guide immunotherapy in colorectal cancer patients.

T cell exhaustion is critical in tumor immunotherapy [20]. Studies show that severe exhaustion of tumor-infiltrating T cells in microsatellite stabilized (MSS) colorectal cancer is one of the important mechanisms by which patients develop resistance to PD-1 inhibitors [21]. Tumor-derived exosomes (TDEs) are involved in various processes of cancer formation and development, including tumor microenvironment remodeling, angiogenesis, invasion, metastasis, and drug resistance [22]. It has been found that TDEs can promote the occurrence of liver metastasis in CRC by regulating the crosstalk between tumor cells and macrophages [23]. To our knowledge, this study is the first to systematically evaluate T cell exhaustion-related exosome genes in CRC. Firstly, a risk score model containing 16 signature genes was established through the TCGA cohort, and a nomogram was constructed to predict the OS of CRC patients, while the model was validated using an external data set. Secondly, the GSE132465 data set was used to characterize the features of the scores of signature genes in single cells. Finally, we also found that signature gene scores were associated with other clinicopathologic features and tumor immune microenvironment in patients with colorectal cancer.

Previous studies have shown that non-coding RNA (ncRNA) and circular RNA (circRNA) derived from exosome can affect the tumor immune microenvironment, including inducing T cell exhaustion and promoting the overexpression of programmed death ligand-1 (PD-L1) [24–26]. However, considering that the cross-talk between mRNA in exosomes and T cell exhaustion on CRC needs to be explored, we identified 179 TEX-exo genes, which are not only DEGs of CRC exosomes versus normal exosomes, but also DEGs closely related to T cell exhaustion. Eight exosomal genes, including DOCK2, LSP1, HCLS1, NCKAP1L, GIMAP1, ARHGEF6, GYPC and ARHGAP25, were significantly positively correlated with T cell exhaustion (R > 0.8, p < 0.05). It has been reported that cholesterol sulfate synthesized by SULT2B1 in hepatocellular carcinoma inhibit DOCK2 activity in T cells and promote effector T cell exhaustion [27]. In addition, a combination of RNA-seq and proteomics analysis of human NCKAP1L deficiency cases reported by Castro et al. [28] showed characteristics of T cell exhaustion. 179 TEX-exo genes were mainly concentrated in immune-related molecules or signaling pathways, such as MHC class II protein complex, lymphocyte proliferation and Intestinal immune network for IgA production. Kilian et al. [29] showed that MHC Class II restricted antigen presentation is required to prevent cytotoxic T cell dysfunction in brain tumors. In addition, despite the prevalence of CNV in 179 TEX-exo genes, the relationship between T cell exhaustion and CNV remains unclear, further studies are needed.

We constructed a signature gene scoring model based on 16 T cell exhaustion-related exosome genes. Signature genes scores in both the TCGA cohort and GSE41258 showed excellent prognostic ability and were identified as an independent risk factors for OS (P < 0.05). Compared with a single risk score model, a nomogram with signature genes scores and clinicopathological features is able to predict patients' OS more accurately. In addition, the current signature genes risk score is associated with some clinicopathologic features, including TNM staging, lymph node invasion, and vascular invasion.

Among the 16 signature genes, RBPMS2, AKAP12, TIPM1, HSPA1A, RBMX2 and C1orf35 were negatively correlated with CRC patients' OS (HR > 1), while MS4A1, POU2AF1, FKBP5, DNASE1L3, CTNND1, FAM177B, SULT1B1, CEP70, UQCRFS1 and PSRC1 were positively correlated with OS in patients with colorectal cancer (HR < 1). He et al. [30] reported that down-regulation of AKAP12 can inhibit the progression and migration of CRC through the PI3K/AKT signaling pathway. It was reported that TIMP1 expression was significantly associated with regional lymph nodes and distant metastases [31]. The down-regulation of HSPA1A inhibits the proliferation and migration of CRC cells, and CRC patients with lower expression levels tend to have longer OS [32, 33]. MS4A1, DNASE1L3 and SULT1B1 were significantly down-regulated in CRC tissues, and their reducing levels were associated with shorter OS in CRC patients [34–36]. UQCRB plays a key role in mitochondrial complex III stability, electron transport, cellular oxygen sensing and angiogenesis, and can be used as an important diagnostic marker for CRC [37]. Although our study showed a positive correlation between CTDNN1 and OS, some basic medicine studies suggested that the activation of CTDNN1 could induce the proliferation and migration of CRC cells, so the experimental results should be treated with caution [38]. In addition, our study found that RBPMS2, RBMX2, C1orf35, POU2AF1, FKBP5, FAM177B, CEP70, and PSRC1 are potential prognotic markers for CRC, however their mechanisms of action have not been investigated.

T cells are critical to the efficacy of current tumor immunotherapy [20]. T cell exhaustion is a state of function diminishing characterized by a progressive loss of T cell effector function and self-renewal [39], which hampers immunotherapy of tumors. Blocking immune checkpoints can eliminate tumors by restoring immunity, thereby restoring dysfunctional/exhausted T cells [40]. In this study, a melanoma cohort receiving immunotherapy was validated [19], we found that the signature gene score model predicted the efficacy of the immune checkpoint inhibitor (ICIs), patients with a higher score having a better prognosis and better efficacy. This may be closely related to immune cells infiltrating and immune checkpoints. These findings indicate that the signature gene score model in this study is promising to be a new indicator for evaluating tumor immune microenvironment and ICIs efficacy. Nevertheless, the association between the signature gene score model and immunotherapy efficacy needs to be further validated in larger samples and the underlying mechanism needs to be further explored.

Inevitably, there are some limitations in our research. Firstly, further experiments in vivo and mechanistic studies are needed to reveal the exact role of each signature gene. Secondly, the potential of the model to predict immunotherapy response was assessed only indirectly, as no mRNA expression data from CRC patients receiving immunotherapy was searched. Thirdly, the predictive power of the model needs to be evaluated based on external validation from prospective and large-scale clinical trials.

In summary, we identified 16 T cell exhaustion-related exosome genes that may play an important role in the development and progression of CRC. The risk score model constructed based on these genes reflects the unique clinicopathological characteristics and immune microenvironment characteristics of CRC patients, and may promote the application of precision medicine in CRC by predicting prognosis and guiding immunotherapy.

CRC

Colorectal Cancer

TME

Tumor Microenvironment

COAD

Colonic Adenocarcinoma

READ

Rectum Adenocarcinoma

TCGA

The Cancer Genome Atlas

GEO

Gene Expression Omnibus

DEGs

Differentially Expressed Genes

ssGSEA

single sample Gene Set Enrichment Analysis

Gene Ontology

KEGG

Kyoto Encyclopaedia of Genes and Genomes

hazard ratio

ROC

Receiver Operator Characteristic

TIM

Tumor Immune Microenvironment

CNV

Copy Number Variation

SNV

Single Nucleotide Variation

TMB

Tumor Mutation Burden

AUC

Area Under ROC Curve

Overall Survival

MSS

Microsatellite Stabilized

TDEs

Tumor-derived Exosomes

ncRNA

non-coding RNA

circRNA

circular RNA

PD-L1

Programmed Death Ligand-1

ICIs

Immune Checkpoint Inhibitor.

Ethics approval and consent to participate

The study was approved by the ethical committee of The First People’s Hospital of Foshan.

Consent for publication

Not applicable.

Availability of data and materials

All data generated or analysed during this study are included in this published article and its supplementary information files.

Competing Interests

The authors declare that there is no conflict of interest regarding the publication of this paper.

Funding

This research was funded by the 14th Five-Year Medical high level key specialty construction project of Foshan (FSGSP145001), The 2023 Foshan Municipal Science and Technology Bureau's Self-Funded Scientific and Technological Innovation Project (2320001006343).

Author Contributions

Yilin Wang: Conceptualization, Software, Writing-original draft. Peizhu Su: Visualization, Writing - review & editing. Qinghua Lu: Formal analysis. Huiwen Huang: Data curation. Zhaotao Li: Project administration, Funding acquisition. All authors read and approved the final manuscript.

Acknowledgements

Not applicable.

Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics. CA Cancer J Clin. 2021; 71:7-33.
Hessvik NP, Llorente A. Current knowledge on exosome biogenesis and release. Cell Mol Life Sci. 2018;75:193-208.
Li X, Corbett AL, Taatizadeh E, Tasnim N, Little JP, Garnis C, et al. Challenges and opportunities in exosome research-Perspectives from biology, engineering, and cancer therapy. APL bioengineering 2019;3:011503.
Sedgwick AE, D'Souza-Schorey C. The biology of extracellular microvesicles. Traffic. 2018;19:319-27.
Ma J, Wang P, Huang L, Qiao J, Li J. Bioinformatic analysis reveals an exosomal miRNA-mRNA network in colorectal cancer. BMC medical genomics. 2021;14:60.
Yao J, Chen Y, Lin Z. Exosomes: Mediators in microenvironment of colorectal cancer. Int J Cancer. 2023;153:904-917.
Zhang N, Nan A, Chen L, Li X, JiaY, Qiu M, et al. Circular RNA circSATB2 promotes progression of non-small cell lung cancer cells. Mol Cancer. 2020;19:101.
Lu Q, Wang X, Zhu J, Fei X, Chen H, Li C. Hypoxic Tumor-Derived Exosomal Circ0048117 Facilitates M2 Macrophage Polarization Acting as miR-140 Sponge in Esophageal Squamous Cell Carcinoma. Onco Targets Ther. 2020;18:13:11883-11897.
Takano Y, Masuda T, Iinuma H, Yamaguchi R, Sato K, Tobo T, et al. Circulating exosomal microRNA-203 is associated with metastasis possibly via inducing tumor-associated macrophages in colorectal cancer. Oncotarget. 2017;8:78598-78613.
Wherry EJ, Kurachi M. Molecular and cellular insights into T cell exhaustion. Nature reviews Immunology. 2015;15:486-99.
Franco F, Jaccard A, Romero P, Yu YR, Ho PC. Metabolic and epigenetic regulation of T-cell exhaustion. Nature metabolism. 2020;2:1001-12.
Zhang Z, Chen L, Chen H, Zhao J, Li K, Sun J, et al. 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. 2022;83:104207.
Smyth GK. limma: Linear Models for Microarray Data. Springer New York. 2005.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC bioinformatics. 2013;14:7.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics-a Journal of Integrative Biology. 2012;16:284-7.
Robert, Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society Series B (Methodological). 1996.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature biotechnology. 2018;36:411-20.
Zeng D, Ye Z, Shen R, Yu G, Wu J, Xiong Y, et al. IOBR: Multi-Omics Immuno-Oncology Biological Research to Decode Tumor Microenvironment and Signatures. Frontiers in immunology. 2021;12:687975.
Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nature medicine. 2018;24:1545-9.
Thommen DS, Schumacher TN. T Cell Dysfunction in Cancer. Cancer cell. 2018;33:547-62.
Kim CG, Jang M, Kim Y, Leem G, Kim KH, Lee H, et al. VEGF-A drives TOX-dependent T cell exhaustion in anti-PD-1-resistant microsatellite stable colorectal cancers. Science immunology. 2019;4:eaay0555.
Mashouri L, Yousefi H, Aref AR, Ahadi AM, Molaei F, Alahari SK. Exosomes: composition, biogenesis, and mechanisms in cancer metastasis and drug resistance. Molecular cancer. 2019;18:75.
Zhao S, Mi Y, Guan B, Zheng B, Wei P, Gu Y, et al. Tumor-derived exosomal miR-934 induces macrophage M2 polarization to promote liver metastasis of colorectal cancer. Journal of hematology & oncology. 2020;13:156.
Xu Z, Chen Y, Ma L, Chen Y, Liu J, Guo Y, et al. Role of exosomal non-coding RNAs from tumor cells and tumor-associated macrophages in the tumor microenvironment. Molecular therapy : the journal of the American Society of Gene Therapy. 2022;30:3133-54.
Yang C, Wu S, Mou Z, Zhou Q, Dai X, Ou Y, et al. Exosome-derived circTRPS1 promotes malignant phenotype and CD8+ T cell exhaustion in bladder cancer microenvironments. Molecular therapy : the journal of the American Society of Gene Therapy. 2022;30:1054-70.
Wang J, Zhao X, Wang Y, Ren F, Sun D, Yan Y, et al. circRNA-002178 act as a ceRNA to promote PDL1/PD1 expression in lung adenocarcinoma. Cell death & disease. 2020;11(1):32.
Wang S, Wang R, Xu N, Wei X, Yang Y, Lian Z, et al. SULT2B1-CS-DOCK2 axis regulates effector T cell exhaustion in hepatocellular carcinoma microenvironment. Hepatology. 2023;78:1064-1078.
Castro CN, Rosenzwajg M, Carapito R, Shahrooei M, Konantz M, Khan A, et al. NCKAP1L defects lead to a novel syndrome combining immunodeficiency, lymphoproliferation, and hyperinflammation. The Journal of experimental medicine. 2020;217:e20192275.
Kilian M, Sheinin R, Tan CL, Friedrich M, Krämer C, Kaminitz A, et al. MHC class II-restricted antigen presentation is required to prevent dysfunction of cytotoxic T cells by blood-borne myeloids in brain tumors. Cancer cell. 2023;41:235-251.
He P, Li K, Li SB, Hu TT, Guan M, Sun FY, et al. Upregulation of AKAP12 with HDAC3 depletion suppresses the progression and migration of colorectal cancer. International journal of oncology. 2018;52:1305-16.
Song G, Xu S, Zhang H, Wang Y, Xiao C, Jiang T, et al. TIMP1 is a prognostic marker for the progression and metastasis of colon cancer through FAK-PI3K/AKT and MAPK pathway. Journal of experimental & clinical cancer research. 2016;35(1):148.
Ding Q, Hou Z, Zhao Z, Chen Y, Zhao L, Xiang Y. Identification of the prognostic signature based on genomic instability-related alternative splicing in colorectal cancer and its regulatory network. Frontiers in bioengineering and biotechnology. 2022;10:841034.
Xing XL, Yao ZY, Xing C, Huang Z, Peng J, Liu YW. Gene expression and DNA methylation analyses suggest that two immune related genes are prognostic factors of colorectal cancer. BMC medical genomics. 2021;14:116.
Mudd TW, Jr., Lu C, Klement JD, Liu K. MS4A1 expression and function in T cells in the colorectal cancer tumor microenvironment. Cellular immunology; 2021;360:104260.
Liu J, Yi J, Zhang Z, Cao D, Li L, Yao Y. Deoxyribonuclease 1-like 3 may be a potential prognostic biomarker associated with immune infiltration in colon cancer. Aging. 2021;13:16513-26.
Lian W, Jin H, Cao J, Zhang X, Zhu T, Zhao S, et al. Identification of novel biomarkers affecting the metastasis of colorectal cancer through bioinformatics analysis and validation through qRT-PCR. Cancer cell international. 2020;20:105.
Kim HC, Chang J, Lee HS, Kwon HJ. Mitochondrial UQCRB as a new molecular prognostic biomarker of human colorectal cancer. Experimental & molecular medicine. 2017;49:e391.
Liu D, Zhang H, Cui M, Chen C, Feng Y. Hsa-miR-425-5p promotes tumor growth and metastasis by activating the CTNND1-mediated β-catenin pathway and EMT in colorectal cancer. Cell cycle. 2020;19:1917-27.
Chow A, Perica K, Klebanoff CA, Wolchok JD. Clinical implications of T cell exhaustion for cancer immunotherapy. Nature reviews Clinical oncology. 2022;19:775-90.
Tsai HF, Hsu PN. Cancer immunotherapy by targeting immune checkpoints: mechanism of T cell dysfunction in cancer immunity and new therapeutic targets. Journal of biomedical science. 2017;24:35.

No competing interests reported.

SupplementaryFigure1.png
Supplementary Figure 1 The differences of CD20 in COAD and READ. (a) The differential expression of CD20 in COAD patients with different clinical features. (b) The differential expression of CD20 in READ patients with different clinical features.
SupplementaryFigure2.tif
Supplementary Figure 2 Differences in active, inhibit and two-side immune checkpoint genes expression between high and low Score groups.
Supplementarytable1.txt
Supplementarytable2.txt

Download PDF

Editor assigned by journal
26 Aug, 2024
Submission checks completed at journal
19 Aug, 2024
First submitted to journal
18 Aug, 2024

You are reading this latest preprint version

T cell exhaustion-related exosome genes for predicting survival and immunotherapy efficacy in colorectal cancer

Status:

Version 1

Abstract

Background

Method

Results

Conclusions

Figures

1. Introduction

2. Materials and Methods

2.1 Data acquisition and preprocessing

2.2 Screening of differentially expressed genes (DEGs)

2.3 Single sample gene set enrichment analysis (ssGSEA)

2.4 Screening of TEX-exo genes

2.5 Functional enrichment analysis

2.6 Construction of the Signature model

2.7 Nomogram construction

2.8 Single-cell analysis

2.9 Evaluation of tumor immune microenvironment (TIM)

2.10 Statistical analysis

3. Results

3.1 Screening of TEX-exo genes

3.2 Functional enrichment analysis

3.3 Building the signature model

3.4 Nomogram construction

3.5 Single-cell analysis

3.6 Evaluation of TME

4. Discussion

5. Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1