Identification of Tregs-specific mRNAs associated with RFS of prostate cancer.
To identify Tregs-specific mRNAs, differential expression analysis of mRNAs was performed between Tregs and other immune cell lines, and 2407 dysregulated mRNAs were identified that were highly expressed in Tregs and downregulated in other immune cell lines (false discovery rate (FDR) < 0.05, LogFC > 2, Table S3). These 2407 dysregulated mRNAs were proposed as Tregs-specific mRNAs. Ngar-Yee Huen et, al. identified 384 genes that were differentially expressed between healthy donors and prostate cancer patients using gene expression profiles excavated from peripheral blood mononuclear cells (PBMCs) of 3 healthy donors and 3 prostate cancer patients (23). To explore mRNAs that were associated with prostate cancer, we further extracted 74 mRNAs out of these 2407 mRNAs which were overlapped with these 384 differentially expressed genes (Table S4). A large body of evidence indicated that recurrent prostate cancer following primary therapy is common with high incidence of BCR (2). Therefore, we utilized univariate Cox regression analysis to explore Tregs-specific mRNAs associated with recurrence-free survival (RFS) using gene expression and clinical profiles of 454 patients with prostate cancer from TCGA database. Subsequently, 18 mRNAs were obtained that were significantly associated with RFS of prostate cancer patients (Table S5, P < 0.01). In these 18 mRNAs, 11 mRNAs were adverse indicators, while 7 mRNAs were protective indicators for prostate cancer, which were identified significantly associated with RFS and defined as candidate mRNAs of the prognostic signature for prostate cancer.
Construction of a prognostic signature which can predict RFS of prostate cancer.
To explore a prognostic signature for monitoring RFS of prostate cancer, multivariate Cox regression models were employed using gene expression profiles and clinical information of 454 patients with prostate cancer obtained from TCGA database. Thus, the prognostic signature, named ‘TILTregSig’, was composed of five Tregs-specific mRNAs (SOCS2, EGR1, RRM2, TPP1 and C11orf54). The risk score system was built as follows: risk score= (-0.332 × expression value of SOCS2) + (-0.111 × expression value of EGR1) + (0.286 × expression value of RRM2) + (-0.609 × expression value of TPP1) + (-0.748 × expression value of C11orf54). 454 PRAD patients were dichotomized into high- and low-risk groups according to the cut-off value of 0.946 as the median value of risk score.
Kaplan-Meier (K-M) curves revealed that patients in high-risk groups tended to suffer recurrence (Log-rank P < 0.0001, Fig. 1A). Patients’ recurrence rate was increased in high-risk group compared to that in low-risk group (Fig. 1B). Correlation analysis indicated that RFS status, T, N, Age, cancer status, treatment response and postoperative RX were significantly associated with the risk score (Fig. 1C, P < 0.05). The strip chart showed that the risk of patient mortality and recurrence rate obviously increased and the TNM staging gradually rose as the risk score increased (Fig. 1D). The chi-squared test demonstrated that the high-risk group significantly tended to recurrence (P < 0.0001, Fig. 1E), higher T stage (P = 0.0167, Fig. 1E), higher N stage (P < 0.0001, Fig. 1E), older (P = 0.0068, Fig. 1E), survival with tumor (P < 0.0001, Fig. 1E), no response to treatment (P = 0.0160, Fig. 1E) and no postoperative RX (P < 0.0001, Fig. 1E) compared to the low-risk group of prostate cancer. These results suggested that highly malignant prostate cancer were associated with high-risk score, and our risk score system based on the TILTregSig had tremendous potential to predict RFS for prostate cancer patients.
Cancer stem was proved to be associated with poor prognosis (28). Tathiane M Malta et al. employed OCLR machine-learning algorithm to TCGA datasets to calculate the stemness indices (mRNAsi and mDNAsi) using transcriptomic and epigenetic signatures (24). Thus, we also involved stemness indices derived from study of Tathiane M Malta et al. to validate the clinical association between the TILTregSig and the prognosis of prostate cancer patients. The results showed a significantly positive correlation between the TILTregSig and the stemness indices both in the mRNA expression levels (P = 0.0006, R = 0.160, Figure S2A) and DNA methylation levels (P < 0.0001, R = 0.296, Figure S2B), which was in consistence with our previous results that high-risk score based on the TILTregSig was markedly associated recurrence of prostate cancer.
The TILTregSig has robustly predictive value for prostate cancer patients in ICGC dataset.
To validate the prediction power of the TILTregSig, we further applied the signature to ICGC database. Then, 25 prostate cancer patients were divided into low- and high- risk groups according to the median value of the risk score in ICGC dataset. As shown in Fig. 1F, K-M curves displayed significantly great utility in predicting RFS of prostate cancer. We also found that patients’ recurrence rate was increased in high-risk group compared to that in low-risk group (Fig. 1G). The strip chart of clinical characteristics of prostate cancer patients indicated that the risk of mortality and recurrence rate increased as the risk score increased (Fig. 1H). The chi-squared test showed that patients in high-risk group markedly tended to suffer recurrence (P = 0.0027, Fig. 1I). These results illustrates that high-risk score is associated with recurrence of prostate cancer patients, which is in consistence with our previous results in TCGA datatset.
The TILTregSig is an independently prognostic indicator for prostate cancer patients.
Next, univariate and multivariate Cox regression analysis was carried out to analyze whether the TILTregSig can be an independent predictor for prostate cancer patients. The risk score and other clinicopathological factors were used as covariates. Stage M was not included in this analysis for the reason of only one patient in stage M1 in the TCGA dataset. The results unveiled that the risk score (Multivariate Cox: HR = 1.170, 95% CI = 1.079–1.267; P < 0.001), stage T (Multivariate Cox: HR = 1.749, 95% CI = 1.235–2.474; P = 0.002) and cancer status (Multivariate Cox: HR = 8.084, 95% CI = 4.557–14.33; P < 0.001) were significantly associated with the RFS and could be seemed as independent RFS prognostic factors for prostate cancer patients (Figure S3A, S3B).
To further investigate the clinical potentiality of the risk score in prostate cancer, stratified analysis based on these clinical characteristics was implemented. The results indicated that the TILTregSig seemed more applicable to predict RFS of prostate cancer patients in the subgroups of T1 and T2, N0, younger than 71, White, response to treatment, did not receive postoperative RX and laterality in bilateral (Figure S3C). Meanwhile, patients in high-risk group had significantly poorer clinical outcomes compared to those in low-risk group. These results illustrate that the TILTregSig can be served as an independently predictor for RFS of prostate cancer patients, and still applicable for patients in some subgroups.
The landscape of genetic variations of the TILTregSig in prostate cancer.
Genetic alterations have been found usually confer susceptibilities to prostate cancer (29). GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) is a webtool which can be used to analyse the genetic variation of the genes, and data in the GSCALite were overlapped with the samples derived from TCGA database. Therefore, in order to comprehensively understand the molecular characteristics of the five genes in the TILTregSig, we examined SNV and CNV status of these genes using GSCALite. It was found that the EGR1 (0.6%) exhibited the highest mutation frequency followed by RRM2 (0.2%) and TPP1 (0.2%), while both SOCS2 and C11or54 did not show any mutations in prostate cancer samples (Fig. 2A). In addition, EGR1 had three effective mutation sites, while both RRM2 and TPP1 had one site, respectively (Fig. 2A). Among alteration types, most were focused on the amplification (SOCS2: 7.72%, RRM2: 4.67%, TPP1: 5.69% and C11orf54: 7.93%, Fig. 2B, Table S6), while EGR1 had a widespread frequency of deletion (EGR1: 6.1%, Fig. 2B, Table S6). The investigation of the correlation between CNV and the expression levels of the five genes indicated a significant positive correlation of SOCS2 (Fig. 2C, 2D; Spearman coefficient: R = 0.18, P = 3.2e-04) and TPP1 (Fig. 2C, 2D; Spearman coefficient: R = 0.19, P = 8.8e-05) expression with CNV, which indicated that patients with high expression of SOCS2 and TPP1 were prone to have high CNV load. The above analysis presents a widespread genetic alteration landscape of the five genes in the TILTregSig in prostate cancer patients, suggesting genetic alterations as the molecular mechanism that high-risk score is related to poor prognosis in prostate cancer.
The five mRNAs involved in the TILTregSig have significantly differential methylation between normal and tumor samples in prostate cancer.
Gene methylation plays a vital role in malignant transformation and can be specific to types of cancers including prostate cancer (30). To get a better understanding of the mechanism of the effect of the genes in TILTregSig on tumorigenesis, we analyzed differential methylation of the five genes using GSCALite. Surprisingly, we found that all of the five genes showed significantly differential methylation between normal and prostate cancer samples (P < 0.05, Fig. 2E). Furthermore, SOCS2 (Fig. 2F, 2G; Spearman coefficient: R = -0.61, P < 0.0001), RRM2 (Fig. 2F, 2G; Spearman coefficient: R = -0.39, P < 0.0001), TPP1 (Fig. 2F, 2G; Spearman coefficient: R = -0.21, P = 2.3e-06) and C11orf54 (Fig. 2F, 2G; Spearman coefficient: R = -0.33, P = 4.7e-14) showed significantly negative correlation between gene methylation and expression in prostate cancer, whereas the expression of EGFR (Fig. 2F; P > 0.05) showed no significantly correlation with gene methylation. These findings could contribute to enhancing our understanding of the potential mechanisms on the predictive ability of the TILTregSig in prostate cancer.
The TILTregSig is a stronger predictor for tumor immunity in prostate cancer.
Because the five genes were initially derived from immune cell lines, we consequently investigated whether the signature was related to the tumor immunity. Therefore, we firstly measured the correlation between the TILTregSig and immune-related factors (chemokines, immunoinhibitors, MHCs and receptors), and found that the risk score was commonly associated with these immune-related factors. Especially, the TILTregSig has the highest positive correlation with CCL17 followed by CCL14 among 39 chemokines (Spearman correlation: CCL17.R = 0.206, CCL14.R = 0.202, P < 0.0001, Figure S4). In 24 immunoinhibitors, the TILTregSig presented the highest positive correlation with LGALS9 followed by TGF-β1 (Spearman correlation: LGALS9.R = 0.218, TGF-β1.R = 0.182, P < 0.0001, Figure S4), while showed the highest negative correlation with CD274 followed by TGFBR1 (Spearman correlation: CD274.R=-0,185, TGFBR1.R=-0.175, P < 0.0001, Figure S4). Moreover, the TILTregSig has the highest negative correlation with B2M among 21 MHCs (Spearman correlation: R=-0.296, P < 0.0001, Figure S4), and CXCL1 among 18 receptors (Spearman correlation: R=-0.244, P < 0.0001, Figure S4). Next, further analysis indicated that type Ⅰ IFN response and anti-inflammatory cytokines were significantly enhanced in patients with high-risk sores, while type Ⅱ IFN response was decreased (P < 0.005, Fig. 3A).
In light of these results, we further conjectured that the TILTregSig was correlated with tumor immunity and might have potential to predict tumor immunity. To test this hypothesis, we firstly employed the ESTIMATE algorithm to quantify the immune score, stromal score, ESTIMATE score and tumor purity of prostate cancer patients in TCGA dataset. The results showed that the high-risk group significantly tended to have higher immune scores (P = 0.0052, Fig. 3B, 3C) and ESTIMATE score (P = 0.0036, Fig. 3B, 3C) and lower tumor purity (P = 0.0286, Fig. 3B, 3C) compared to the low-risk group of prostate cancer patients. However, for stromal score, there was no significant difference between low- and high- risk groups (P = 0.061, Fig. 3B, 3C). These results illustrate that the TILTregSig is significantly correlated with tumor immunity, which suggests that the TILTregSig promises to predict tumor immunity in prostate cancer.
TMB (31) and glycolytic activity (11) have been demonstrated to have predictive ability for immune signatures. Next, to compare the predictive ability to tumor immunity, we involved TMB and glycolytic signature in this analysis. ROC curves indicated that the TILTregSig achieved AUC of 0.665 in predicting tumor immune score for prostate cancer, while TME achieved AUC of 0.610 and glycolysis achieved AUC of 0.509 (Fig. 3D). For validating the predictive potential to tumor immunity, we further involved CYT recognized as an immune signature. The TILTregSig achieved AUC of 0.727 in predicting CYT, while TME achieved AUC of 0.626 and glycolysis achieved AUC of 0.459 (Fig. 3E). Our work strongly indicates that the TILTregSig is significantly correlated with tumor immunity and is a stronger predictor than TMB and glycolytic activity for tumor immunity with moderate predictive potential.
The TILTregSig is significantly associated with tumor-infiltrating Tregs in prostate cancer.
To elucidate the mechanism of the correlation between the TILTregSig and tumor immunity, we further investigated the immune functional annotation using the gene set of “c7.all.v7.4.symbles.gmt [immunologic signature]” by GSEA analysis. The results demonstrated that the TILTregSig is highly associated with many immune cells, such as CD8+ T cells, CD4+ T cells, B cells and Tregs et al. (Fig. 4A). Subsequently, to deeply unveiled the relationship between the TILTregSig and tumor infiltrating immune cells, we then evaluated the infiltration levels of immune cells in high- and low- risk groups in prostate cancer samples using marker genes’ expression analysis, CIBERSORT algorithm and ImmuneCell AI database, respectively. We noticed that the consistent results in marker genes’ analysis (Fig. 4B), CIBERSORT algorithm (Fig. 4C) and ImmuneCell AI database (Fig. 4D) were that patients with high-risk score had higher infiltrations of Tregs compared to those with low-risk score (Mann-Whitney U test, P < 0.05, Fig. 4B, 4C, 4D), which suggested a significant correlation between the TILTregSig and Tregs in prostate cancer. For validating these results, we further involved Tregs’ marker genes (FoxP3 and TGF-β1) into our study, and the correlation analysis demonstrated that the TILTregSig was significantly positive related to FoxP3 and TGF-β1 expression (P < 0.05, Fig. 4E), which was consistent with our previous results (Fig. 4B, 4C, 4D). Additionally, RRM2 showed highest correlation with iTreg cells and nTreg cells among these five genes in the TILTregSig (Fig. 4F), suggesting that RRM2 might be a key gene determined the correlation between the TILTregSig and tumor immunity. These results suggest that the TILTregSig may influence tumor immunity mainly by mediating tumor-infiltrating Tregs, and RRM2 may play a vital role in this section, which needs to be verified by further experiments.
The TILTregSig is a powerful predictor for infiltrations of Tregs in prostate cancer.
In view of the relationships between Tregs and the TILTregSig in the TME, we further conjectured that the TILTregSig had potential for predicting infiltrations of Tregs in prostate cancer. To test this hypothesis, prostate cancer patients in TCGA dataset were firstly split into Treg-low group (n = 227) and Treg-High group (n = 227) based on the median infiltration levels of Tregs. Principal component analysis (PCA) demonstrated that the TILTregSig was able to distinguish Treg-low cluster from Treg-High cluster (Fig. 4G), which preliminary hinted us that the TILTregSig had the potential of predicting the infiltrations of Tregs in the prostate cancer.
Additionally, we found that patients in Treg-high group had significantly higher risk scores compared to those in Treg-low group (P < 0.001, Fig. 4H). Particularly, the expression of SOCS2, EGR1, TPP1 and C11orf54 was markedly lower in patients with high infiltration levels of Tregs compared to patients with low infiltration levels of Tregs, while RRM2 was markedly higher (P < 0.0001, Fig. 4I). These results reconfirmed that the predictive potential of the TILTregSig for Tregs’ infiltrations in prostate cancer.
Therefore, to evaluate the predictive power of the TILTregSig for Tregs’ infiltrations, ROC curves were employed in further analysis. Additionally, FoxP3 has been proved as a classic indicator of Tregs (32). Therefore, we also involved FoxP3 in this analysis to compare the predictive ability to the TILTregSig. Surprisingly, we found that the TILTregSig represented high potential as an indicator of Tregs (AUC = 0.897, Fig. 4J), as compared to FoxP3 being a predictor with moderate potential in prostate cancer (AUC = 0.763, Fig. 4J). Moreover, we further constructed a combined model consist of the TILTregSig and FoxP3, and found that the combined model had weaker predictive power compared with the TILTregSig (AUC = 0.767, Fig. 4J). These results suggest that the TILTregSig is a robust and accurate predictor for Tregs in prostate cancer, and its predictive power is stronger than FoxP3.
The TILTregSig has predictive potential as an indicator of response to CIT.
Accumulated evidence demonstrated patients with low infiltrations of Tregs presented a durable clinical response to CIT (33). Our previous data demonstrated that the TILTregSig was associated with Tregs, suggesting that the TILTregSig might be a pivotal factor that mediated the clinical response to CIT. The correlation between our TILTregSig and check-point genes expression indicated that the risk score was markedly correlated with check-point genes expression, and most were positively correlated (P < 0.05, Fig. 5A). Subsequently, for investigated whether the TILTregSig could predict patients’ response to CIT, we utilized five CIT response associated datasets (GSE19423, GSE111636, GSE67501, GSE53922 and Miao D, et al. (22)). We found that the significant therapeutic advantages and clinical response to CIT in patients with high-risk score compared to those with low-risk score were confirmed in bladder urothelial carcinoma (BLCA), kidney renal clear cell carcinoma (KIRC), prostate carcinoma (PRCA) and renal cell carcinoma (RCC). However, in bladder adenocarcinoma (BLAD), patients with high-risk score had lower immunotherapy response rate compared to those with low-risk score (Fig. 5B). Tregs and TME stroma activity, usually mediated immune tolerance of tumors, was also assessed (12, 34). We found that Tregs were significantly activated in tumors with high-risk score in BLAD-GSE19423 dataset, while inhibited in BLCA-GSE111636 and PRCA-GSE53922 datasets (P < 0.05, Fig. 5C). Additionally, TME stroma activity was significantly enhanced in patients with high-risk sore such as the activation of epithelial-mesenchymal transition (EMT) in BLAD-GSE19423 dataset, while was decreased in tumors with high-risk score in BLCA-GSE111636 and PRCA-GSE53922 datasets (P < 0.05, Fig. 5C). The above data imply that the TILTregSig is correlated with CIT response and might have potential for predicting CIT response. To assess the predictive power of the TILTregSig for CIT response, ROC curves were employed in this work. The results showed that the TILTregSig achieved AUC of 0.710, 0.677, 0.799, 0.614 and 0.695 for BLAD-GSE19423, BLCA-GSE111636, KIRC-GSE67501, PRCA-GSE53922 and RCC respectively in predicting the response to CIT, which implied that the TILTregSig was a potential and robust biomarker for response assessment of CIT with moderate predictive potential (Fig. 5D). In summary, our work strongly indicated that the TILTregSig was significantly correlated with immunotherapy response, and the established TILTregSig would contribute to predicting the response to CIT.
The TILTregSig is a promising marker of therapeutic resistance in pan-cancers.
Considering Treg-cells have been indicated that usually promotes resistance to cancer therapy (35), we investigated whether our TILTregSig was correlated with cancer therapeutic resistance. To unveil the relationship between the TILTregSig and therapeutic resistance, we utilized GSEA analysis. As shown in Fig. 6A, GSEA predicted that the TILTregSig was significantly associated with resistance to different therapies, including salirasib, endocrine therapy, doxorubicin and SB216763 in prostate cancer (P < 0.05, Fig. 6A). Next, a landscape plot was generated by GSCALite to depict the relationships between the five genes expression in the TILTregSig and drug responses in pan-cancers. The bubble heat map showed that some genes exhibited significant correlations with lower half inhibitory centration (IC50) data. In detail, EGR1, TPP1 and SOCS2 conferred drug resistance, while RRM2 and C11orf54 exhibited drug sensitivity in pan-cancers (Fig. 6B). These results imply the TILTregSig as a promising indicator for therapeutic resistance in pan-cancers.