Single-cell transcriptomes revealed that TRMs are more immunoreactive than non-TRMsWe obtained scRNA-seq data from 38,572 T cells from 11 CC and 5 normal cervical samples (Fig. 1a). T cells were classified as TRMs (12,945 cells; ITGAE+ CD3D+) or non-TRMs (25,627 cells; ITGAE− CD3D+) on the basis of ITGAE expression (Fig. 1b), consistent with TRMs described in previous researches21,22. In addition, we compared the expression of ITGAE in TRMs, non-TRMs, and T cells in peripheral blood mononuclear cells (PBMC), further validating the specificity of ITGAE for defining TRMs (Supplementary Fig. 1d). Differentially expressed genes (DEGs) analysis revealed that TRMs expressed higher levels of genes associated with tissue residence (TGFB1 and ITGB7), cytotoxicity (GZMB and PRF1), and chemokine production (CXCL13 and CCL5) than non-TRMs (Fig. 1c, d). The gene ontology (GO) term analysis showed that DEGs related to the cell adhesion pathway were enriched in TRMs than non-TRMs, which is associated with the retention of T cells in tissues. In addition, pathways associated with T cell activation, lymphocyte-mediated immunity, and antigen processing and presentation were also enriched in TRMs than non-TRMs (Fig. 1e, f and Supplementary Fig. 1c). Moreover, TRMs were more abundant in CC tissues than in normal cervical tissues (p = 0.0053; Fig. 1g, h). Subsequently, we confirmed that the infiltration degree of TRMs in CC was higher than that in normal cervical tissues by performing immunofluorescent labeling on formalin-fixed and paraffin-embedded (FFPE) samples of CC and normal cervical tissues (p = 0.0113; Fig. 1i). These findings indicate that TRMs are more abundant and immunoreactive than non-TRMs in the TME of CC.Fig. 1: TRMs are more immunoreactive than non-TRMs.a Overview of the experimental workflow. b t-Distributed Stochastic Neighbor Embedding (tSNE) plots showing subcluster, sample origin, and marker gene expression. The color key shows the gradient of normalized expression. c Volcano plot displaying differentially expressed genes (DEGs) between TRMs (red) and non-TRMs (blue). d Violin plots illustrating the expression of the indicated genes in TRMs (red) and non-TRMs (blue). *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (two-sided Wilcoxon test). e Gene ontology (GO) term analysis in TRMs. The intensity of color indicates p-value. f Gene Set Enrichment Analysis (GSEA) showing that TRMs were enriched in gene sets associated with the immune response. NES, normalized enrichment score. g Bar plots depicting the proportions of TRMs (red) and non-TRMs (blue) in each sample. h Box plots illustrating the proportions of cell subclusters in normal cervix and cervical cancer (CC) tissues. i Representative immunofluorescent labeling of CD103 (green), CD3 (red), and 4,6-diamidino-2-phenylindole (DAPI, blue) in sections from CC (up) and normal cervical tissues (down). Box plots showing the percentage of CD103+ CD3+ TRMs and CD103- CD3+ non-TRMs in CC and normal cervical tissues based on immunofluorescent labeling results. Scale bar in figure, 80 μm. The p-values were generated using the Wilcoxon test.Single-cell transcriptome analysis revealed that TRMs are more activated in CC tissues than in normal cervical tissuesBased on the higher infiltration of TRMs in CC tissues than in normal ones, we next explored the transcriptomes of TRMs in CC tissues and in normal cervical tissues (Fig. 2a). We first identified DEGs between TRMs from CC tissues and those from normal cervical tissues (Fig. 2b). Functional gene set analysis revealed that TRMs from CC tissues not only expressed more genes related to cytotoxicity and proliferation, but also expressed higher levels of genes encoding immune checkpoint molecules than TRMs from normal cervical tissues. By contrast, the TRMs from normal cervical tissues had a more naïve-like T cell signature (Fig. 2c). The GO term analysis of DEGs showed that TRMs in CC tissues were enriched in pathways related to T cell activation, antigen processing and presentation of peptide antigen via MHC class II, and lymphocyte-mediated immunity. They were also enriched in pathways associated with defense response to virus (Fig. 2d and Supplementary Fig. 2). The gene set enrichment analysis (GSEA) results confirmed the enrichment of gene sets related to T cell activation, antigen processing and presentation of peptide antigen via MHC class II, as well as the response to type I interferon and defense response to virus (Fig. 2e). These results suggest that TRMs residing in the TME of CC are more activated than those localized to normal cervical tissues; moreover, TRMs may play a crucial role in the anti-viral immune response in the cervix.Fig. 2: Characteristics of TRMs in the TME of CC.a tSNE plots showing distribution of TRMs from CC tissues (blue) and those from normal cervical tissues (pink). b Volcano plot showing DEGs between TRMs from CC tissues and those from normal cervical tissues. c Heatmap showing the expression of functional signature genes, whereby color intensity indicates average gene expression. d GO term analysis in TRMs from CC tissues. The color intensity indicates p-value. e GSEA showing differences in the activation of various immune cell pathways in TRMs from CC tissues and in those from normal cervical tissues.Single-cell transcriptome analysis of TRMs heterogeneityTo further dissect the molecular heterogeneity of TRMs, we extracted data from all the TRMs for use in a sub-clustering analysis. We obtained seven TRM subclusters: six CD8+ TRM subclusters, differentiated by CXCL13, PLAC8, ITM2C, STMN1, FCER1G, or XIST expression; and one CD4+ TRM cluster (Fig. 3a). These TRM subclusters exhibited widely divergent gene expression profiles (Fig. 3b and Supplementary Fig. 3a). Among them, CXCL13+ CD8+ TRMs demonstrated elevated expression of genes encoding inhibitory molecules (LAG3, TIGIT, PDCD1, HAVCR2, and CTLA4) and had a higher inhibitory score than other TRM subclusters. Additionally, CXCL13+ CD8+ TRMs also expressed genes encoding cytotoxic molecules (GZMB and PRF1; Fig. 3c and Supplementary Fig. 3b). In addition, PLAC8+ CD8+ TRMs demonstrated elevated expression of genes encoding cytotoxic molecules (GZMA, GNLY, NKG7, and PRF1; Fig. 3c) and had a higher naïve score than other CD8+ TRM subclusters (Supplementary Fig. 3b). The relative abundances of the seven major TRM subclusters in the CC tissues and normal cervical tissues were shown in Fig. 3d. Further comparisons revealed that the proportion of infiltrating CXCL13+ CD8+ TRMs was significantly higher in CC tissues than in normal cervical tissues (p = 0.0022), whereas the proportion of infiltrating PLAC8+ CD8+ TRMs was significantly lower in CC tissues than in normal cervical tissues (p = 0.00092; Fig. 3e).Fig. 3: Single-cell transcriptome profiles of 12,945 TRMs.a tSNE plots showing the TRM subclusters and the expression of marker genes; the color key shows the gradient of normalized expression. Heatmaps displaying the relative expression of the top 3 DEGs (b) and functional signature genes (c) in the seven TRM subclusters; the color intensity reflects the average level of gene expression. d The proportions of the seven TRM subclusters in each tissue and type (CC tissues vs. normal cervical tissues), colored by cell types. e Box plots illustrating the proportions of TRM subclusters in normal cervical (pink) and CC (red) tissues. The p-values were calculated using the paired Wilcoxon test.
CXCL13
+ CD8+ TRMs are enriched in CC tissues and express transcripts related to antigen processing and presentation and defense responses to the virusTo determine the unique molecular signature of CXCL13+ CD8+ TRMs, we identified DEGs between CXCL13+ CD8+ TRMs and the other TRM subclusters. The result exhibited that CXCL13+ CD8+ TRMs showed high expression of cytotoxic genes (GZMB, GNLY, and PRF1) and inhibitory genes (HAVCR2, TIGIT, LAG3, CTLA4, and PDCD1) (Fig. 4a, b). The GO term analysis revealed that CXCL13+ CD8+ TRMs were enriched in pathways associated with T cell activation, cytokine production, as well as response to virus and response to type I interferon (Fig. 4c). Moreover, CXCL13+ CD8+ TRMs expressed higher level of genes related to defense response to virus, including IFI6, IFIT3, BST2, MX1, PLSCR1, ISG15, OAS1, STAT1, HERC5, and IRF7, than the other TRM subclusters (Fig. 4d). The GSEA results also confirmed that CXCL13+ CD8+ TRMs were enriched in gene sets associated with defense response to virus and response to type I interferon (Fig. 4e and Supplementary Fig. 4). These findings suggest that CXCL13+ CD8+ TRMs may exert anti-viral effects upon stimulation by type I interferon.Fig. 4: Characteristics of CXCL13+ CD8+ TRMs.a tSNE plot showing the distribution of CXCL13+ CD8+ TRMs (red) relative to that of the other TRM subclusters (blue). b Volcano plot depicting DEGs between CXCL13+ CD8+ TRMs (red) and the other TRM subclusters (blue). c GO term analysis of CXCL13+ CD8+ TRMs. The color intensity indicates p-value magnitude. d Violin plots illustrating the expression of the indicated genes in CXCL13+ CD8+ TRMs and the other TRM subclusters. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (two-sided Wilcoxon test). e GSEA showing differences in pathway activity between CXCL13+ CD8+ TRMs and the other TRM subclusters. f Volcano plot depicting the DEGs between CXCL13+ CD8+ TRMs from CC tissues (red) and those from normal cervical tissues (blue). g Violin plots displaying the expression of the indicated genes in CXCL13+ CD8+ TRMs from CC tissues and in those from normal cervical tissues. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (two-sided Wilcoxon test). h GO term analysis of CXCL13+ CD8+ TRMs from CC tissues; color intensity indicates p-value magnitude. i GSEA results showing differences in the pathway activities of CXCL13+ CD8+ TRMs from CC tissues and those from normal cervical tissues. j Kaplan–Meier survival analysis of patients with CC following radiotherapy from TCGA database, stratified based on high or low signature gene set score of CXCL13+ CD8+ TRMs. The p-values were calculated using the two-sided log-rank test.Because CXCL13+ CD8+ TRMs were predominantly found in CC tissues, we explored the differences in the signatures of CXCL13+ CD8+ TRMs from CC tissues and normal cervical tissues. DEGs analysis revealed that MHC class II molecule genes (HLA-DRB5 and HLA-DRA) and chemokine-related genes (CCL4L2 and CCL5) were highly expressed in CXCL13+ CD8+ TRMs from CC tissues (Fig. 4f, g). The GO term analysis revealed that CXCL13+ CD8+ TRMs from CC tissues were enriched in pathways related to antigen processing and presentation, cell adhesion, and chemotaxis (Fig. 4h); this was confirmed by the results of the GSEA (Fig. 4i). We therefore next evaluated the relationship between CXCL13+ CD8+ TRM signature and the survival outcomes of patients with CC following radiotherapy using bulk RNA-seq data obtained from The Cancer Genome Atlas (TCGA). We found that the signature gene set score of CXCL13+ CD8+ TRMs was positively correlated with the overall survival (OS) and progression-free survival (PFS) of patients with CC following radiotherapy (Fig. 4j).In summary, these findings indicate that CXCL13+ CD8+ TRMs exhibit cytotoxic and inhibitory characteristics, enriched in pathways associated with defense response to virus and the type I interferon response than the other TRM subclusters, and are associated with the improved outcome of patients with CC following radiotherapy.
PLAC8
+ CD8+ TRMs are less enriched in CC tissues than in normal cervical tissues but exhibit highly cytotoxic characteristicsWe next dissected the molecular signature of PLAC8+ CD8+ TRMs, the subcluster which was less abundant in CC tissues than in normal cervical tissues (Fig. 5a). DEGs analysis exhibited that PLAC8+ CD8+ TRMs showed high expression of cytotoxic genes (GZMH, GNLY, GZMM, KLRD1, KLRF1, KLRG1, NKG7, and PRF1) and low expression of genes encoding inhibitory molecules (PDCD1, CTLA4, TIGIT, and HAVCR2) (Fig. 5b, c), exhibiting similar characteristics to cytotoxic T cells. Further GO term analysis revealed that PLAC8+ CD8+ TRMs were enriched in pathways associated with T-cell-mediated immune effector processes, and also enriched cytoplasmic translation and ribosome biogenesis pathways (Fig. 5d, e).Fig. 5: Characteristics of PLAC8+ CD8+ TRMs.a tSNE plot showing the distribution of PLAC8+ CD8+ TRMs (red) relative to that of the other TRM subclusters (blue). b Volcano plot showing DEGs between PLAC8+ CD8+ TRMs (red) and the other TRM subclusters (blue). c Violin plots displaying the expression of the indicated genes between PLAC8+ CD8+ TRMs and the other TRM subclusters. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (two-sided Wilcoxon test). d GO term analysis of PLAC8+ CD8+ TRMs; the color intensity indicates p-value magnitude. e Violin plots showing differences in pathway activity between PLAC8+ CD8+ TRMs (red) and the other TRM subclusters (blue). f Volcano plot showing DEGs in PLAC8+ CD8+ TRMs between CC tissues (red) and normal cervical tissues (blue). g Violin plots illustrating the expression of indicated genes in PLAC8+ CD8+ TRMs between CC tissues and normal cervical tissues. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 (two-sided Wilcoxon test). h GO term analysis in PLAC8+ CD8+ TRMs of CC tissues. The intensity of color indicates p-values. i GSEA showing differences in the pathway activities of PLAC8+ CD8+ TRMs from CC tissues and those from normal cervical tissues. j Kaplan–Meier survival analysis of patients with CC following radiotherapy from TCGA database, stratified based on high or low signature gene set score of PLAC8+ CD8+ TRMs. The p-values were calculated using a two-sided log-rank test.In addition, we explored signature differences between PLAC8+ CD8+ TRMs from CC tissues and those from normal cervical tissues. DEGs analysis revealed that HLA-DRB5, CD74, HLA-DRA, and HLA-DQA2 were more highly expressed in PLAC8+ CD8+ TRMs from CC tissues than in those from normal cervical tissues (Fig. 5f, g). The GO term analysis demonstrated that PLAC8+ CD8+ TRMs from CC tissues were enriched with immune response-related pathways, including T cell activation, cytokine-mediated signaling pathway, antigen processing and presentation, and cell adhesion (Fig. 5h and Supplementary Fig. 5). The GSEA confirmed that PLAC8+ CD8+ TRMs were enriched in gene sets associated with antigen processing and presentation, as well as cell adhesion (Fig. 5i). Moreover, survival analysis revealed a positive correlation between the signature gene set score of PLAC8+ CD8+ TRMs and the OS and PFS of patients with CC following radiotherapy (Fig. 5j). These findings suggest that although PLAC8+ CD8+ TRMs were not abundant in the TME of CC, they were a highly cytotoxic, readily activated, clinically important TRM subset.
CXCL13
+ CD8+ TRMs and PLAC8
+ CD8+ TRMs interact with epithelial cells in the TME of CCTo further investigate the interaction between epithelial cells and TRMs, we used CellChat to interrogate the scRNA-seq data and gather evidence of cell-cell communication through the analysis of ligand-receptor complexes. We found that the interactions between epithelial cells and TRMs were more robust in CC tissues than in normal cervical tissues, both in terms of their number and strength (Fig. 6a). Specifically, epithelial cells, as senders, showed more active interactions with CXCL13+ CD8+ TRMs and PLAC8+ CD8+ TRMs than the other TRM subclusters in CC tissues, both in terms of number and strength (Fig. 6b and Supplementary Fig. 6a).Fig. 6: Crosstalk between epithelial cells and TRMs.a Bar graphs displaying the number and strength of interactions between epithelial cells and TRMs in normal cervical tissues and CC tissues, respectively. b CellChat analysis of scRNA-seq data to identify cell-cell communications between epithelial cells and TRM subclusters in CC tissues. The color intensity indicates the number (top panel) and strength (bottom panel) of the cell-cell interactions. c Bubble plots showing differences in ligand-receptor pairs involved in interactions between epithelial cell-CXCL13+ CD8+ TRM and epithelial cell-PLAC8+ CD8+ TRM in normal cervical tissues and CC tissues. Bubble size and color intensity indicate p-value magnitude and communication probability, respectively. d Circle plots showing the indicated ligand-receptor pairs involved in interactions between epithelial cells and either CXCL13+ CD8+ TRMs, or PLAC8+ CD8+ TRMs in CC tissues (left) and normal cervical tissues (right); color-coded by cell type. e Representative immunofluorescent labeling of Pan-CK (yellow), CXCL13 (red), CD103 (green), CD8 (white), and DAPI (blue) to identify epithelial cells, CXCL13+ CD8+ TRMs, and their spatial relationship in CC sections. Scale bar in figure, 20 μm.Next, we explored the specific ligand-receptor pairs which participated in the interactions between epithelial cells and either CXCL13+ CD8+ TRMs or PLAC8+ CD8+ TRMs both in CC tissues and in normal cervical tissues (Fig. 6c). We found that MHC class I molecules-related pairs (HLA-A/B/C/E/F – CD8A/ CD8B), T cell proliferation and activation-related pairs (MIF-(CD74+CD44)(CD74+CD44) and MDK-NCL) were more active both in epithelial cell-CXCL13+ CD8+ TRM and epithelial cell-PLAC8+ CD8+ TRM at CC tissues than normal cervical tissues (Fig. 6c and Supplementary Fig. 6b). This result suggested that tumor antigens may be presented to TRMs via MHC class I molecules in the TME of CC, leading to the rapid activation of this T cell subset and the induction of anti-tumor immune responses. In addition, HLA-E – KLRC1 and HLA-E – CD94:NKG2A pairs were enriched in epithelial cell-CXCL13+ CD8+ TRM at CC tissues, rather than in epithelial cell-PLAC8+ CD8+ TRM (Fig. 6c, d). Of note, both of these interactions are associated with T cell inhibition and increased cancer growth23,24. We confirmed the presence and spatial relationship of epithelial cells and CXCL13+ CD8+ TRMs within the TME of CC by multiplex immunofluorescent staining by labeling for pan-cytokeratin (Pan-CK), CXCL13, CD103, and CD8. The proximity observed between epithelial cells and CXCL13+ CD8+ TRMs on tissue sections supports their potential cellular interactions (Fig. 6e).In summary, we detected active interactions between epithelial cells and either CXCL13+ CD8+ TRMs, or PLAC8+ CD8+ TRMs in the TME of CC. Moreover, the immunosuppressive HLA-E – KLRC1 and HLA-E – CD94:NKG2A pairs were enriched in CXCL13+ CD8+ TRMs interacting with epithelial cells at CC tissues.
Cellular heterogeneity and key subsets of tissue-resident memory T cells in cervical cancer
