For our independent cohort (NanTong cohort), a total of 94 prostate cancer patients were recruited. All of patients underwent prostate biopsy for diagnosis of prostate cancer. The follow-up protocol involved conducting telephone calls subsequent to the initial diagnosis. Prostate cancer tissues were obtained during tumorectomy procedures and were immediately frozen at -80°C for subsequent analyses. The extraction of tissue RNA was performed in accordance with the manufacturer’s instructions. Clinicopathological findings were assessed based on the tumor–node–metastasis (TNM) classification system. Informed consent was obtained from all patients.

In this study, two publicly available databases were utilized: The Cancer Genome Atlas Prostate Adenocarcinoma (TCGA PRAD) (https://portal.gdc.cancer.gov/) and the Memorial Sloan Kettering Cancer Center (MSKCC) (http://cbio.mskcc.org/cancergenomics/prostate/data/). The TCGA PRAD and MSKCC cohorts were acquired for subsequent analyses. Within the TCGA database, data included transcripts per million (TPM) of RNA sequencing and matched somatic mutation datasets of PCa, were obtained using the TCGA biolinks package in the R software. Additionally, matched tumor purity estimated by immunohistochemistry (IHC) data was downloaded (15), and matched batch information was obtained from https://bioinformatics.mdanderson.org/MQA/. For cases where a gene symbol had multiple expression measurements, the measurement with higher expression was retained. Due to significant RNA degradation in a portion of TCGA PRAD samples, 333 cases were utilized. Ultimately, 282 samples with survival information related to biochemical recurrence were included in our analysis. For the MSKCC cohort sequenced by microarray, normalized log2 mRNA expression data was downloaded.

Referring to the previous study (16), we selected principal component analysis to analyzed and visualize batch effects of TCGA PRAD.

The CIBERSORT algorithm, a computational method used to estimate the composition of different immune cell types in a tissue sample from bulk gene expression profiles, was applied to the TCGA PRAD and MSKCC dataset (17). Analyses were performed using 1,000 permutations and default statistical parameters in reference to LM22 matrix. The threshold for categorizing CD8+ T cell infiltration into high and low groups was established by comparing the differences in biochemical recurrence. According to different sequencing methods in this study, the cutoff value was calculated for each data, respectively. The cutoff value selected was the one that yielded the lowest p-value. The cytolytic score, which serves as an indicator of local immune cytolytic activity, was calculated as the geometric mean of gene expression values for granzyme A and perforin (18). To estimate T-cell exhaustion level, the murine T-cell exhaustion signature was obtained (19, 20). The murine genes were manually converted to their corresponding human gene equivalents. The degree of T-cell exhaustion was assessed by calculating the mean expression of up-regulated genes minus the mean expression of down-regulated genes.

The somatic mutations were analyzed using the “maftools” R package. The tumor mutation burden (TMB) score was calculated for each patient in accordance with established methods. Oncoplots were utilized to visually present the somatic mutation signatures. The identification of cancer driver genes was achieved through the implementation of the oncodrive CLUST algorithm. Differentially mutated genes were identified by Fisher’s exact test.

Before differential expressed gene analysis, the genes with TPM equal to 0 among more than 10% samples were filtered in TCGA PRAD. Differential gene analysis was performed using multiple statistical approaches, including the Wilcoxon rank sum test and signed rank test, DESeq2 and edgeR (21). The genes with an adjusted p-value less than 0.05 and an absolute log2 scaled fold change more than 0.5849625 were further analyzed. The false discovery rate (FDR) method was employed for adjusting p-values.

The Least Absolute Shrinkage and Selection Operator (LASSO) method was employed to identify stable prognostic candidate genes using biochemical recurrence as the endpoint. The prognostic candidate genes were subsequently validated using the log-rank test, as well as univariate and multivariate Cox proportional hazards regression analysis, with biochemical recurrence as the endpoint. To compute the risk score for each patient, the following formula was utilized:

Risk scores were computed by Cox regression coefficients of the adjusted covariates in both the TCGA and MSKCC cohorts. Additionally, the predictive value of the risk score was evaluated using the area under the receiver operating characteristic curve.

Gene annotation enrichment analyses were conducted on the DEGs between the low and high CD8+ T cell groups using the R package clusterProfiler (22). The analysis included identification of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) terms. Statistical significance was determined using an adjusted p-value cutoff of< 0.05. Additionally, gene set enrichment analysis (GSEA) was performed to identify consistent biological differences between the high and low CD8+ T cell groups, with an adjusted p-value cutoff of< 0.05 indicating statistical significance.

The MLXIPL expression in NanTong cohort was measured by qRT-PCR. The relative expression of MLXIPL mRNA was calculated by 2^−ΔCt method with the normalization to GAPDH. Primer sequences of MLXIPL were F: AAGATCCGCCTGAACAACG and R: CACTTGTGGTATTCCCGCATC. Primer sequences of GAPDH were F: CTGGGCTACACTGAGCACC and R: AAGTGGTCGTTGAGGGCAATG.

The immunohistochemistry (IHC) analysis was conducted using rabbit anti-ChREBP (1:200, ab101500, Abcam, Cambridgeshire, England) following the manufacturer’s instructions. Immunostaining intensity was categorized into four grades: 0 (no expression), 1 (mildly positive), 2 (moderately positive), and 3 (markedly positive). The proportion of positive-staining cells was assessed and categorized into five grades: 0 (0%), 1 (1–25%), 2 (26–50%), 3 (51–75%), and 4 (>75%). To generate the IHC score, the percentage of tumor cells showing positivity and the staining intensities were multiplied.

For comparisons between two subtypes, the Wilcoxon rank sum and signed rank tests were employed. Discrete data comparisons were conducted using Fisher’s exact test. Spearman’s correlation analysis was utilized to explore the relationships. All statistical tests were two-sided, and p< 0.05 was considered statistically significant unless otherwise stated. The thresholds for p-values were set at 0.05, 0.01, and 0.001 (*, ** and ***, respectively). All statistical analyses were performed using R software, version 4.3.1.

We investigated the batch effect of TCGA PRAD dataset ( Supplementary Figure S1 ). Principal component analysis indicated no batch effects. The baseline characteristics of these two cohorts are present in Supplementary Table S1 . In line with previous studies (11, 12), patients with higher levels of infiltrated CD8+ T cells exhibited poorer prognosis trend ( Figure 1A ). MSKCC cohort supported this result ( Figure 1B ).

The presence of higher densities of CD8+ T cells has been suggested to be associated with more advanced tumors (10). We examined clinicopathological characteristics, namely age, tumor purity, histopathological subtype, T and N stages, and Gleason score, in relation to the level of CD8+ T cells in TCGA PRAD ( Supplementary Figures S2A–F ) and clinicopathological characteristics, namely age, prostate specific antigen (PSA) level, T and N stages, and Gleason score and ERG-fusion status, in relation to the level of CD8+ T cells in MSKCC ( Supplementary Figures S2I–O ). The results indicated that compared with other histopathological type, less CD8+ T cells infiltrated in acinar PCa and age, tumor purity, PSA level, T, N stages, Gleason grade and ERG-fusion are not correlated to CD8+ T cell infiltration level.

Given the negative association of exhausted T cells with ccRCC prognosis (13, 14), cytolytic score and T cell exhaustion levels were calculated. Greater immune cytolytic activity was evident in the high CD8+ T cell group ( Supplementary Figures S2G, P ). Comparable levels of exhausted T cells were observed in the TCGA PRAD cohort. In the MSKCC cohort, levels of exhausted T cells were elevated in the high CD8+ T cell group ( Supplementary Figures S2H, Q ). This suggests that T cell exhaustion may partially elucidate the poor prognosis associated with CD8+ T cells.

In summary, our study reaffirmed the observation that CD8+ T cells are linked to an unfavorable prognosis based on available data. Our findings suggest that CD8+ T cells are associated with a poorer prognosis, rather than a higher degree of malignancy resulting in increased infiltration of these T cells.

To gain further insights into how CD8+ T cells contribute to the poor prognosis of PCa, we explored whether this involvement is mediated through the TME. Within the TME, immunosuppression can be generated through two possible mechanisms: (1) the presence of differential mutations, which may modulate the immune response in distinct ways; and (2) the presence of tumor-infiltrating T cells that can be suppressed through feedback-induced expression of checkpoint molecules and recruitment of immunoregulatory cells (18).

To investigate the potential role of mutations in mediating the association between CD8+ T cells and poor prognosis in PCa, we conducted somatic mutation frequency analysis. Both the low and high CD8+ T cell groups exhibited low TMB status and similar mutational patterns ( Figure 2A ). High rate of SPOP, TP53, TTN and FOXA1 mutations were found. Missense mutations were found to be the most prevalent variant classification. Notably, SPOP were identified as driver genes despite CD8+ T cell infiltration level ( Supplementary Figure S3A ). SNPs emerged as the most frequent variant type ( Supplementary Figure S3B ). The results showed comparable frequencies of transitions and transversion between the low and high CD8+ T cell groups ( Figure 2B ). Results of Fisher’s exact test confirmed that CD8+ T cell infiltration level was not implicated in gene mutations ( Supplementary Table S2 ). Collectively, these results indicate that the association between CD8+ T cells and poor prognosis in PCa is not driven by somatic mutations.

Subsequently, we explored the correlation between immune checkpoint genes and the level of CD8+ T cell infiltration. After filtering the genes utilized in the LM22 matrix, six genes—specifically, ITGAL, CD74, HAVCR2, CD274, SIGLEC15, and TIGIT—were selected based on a review of the literature (23, 24). Increased ITGAL, CD74 and TIGIT were observed in the high-density CD8+ T cell subgroup ( Figure 2C ). The relative abundances of 22 infiltrating immune cells were investigated. Notably, elevated immunosuppressive regulatory T cells (Tregs) were observed ( Figure 2D ). These findings suggested that CD8+ T cells may facilitate an unfavorable prognosis through a feedback mechanism involving the abnormality of immune checkpoint gene expression and recruitment or differentiation of immune cells specialized in immune suppression.

We hypothesized that protein-coding genes influenced by CD8+ T cells should (1) alter in response to CD8+ T cell infiltration and (2) play a significant role in PCa prognosis. Initially, a differential expression analysis was conducted to compare the low and high CD8+ T cell groups. Three distinct methods, Wilcoxon rank sum test, DESeq2, and edgeR, were utilized. Genes with an adjusted p-value below 0.05 and an absolute log2-scaled fold change exceeding 0.5849625 were identified as DEGs. In total, 74, 233 and 356 DEGs were identified using each respective method. In total, 55 genes exhibited an increase, while 2 genes showed a decrease in the high CD8+ T cell group ( Figures 3A, B ). Among them, 34 genes were involved in the LM22 matrix. Furthermore, we conducted a screening of genes associated with the prognosis of PCa. Among 55 DEGs, MLXIPL was identified using LASSO Cox analysis with the TCGA dataset ( Figure 3C ). Specifically, MLXIPL showed elevated expression in the high CD8+ T cell group (p = 2.06E-4, Supplementary Figure S4A ) and exhibited a significant correlation with the level of CD8+ T cells (rho = 0.25, p = 2.92E-5, Supplementary Figure S4B ). These findings were confirmed using the MSKCC dataset (p = 3.81e-10, Supplementary Figure S4C ; rho = 0.67, p = 9.54E-20, Supplementary Figure S4D ).

We conducted a functional enrichment analysis to explore the biological pathways that exhibit alterations in response to CD8+ T cell infiltration. Based on KEGG database, we identified a total of 23 significantly enriched pathways majorly related to such as cell communication, immune responses, and immune checkpoint pathways ( Figure 3D ). The most of significant altered pathways were linked to functions cell communication, antigen presentation, and immune responses. Based on KEGG, GSEA revealed 31 upregulated and 11 downregulated pathways. Consistent with the results of the enrichment analysis, these altered pathways were associated with upregulated cell communication, antigen presentation, and immune responses ( Supplementary Figure S5 ). Collectively, these findings suggested that changes in cell communication occur in response to CD8+ T cell infiltration, which may contribute to an unfavorable prognosis of PCa.

Elevated expression of MLXIPL mediated by CD8+ T cells was correlated with an increased rate of biochemical recurrence (log-rank p = 2.30E-02, Figure 4A ). Univariate and multivariate Cox proportional hazards regression models validated that MLXIPL significantly contributed to a poor prognosis (crude p = 4.20E-03, adjusted p = 6.70E-02, Supplementary Figure S6 ). Moreover, the risk score, computed using MLXIPL, age, T stage and Gleason score, exhibited strong predictive capabilities for PCa prognosis (1-year AUC = 0.70; 3-year AUC = 0.72; 5-year AUC = 0.82, Figure 4B ). The C-index value of the model was 0.72 (95% CI: 0.64-0.80).

Consistently, MLXIPL remained an independent predictor, unaffected by variables such as age, PSA levels, tumor purity, T and N stages, and Gleason score. In the high MLXIPL group, we observed a decreased acinar adenocarcinoma ratio and an increased number of nodes ( Supplementary Figures S7A–F ). Additionally, reduced mutation frequency and TMB were observed in high MLXIPL group ( Supplementary Figure S7G ). All TMBs were less than 10 mutations per megabase (MB), indicating a low TMB level. Moreover, somatic mutation frequency analysis revealed no mutation was associated with MLXIPL expression ( Supplementary Table S3 ). Consistent with the results above, ITGAL, CD74, and TIGIT showed elevated expression in response to MLXIPL ( Figure 4C ). Furthermore, MLXIPL was associated with the infiltration fraction of several immune cells ( Figure 4D ). Specifically, the high MLXIPL group demonstrated an increased infiltration of regulatory T cells.

To validate the findings, we investigated the role of MSKCC in PCa prognosis using the MSKCC cohort. High MLXIPL was associated with poor prognosis (log-rank p = 1.6E-02, Figure 5A ). Univariate and multivariate Cox proportional hazards regression models confirmed that MLXIPL was associated with poor prognosis (crude p = 1.46E-02, adjusted p = 4.23E-01, Supplementary Figure S8 ). The AUC values of the risk score, calculated based on MLXIPL, T stage, and Gleason score, for predicting 1-, 3-, and 5-year biochemical recurrence were 0.94, 0.90, and 0.81, respectively ( Figure 5B ). The C-index value of the model was 0.85 (95% CI: 0.80-0.91). MLXIPL was implicated in N stage, while not correlated to age, PSA level, T stage, Gleason score, metastasis and ERG-fusion status ( Supplementary Figure S9 ). Four immune checkpoint genes, namely ITGAL, HAVCR2, SIGLEC15 and TIGIT, increased and CD74 decreased ( Figure 5C ). Among them, ITGAL and TIGIT also elevated in TCGA. Additionally, MLXIPL were related to infiltration fraction of several immune cells ( Figure 5D ). As anticipated, Tregs infiltration increased in the high MLXIPL group.

In summary, these results suggest that CD8+ T cells may modulate MLXIPL, thereby affecting PCa prognosis by upregulating ITGAL and TIGIT and recruiting immunosuppressive Tregs.

Finally, we examined the role of MLXIPL in our own cohort, comprising 94 PCa patients with follow-up information ( Supplementary Table S4 ). Briefly, the mean ages were 60.02 ± 7.08 and 61.04 ± 7.37 in the low and high MLXIPL groups, respectively. The results showed that MLXIPL expression was not correlated to age, T stage, N stage and Gleason score ( Supplementary Figure S10 ).

Consistent with previous findings, MLXIPL promoted to poor prognosis (log-rank p = 7.2E-04, Figure 6A ; crude HR = 2.22, 95% CI: 1.33-3.72, p = 2.42E-03; adjusted HR = 2.57, 95% CI: 1.42-4.65, p = 1.76E-03, Supplementary Figure S11 ). Moreover, when combined with clinicopathologic characteristics, MLXIPL demonstrated high predictive performance. The AUC values of the risk score, calculated based on MLXIPL, age, T stage, and Gleason score, for predicting 1-, 2-, and 3-year overall survival were 0.77, 0.75, and 0.80, respectively ( Figure 6B ). The C-index value of the model was 0.76 (95% CI: 0.65-0.86). To further improve prognostic prediction, a nomogram model was established using multivariable Cox regression in the NanTong cohort to estimate the 1-, 2-, and 3-year biochemical recurrence, incorporating MLXIPL expression, age, T stage, and Gleason grade as variables ( Figure 6C ).

Furthermore, we investigated the protein levels of MLXIPL in human prostate tissue microarray to understand its role in PCa tumorigenesis. The baseline characteristics of 59 PCa patients, including 55 pairs of normal adjacent tissue and tumor samples, are presented in Supplementary Table S5 . In summary, the mean age of the 59 men was 66.83 ± 5.50. IHC results revealed an elevation in MLXIPL protein levels in the tumor tissue (unpaired p = 1.03E-2, paired p = 3.27E-2, Supplementary Figure S12 ).