Different metabolic paradigms and distribution of regulatory T cells between primary and lymph node metastasis prostate cancer

ScRNA-seq and data analyzation

The research was approved by the Ethics Committee of Shanghai East Hospital. PCa tissue and LN with cancer metastasis were acquired from 4 patients who were diagnosed with PCa and underwent radical prostatectomy. These samples were then sent for cell isolation and scRNA-seq using ×10 chromium single-cell platform (×10 genomics). Further data standardization and dimensionality reduction were completed by R software (version 4.1.3, R Foundation for Statistical Computing, Vienna, Austria).

Characterization of regulatory T cell and malignant cell in ScRNA-seq

R package Seurat^[31] (version 4.0, Satija Lab, New York, USA) was utilized to analyze the scRNA-seq data. Doublets were removed using the DoubletFinder_v3 algorithm. Data preprocessing included the following quality control filters: (1) genes detected in ≥3 cells, (2) cells containing ≥250 expressed genes, (3) cells with ribosomal RNA counts between 100 and 5,000, and (4) cells exhibiting mitochondrial gene content <25%. The filtered dataset was log-normalized using R.^[32] For dimensional reduction, we identified 2,000 highly variable genes based on mean expression-dispersion relationships. Principal component analysis (PCA) was performed for feature selection with the top 30 principal components (PCs) subsequently used for t-distributed stochastic neighbor embedding (t-SNE) analysis (resolution = 0.4) to perform unsupervised cell clustering.

Immune cell clusters were further annotated using the Human Primary Cell Atlas through the SingleR package (version 1.8.1, Aaron Lun, Ontario Institute for Cancer Research, Canada) as well as identification of canonical immune cell markers from established literature. Especially, clustering of T cells has proceeded to characterize the subtypes of Treg. Moreover, inferCNV algorithm (https://github.com/broadinstitute/inferCNV) was performed to locate malignant cells. For visualization, the top 5 marker genes of each cluster were displayed in dot plots, while cellular composition was represented using bar plots.

Resolving the metabolic and transcriptional features of 3 subtypes of Treg

To explore the functional difference between 3 Treg subtypes, we applied R package scMetabolism (https://github.com/wu-yc/scMetabolism) to quantify the activities of 33 typical metabolic pathways. 3 algorithms were integrated with this software to calculate the enrich score of individual pathways, including VISION,^[33] AUCell,^[34] and ssGSEA. Moreover, Monocle3 (version 0.2.0, Trapnell Lab, Seattle, USA) was employed to infer the time evolution of the subtypes of Tregs (https://cole-trapnell-lab.github.io/monocle3). This software executed trajectory inference and pseudotime analysis to predict the evolution of each cell type by comparing the difference in genomic profiling. Following SCENIC^[34] (version 0.1.5, Aerts Lab, VIB-KU Leuven, Belgium) analysis further investigated the transcriptional features of 3 Treg subtypes. Master transcription factor (TF) associated with the 3 cell states was identified based on the co-expression modules of this package. We obtained motif databases from https://resources.aertslab.org/cistarget/databases/. Cell–cell communication analysis was conducted using CellphoneDB (version 2.1.7), with edge weights normalized for cross-network comparison.^[35] Results were visualized through circular plots depicting the interactions of cell populations, and a bubble plot was illustrated to count all important ligand-receptor pairs during intercellular signaling.

The correlation between metabolic pathways and Treg cells

To further estimate the correlation between metabolic and Treg cells, the differential analysis between primary and metastatic LN samples from scRNA-seq, N0 and N1 samples from The Cancer Genome Atlas (TCGA)^[36] cohort was conducted. Some pathways with common differential trends were selected for further clinical correlation analysis. Using Gene Set Variation Analysis (GSVA),^[37] we quantified the activity levels of these metabolic pathways in each tumor sample. Finally, we generated correlation scatter plots to evaluate the relationship between metabolic pathway activity and Treg subtype abundance.

Cluster analysis based on inhibitory regulatory T (iTreg) and effector regulatory T (eTreg) cells

Given the significant association between iTreg, eTreg, and clinical characteristics in PCa patients, we performed cluster analysis on these Treg subsets. After removing duplicate genes, we identified 225 marker genes for univariate Cox regression analysis. A heatmap was generated to visualize the expression patterns of these 225 genes in N0 versus N1 samples. Using the R package ConsensusClusterPlus,^[38] we conducted unsupervised hierarchical clustering on 499 tumor samples based on the messenger RNA expression levels of prognostic genes. Then, the survival analysis was conducted, and Kaplan-Meier (KM) curve illustrated prognostic difference between 2 clusters. To further investigate the relationship between these clusters, we generated a heatmap to display differential gene expression and performed correlation analysis to assess clinical associations.

Functional difference and immune infiltration analysis

To elucidate the functional difference between 2 sub-clusters, Gene Ontology (GO)^[39] and Kyoto Encyclopedia of Genes and Genomes (KEGG)^[40] enrichment analysis were employed to depict the functional difference between 2 clusters with “clusterProfiler”^[41] package.

Furthermore, the ESTIMATE^[42] algorithm was utilized to calculate the StromalScore, ImmuneScore, ESTIMATEScore, and TumorPurity of each sample, to discover the TIME difference between 2 clusters. These scores represented the stromal cell infiltration, immune cell infiltration, and tumor purity. Furthermore, ImmuCellAI database^[43] was employed to explore the difference in TIME between 2 clusters, which obtained accurate data on the abundance of 24 immune cell types, including 18 types of T cells and 6 other immune cells (B cells, natural killer cells, monocytes, macrophages, neutrophils, and dendritic cells).^[43]

Clinical sample and data collection

Clinical samples (44 pairs LN metastases and 15 pairs primary PCa samples) were collected in the form of paraffin section from patients undergoing prostate surgery at Shanghai East Hospital from 2019 to 2023. PCa samples were obtained by RP, including prostate and LN tissues. Each sample was histologically evaluated on hematoxylin/eosin-stained tissue sections by two experienced uropathologists based on the National guidelines for diagnosis and treatment of prostate cancer 2022 in China.^[44] Each tumor sample must contain at least 70% tumor tissue. The clinical data were obtained from their hospitalization information and follow-up, including age, Gleason score, Prostate-specific antigen (PSA), tumor stage, biochemical recurrence (BCR), etc.

The validation of Treg cells infiltration proportion

IHC staining analysis of Foxp3 and CD45 in clinical samples was used to assess Treg infiltration and distribution. Foxp3 encodes a TF, which is a key regulatory gene in the development of Treg cells and plays a leading role in the development, differentiation, and immune suppressor function.^[18] Foxp3 in naturally occurring Treg specificity expression can be used as a specific marker to distinguish Treg cells from other T cells.^[19] In addition, CD45 (leukocyte common antigen, LCA), which was expressed on all leukocytes, was utilized to exclude the effect of the minimal expression of Foxp3 in other cells.^[45]

Foxp3 and CD45 antibody was acquired from Abcam (ab20034, ab40763). The secondary antibody used was horseradish peroxidase goat anti-rabbit immunoglobulin G H&L (ab97051). Immunostaining intensity was assessed by two independent pathologists blinded to clinical data, following the World Health Organization histological classification criteria for PCa.^[46] Staining intensity was scored as 0 (negative), 1 (weak), 2 (medium), 3 (strong), or 4 (very strong). The positive rate of target cells was evaluated at least 5 high-power fields, Foxp3: 1 NC (The average number of cells per high-power field) <5 or PC (the percentage of positive cells) < 1%), 2 (5 ≤ NC < 10 or 1% ≤ PC < 5%), 3 (10 ≤ NC < 20 or 5% ≤ PC < 10%), 4 (20 ≤ NC < 50 or 10% ≤ PC < 25%), 5 (NC ≥ 50 or PC ≥ 25%); CD45: 1 (NC < 10), 2 (10 ≤ NC < 50), 3 (50 ≤ NC < 100), 4 (100 ≤ NC < 200), 5 (NC ≥ 200).

Statistic and software

Data processing and bioinformatic analyses were accomplished by R software (version 4.1.3), Statistical Package for the Social Sciences (version 26.0, IBM Corp., Armonk, NY, USA) and GraphPad Prism (version 9.0.0, GraphPad Software, San Diego, CA, USA). Packages limma (version 3.50.0, Walter and Eliza Hall Institute, Melbourne, Australia),^[47] timeROC (version 0.4, Paul Blanche, University of Copenhagen, Denmark),^[48] survminer (version 0.4.9, Alboukadel Kassambara, Montpellier, France),^[49] Seurat,^[31] GSVA (version 1.42.0, Robert Gentleman, Fred Hutchinson Cancer Center, USA), glmnet (version 4.1-3, Trevor Hastie & Jerome Friedman, Stanford University, USA), etc. were employed for analyses with proper citation. Wilcox or Kruskal–Wallis test was applied for comparisons between two and more groups. Pearson and Spearman’s rank correlation were adopted to estimate the statistical correlation of parametric or non-parametric variables. Log-rank test was utilized for survival analysis. Two-sided P < 0.05 was considered significant threshold for all statistical tests.

T cell subtypes and their infiltration in primary and LN tissue of PCa

A total of 23097 genes and 33809 cells were detected in the scRNA-seq data of PCa. 9 cell types and 28 cell clusters were identified in the t-SNE analysis, distributing evenly among 4 PCa samples [Figure 1a-c]. Of them, CD8 T cell took the biggest part, followed by CD4 T cell and B cell, homogeneously infiltrated in the 4 samples [Figure 1d]. Cancer cell was located at the bottom left of the t-SNE plot, having higher stemness than other cell clusters [Figure 1e and f]. Further analysis characterized the marker genes for each cell type. Malignant cells were featured by the expression of EPCAM and KRT8 [Figure 1g]. The infiltration of abundant immune cells suggested a potential exclusive TIME in PCa.

Close

Figure 1:

The identification of cell type based on scRNA-seq. (a) 23097 genes and 33809 cells in the PCa scRNA-seq; (b and c) 9 cell types and 28 cell clusters identified with t-SNE; (d) The distribution of each cell cluster in 4 samples; (e and f) Malignant cells and their stemness analysis in PCA samples; (g) Marker genes for each cell type (PTPRC: Leukocytes, CD3D: T cells, CD4: T helper cells, regulatory T cells, mononuclear cells/macrophages, CD8A and CD8B: cytotoxic T cell, CD79A: B cells, NKG7: Natural killer cells and activated T cells, LYZ: Myeloid cells, PLVAP: Endothelial cells, ACTA2: Smooth muscle cells and myofibroblasts, EPCAM: Epithelial cell, KRT8: Monolayer epithelial cell). (pca01, 03: Primary PCa tissue, pca02, 04: LN metastasis tissue of paired PCa). scRNA-seq: Single-cell RNA sequencing, PCa: Prostate cancer, PCA: Principal component analysis, t-SNE: t-distributed stochastic neighbor embedding.

Export to PPT

There were 15559 T cells and they were then divided into 8 subtypes, including Naive_CD8+_T, Naive_CD4+_T, Treg, CD8+_Tef, CD8+_T_EMRA, and so on [Figure 2a and b]. Of them, CD8_Tef occupied almost half part of the whole tissue in primary PCa samples but was with much less proportion in the metastasized LN. On the contrary, Treg had a rather higher percentage in the positive LN samples than primary PCa samples [Figure 2c]. This difference implied a potential inhibitory TIME in the metastasized LN samples. Not surprisingly, the top 5 marker genes of Treg contained RTKN2, LAIR2, CTLA4, TNFRSF18, and TNFRSF4, well-known immune checkpoints to block immune response [Figure 2d].

Close

Figure 2:

T cell subtypes in PCa samples. (a and b) Dimensionality reduction clustering and annotation of T cells; (c) The distribution of T cell subtypes in PCa samples; (d) Expression of marker genes for each T cell subtype. UMPA: Uniform Manifold Approximation and Projection, PCa: Prostate cancer.

Export to PPT

Treg subtypes and their metabolic feature in primary and LN tissue of PCa

3 clusters of Treg were characterized in the dimensionality reduction analysis [Figure 3a and b], including iTreg, eTreg, and DP-Treg. eTreg was featured by the high expression of NKG7, GZMA, and GZMK, related to cytotoxic function. iTreg was seen with high expression IL7R and CCR7, possibly inducing an inhibitory immune response. DP-Treg, however, had the distinction of both eTreg and iTreg [Figure 3c].

Close

Figure 3:

Treg cell subtypes in PCa samples. (a and b) Dimensionality reduction clustering and annotation of Treg cells; (c) Expression of marker genes for iTreg cells, eTreg cells, and DP-Treg cells. PCa: Prostate cancer, iTreg: Inhibitory regulatory T, eTreg: Effector regulatory T, DP-Treg: double positive regulatory T.

Export to PPT

As expected, LN tissue of PCa showed increased infiltrating proportion of Treg and Naive_CD4+_T but less percentage of CD8+_Tef, suggesting a probable inhibitory TIME in metastasized tissue. Correspondingly, LN tissue was seen with evidently high infiltration of iTreg, but less existence of eTreg [Figure 4a and b]. Further analysis found that they were totally different in metabolic function. While iTreg had higher pathway activities of Pyruvate metabolism, Oxidative phosphorylation, and Pentose phosphate, eTreg was observed to be more active in sulfur metabolism, Primary bile acid biosynthesis, Nitrogen metabolism, and Linoleic acid metabolism [Figure 4c]. The metabolic profile of these different Treg cell types may be related to their immune function and, therefore, we may regulate the subtype of Treg and create a positive TIME for tumor by changing the key metabolic enzymes of these pathways.

Close

Figure 4:

Different distribution and metabolic patterns of Treg subsets. (a) The different distribution of T cell subtypes between primary and metastasis LN PCa tissues; (b) The different distribution of Treg cell subtypes between primary and metastasis LN PCa tissues; (c) Metabolic differences among different subsets of Treg. Treg: Regulatory T, LN: Lymph node, PCa: Prostate cancer.

Export to PPT

Evolution of 3 subtypes of Treg and cell-cell communication

To investigate the developmental trajectory of Treg subtypes, we performed Monocle3 pseudotemporal analysis. Consistent with previous findings, effector Tregs (eTregs) showed significantly elevated expression of CCL4 and GZMA, while iTregs were characterized by EEF1A1 expression. DP-Treg, however, had the features of both of them [Figure 5a]. Subsequent pseudo-time analysis found DP-Treg located at the sequential starting point. This indicated that DP-Treg may be the precursor of eTreg and iTreg, evolving into these two completely distinct descendants [Figure 5b and c]. Moreover, we found that Treg cell had close cell–cell communications with other cell types. The most frequently communicated CD8+_Tex_Proliferating cell possessed 101 interaction pairs with Treg, followed by NKT_Cytotoxic cell and Naive_CD4T cell with an interaction number of 76 and 60, respectively [Figure 5d and e]. These findings suggest that Treg-mediated immunosuppression may primarily occur through interactions with proliferating exhausted CD8+ T cells, highlighting the need to further elucidate these regulatory mechanisms to modulate the immunosuppressive tumor microenvironment.

Close

Figure 5:

Association between Treg subtypes. (a) Molecular evolution of 3 Treg subtypes; (b and c) Pseudo-time analysis of Treg subtypes; (d and e) Cell–cell communications between Treg and other cells. Treg: Regulatory T.

Export to PPT

The metabolic feature and clinical correlation

To identify metabolic pathways associated with LN metastasis in PCa, we performed differential expression analysis comparing: (1) primary versus metastatic LN samples from scRNA-seq data, and (2) N0 versus N1 samples from the TCGA cohort. Our analysis revealed significant suppression of three metabolic pathways in iTreg cells and LN metastatic samples: Alpha-Linolenic acid metabolism, Linolenic acid metabolism, and Arachidonic acid metabolism [Figure S1a and b]. Based on the enrich score of each pathway and BCR time, KM curves of 3 pathways illustrated that these pathways negatively correlated with the prognosis of PCa patients [Figure S1c], (P < 0.001). The activities of 3 pathways were related to the PCa pathological stage, which were high activated in advanced patients with LN metastasis or T3/4 stage [Figure S1d]. In addition, different correlation of metabolism pathways with iTreg or eTreg cells is illustrated in Figure S2.

Two sub-clusters based on Treg cells

We identified 225 marker genes differentially expressed in iTreg and eTreg cells (DP-Treg cells had no significant difference in LN metastasis samples) for univariate Cox regression analysis. 43 genes were identified as prognostic genes with P < 0.05 [Figure S3a]. The expression difference between N0 and N1 patients of these genes was demonstrated in a heatmap [Figure S3b]. The samples from TCGA were clearly divided into 2 clusters based on 43 prognostic genes [Figure S3c and d]. The prognosis of patients in cluster 1 was significantly poorer than those in cluster 2 with P = 0.015 [Figure S3e]. The expression of 43 genes in cluster 1 and 2 was illustrated in a heatmap. Most genes showed clear differences in expression between the two groups, such as CD4, which was highly expressed in iTreg, was highly expressed in cluster1 [Figure S4a]. The difference in clinical features suggested that patients in cluster 1 may have higher Gleason score, T stage, and N stage [Figure S4b-e]. These all suggested that cluster 1 represents a more aggressive tumor phenotype with molecular characteristics resembling iTreg cells, corresponding to poorer clinical outcomes.

The function of enrichment and immune microenvironment analysis

The marker genes that were used to classify patients into 2 clusters were also utilized for KEGG and GO enrichment analysis. The results showed that these genes were significantly involved in immune response-activating cell surface receptor signaling pathway, immune response-activating signal transduction, lymphocyte-mediated immunity, Th17 cell differentiation, Th1 and Th2 cell differentiation, etc, which were all closely related with tumor immune infiltration and microenvironment [Figure S5a and b].

The StromalScore, ImmuneScore, ESTIMATEScore, and TumorPurity of each sample were calculated with ESTIMATE algorithm, and the difference between 2 groups was illustrated under the violin diagram. The result showed that cluster 1 had higher StromalScore, ImmuneScore, and ESTIMATEScore, while cluster 2 had higher TumorPurity [Figure S5c], which suggested that patients in cluster 1 had a higher infiltration of immune cells.

To further investigate tumor-infiltrating immune cell subsets of sub-clusters, ImmuCellAI was utilized to estimate the difference in the infiltration of 24 subsets in TCGA. The results of differential analysis of immune cells illustrated that cluster 1 had more infiltration of immune cells, such as CD4 T cell, CD8 T cell, killer T cell, helper T cell 1, Follicular helper T cell, and Follicular helper T cell [Figure S6a and b]. In particular, iTreg was also highly infiltrated in cluster 1, which was consistent with the result of the previous analysis.

Distribution and clinical significance of Treg cell in PCa

Clinical data were collected from Shanghai East Hospital. To evaluate the differential infiltration proportion of Treg cell in PCa across disease progression stages, IHC for Foxp3 and CD45 was performed on paraffin-embedded sections from 44 paired LN metastases and 15 paired primary PCa samples. Although the IHC staining intensity for CD45 and Foxp3 was consistently moderate to strong (scores 2–3), their positive incidence varied significantly [Figure 6 and Figure S7]. Notably, most Foxp3-positive cells localized within CD45-positive regions, confirming Foxp3 as a specific marker for Tregs in T cells [Figure 6 and Figure S7]. IHC analysis revealed significantly higher Treg infiltration in the primary tumor, para-carcinoma, and stroma of PCa patients with LN metastasis compared to those without LN metastasis [Figure 6a]. Tregs were also more abundant in adjacent tissues and mesenchyme. In LN tissues, Treg infiltration was more frequent in metastatic LN than in non-metastatic LN. Moreover, negative LNs from metastasis-positive patients showed elevated Treg levels relative to LNs from metastasis-negative patients [Figure 6b]. In addition, as the progression of PCa (such as Gleason score and pathological T stage), Treg infiltration proportion in tumor tissue was also increased [Figure 6c]. These findings underscore the immunosuppressive role of Tregs and align with prior results.

Close

Figure 6:

IHC and statistical analysis based on clinical PCa samples. (a and b) IHC of Foxp3 based on primary and LN metastatic PCa tissue; (c) The different distribution of Treg cells in different location and stage of PCa. Scale bars: 500 μm (left); 100 μm (right). Foxp3: 1 (NC (The average number of cells per high-power field) <5 or PC (the percentage of positive cells) < 1%), 2 (5 ≤ NC < 10 or 1% ≤ PC < 5%), 3 (10 ≤ NC < 20 or 5% ≤ PC < 10%), 4 (20 ≤ NC < 50 or 10% ≤ PC < 25%), 5 (NC ≥ 50 or PC ≥ 25%). IHC: Immunohistochemical, LN: Lymph node, PCa: Prostate cancer, Treg: Regulatory T.

Export to PPT