The characteristics of eight melanoma patients, including three who developed irAEs (AE group) following ICI therapy, three who did not develop irAEs (NAE group) post-ICI therapy, and two control patients who did not receive ICI therapy (UNT group), are summarized in Table 1. The specific irAEs observed included nausea, vomiting and fatigue in one patient after eight cycles of pembrolizumab, skin rash and vitiligo in another after four cycles of toripalimab, and erythematous papules in a third patient after two cycles of toripalimab. Patients in non-irAEs group were all treated with toripalimab.

Single-cell RNA sequencing (scRNA-seq) was performed on peripheral blood mononuclear cells (PBMCs) to assess the composition and functional changes of various immune cells in melanoma patients undergoing anti-PD-1 therapy. After excluding low-quality cells, a total of 54,793 cells were analyzed and mapped to a multimodal PBMC reference dataset, and the results were visualized using Uniform Manifold Approximation and Projection (UMAP) plots (Figure 1A). The relative abundance of these broad immune cell categories across different patient groups is depicted in Figure 1B. We observed substantial alterations in the immune cell landscape: there was a decrease in the frequency of CD4 T cells and an increase in CD8 T cells in the AE group compared to both the NAE group and the UNT group (Figure 1C). Notably, the NAE group exhibited a higher frequency of NK cells and a lower frequency of monocytes than either the AE group or the UNT group (Figure 1C). Meanwhile, the frequency of DC was slightly higher in AE group, whereas the frequency of B cells remained comparable across all groups (Figure 1C). When the data were normalized to total T cells, the frequency of CD8 T cells was found to be highest in the AE group and lowest in the UNT group (Figures 1D, E). Conversely, the frequency of CD4 T cells was lowest in the AE group and highest in the UNT group, with the NAE group displaying intermediate frequencies (Figures 1D, E). To further assess the impact of anti-PD-1 treatment, we examined the expression levels of PD-1 and other immune checkpoint receptors. A dot plot analysis revealed a decreased expression of PDCD1 (PD-1) and an increased expression of CTLA4, TIGIT, LAG3, HAVCR2 (TIM-3), and the transcriptional factor TOX in T cells post-treatment compared to the UNT group (Figure 1F). These changes indicate the effectiveness of anti-PD-1 therapy in both the AE and NAE groups and suggest a compensatory regulation of immune checkpoint receptors during treatment, with these effects being more pronounced in the AE group. In summary, anti-PD-1 therapy was associated with a shift in the balance of T cell subpopulations within PBMCs, characterized by an increased frequency of CD8 T cells and a reduced frequency of CD4 T cells. This shift was most pronounced in the AE group, underscoring the potential link between T cell dynamics and the development of irAEs during therapy.

To determine the specific T cell subsets associated with the development of adverse events in melanoma patients undergoing anti-PD-1 therapy, all T cells were extracted for sub-analysis. UMAP visualization categorized T cells into 11 distinct clusters: three CD4 T cell clusters, three CD8 T cell clusters, T regulatory cells (Tregs), gamma delta T cells (gdT), Mucosal-associated invariant T cells (MAIT), proliferating T cells, and other unspecified T cells (Figure 2A). Using canonical markers, we defined CD4 T_1 and CD8 T_1 as naïve T cells (TCF7+), CD4 T_2 as activated memory CD 4+ T cells (TNFRSF4+), CD8 T_2 as effector-memory-like CD8+ T cells (GZMK+), and CD8 T_3 as terminal effector/cytotoxic CD8+ T cells (CX3CR1+). The CD4 T_3 cluster displayed a distinct profile, co-expressing platelet-related genes (e.g., PPBP, PF4) and a strong interferon-stimulated gene signature (e.g., ISG15, IFIT1). This pattern may reflect a state of tumor-educated, hyperactivated CD4+ T cells, though a contribution from potential cell doublets cannot be excluded. The expression and distribution of these markers are visualized as DotPlot and FeaturePlot in Figures 2G, H. The relative abundance and cell cycle distribution of each cluster were analyzed (Figures 2B, C). Notably, the AE group exhibited a slightly higher percentage of cells in the G2/M phase, indicating enhanced cellular proliferation. Consistently, the frequency of proliferating T cells was increased in the AE group compared to NAE group (Figure 2F). While overall T cell cycle activity shows a trend in NAE group, the clonal expansion of a dedicated proliferating T cell population is a distinctive feature of the AE group. Following anti-PD-1 therapy, the frequency of CD4 T_2 and CD8 T_2 cells decreased, while the frequency of CD8 T_3 cells increased in the AE group compared to both the untreated and NAE groups (Figures 2D, E). Meanwhile, the combined effector CD8 T cells (CD8 T_2 and CD8 T_3) were increased in AE group (Figure 2E). However, these changes did not reach statistical significance, likely due to the limited sample size. Though not statistically significant, the ratio of CD8 T_2 to CD8 T_3 cells was considerably decreased in the AE group compared to the UNT group (Figure 2E). Additionally, the frequencies of CD4_1 cells and Treg cells remained unchanged across all groups (Figures 2D, F). In summary, the analysis revealed that the frequency of CD8 T_3 (cytotoxic T) cells and proliferating T cells were elevated in the AE group. This pattern suggests a mechanistic link between the enhanced propensity of T cells toward a cytotoxic and proliferative phenotype and the development of irAEs in patients receiving anti-PD-1 therapy.

To further delineate the molecular characteristics of each T cell cluster, differentially expressed genes (DEGs) were calculated using FindAllMarkers function and heatmap was plotted with top 10 genes in each cluster (Figure 3A). The CD8 T_3 cluster demonstrated enrichment of cytotoxic T cell markers, including GNLY, FGFBP2, GZMH, GZMB, NKG7, and CX3CR1, indicating a terminally differentiated cytotoxic phenotype. These markers were also present in the CD8 T_2 cluster, albeit at lower levels (Figure 3A). Notably, CD8 T_2 cluster was marked by GZMK, CCL4, CCL5, DUSP2, and CD160, indicating a effector memory phenotype. To assess the functionality of these clusters, Gene Set Variation Analysis (GSVA) scores for pathways related to immunotherapy were calculated. Both CD8 T_2 and CD8 T_3 clusters displayed elevated cytotoxic scores compared to other clusters, with CD8 T_3 exhibiting the highest scores (Figures 3B, C). In contrast, CD8 T_2 was marked by the highest expression of GZMK+ transitional genes (PMID: 33271118), which was markedly lower in CD8 T_3 (Figures 3B, C). Notably, the exhaustion score was similar between CD8 T_2 and CD8 T_3. Additionally, both the PD-1 high and PD-1 low signatures were modestly enriched in CD8 T_2. However, CD8 T_3 exhibited the highest PD-1 low signature and the lowest PD-1 high signature, suggesting that CD8 T_3 was the primary cluster affected following anti-PD-1 therapy (Figure 3C). Upon comparing treatment groups, it was evident that anti-PD-1 therapy enhanced the expression of cytotoxicity and exhaustion-related pathways while diminishing the expression of pathways associated with GZMK+ transitional genes, naïve cell markers, and high PD-1 expression (Figure 3D). Given the established link between GZMK and immune aging (20), we assessed its correlation with chronological age in our cohort. Analysis revealed no significant association between patient age and either GZMK expression or the GSVA score for the GZMK-transitional gene signature in T cells (Supplementary Figures 1A, B). To further examine the biological differences between CD8 T_2 and CD8 T_3, additional DEGs were identified and the top 25 genes in each cluster were displayed in a heatmap (Figure 3E). Subsequent gene ontology (GO) analysis revealed that the CD8 T_2 was more enriched in pathways associated with T cell activation and differentiation (Figure 3F), whereas the CD8 T_3 was predominantly enriched in pathways linked to cytotoxicity (Figure 3G). Taken together, our findings highlight distinct molecular and functional profiles within the CD8 T cell subpopulations in melanoma patients receiving anti-PD-1 therapy.

To gain a deeper understanding of immune changes in CD8 T cells during anti-PD-1 therapy, all CD8 T cells were extracted for further analysis. UMAP visualization categorized these cells into five clusters, naïve, transitional (GZMK+), cytotoxic_1, cytotoxic_2, and proliferating, the former three corresponding to CD8 T_1, CD8 T_2, and CD8 T_3 in Figures 2A, respectively (Figure 4A). Violin plots displaying GSVA scores revealed that transitional (GZMK+) cluster represents an intermediate state, showing moderate expression of naïve and cytotoxic signatures and a peak in the GZMK-transitional signature (Figure 4B). Following anti-PD-1 therapy, a substantial increase in cytotoxic_1 CD8+ T cells and a corresponding decrease in naïve and transitional (GZMK+) cells were observed in the AE group, meanwhile the low-frequency cytotoxic_2 subset followed a similar trend with cytotoxic_1 subset. However, these changes did not reach statistical significance due to the limited sample size (Figures 4C, D). Additionally, the ratio of transitional-to-cytotoxic T cells was considerably decreased in the AE group compared to the UNT group (Figure 4D). In contrast, the NAE group showed a slightly higher percentage of naïve CD8 T cells and a comparable percentage of cytotoxic cells relative to the UNT group (Figure 4D). Additionally, the AE group demonstrated increased proliferative capacity compared to both the NAE and UNT groups (Figure 4D). To delineate differences among the groups, DEGs were identified and the top 10 genes in each group were visualized through a heatmap (Figure 4E). Further GO analysis showed that the AE group was specifically enriched in cell-killing-associated pathways, while the UNT and NAE groups were more associated with metabolic and cell-activation pathways (Figure 4F). In summary, these results further indicate a marked increase in cytotoxic CD8+ T cells in the AE group following anti-PD-1 treatment, accompanied by a decreased ratio of transitional to cytotoxic CD8 T cells.

To elucidate the dynamics of CD8 T cell subpopulations during anti-PD-1 therapy, we conducted a trajectory analysis on CD8 T cells. This analysis identified two principal differentiation pathways originating from naïve CD8 T cells: one leading to transitional CD8 T cells and the other progressing to terminal cytotoxic CD8 T cells, as depicted in Figures 5A–C. Transitional CD8 T cells are predominantly located at the midpoint of this trajectory, indicating a crucial transitional phase prior to their full differentiation into the cytotoxic phenotype. Distribution analysis of these cells across the identified states revealed that naïve CD8 T cells were primarily found in state 8, whereas transitional CD8 T cells were concentrated in states 7 and 8, and cytotoxic CD8 T cells were distributed across other states, as shown in Figure 5D. To visualize the dynamic shifts in gene expression associated with these subpopulations, we plotted the expression levels of key genes over the pseudotime continuum in Figure 5E. Here, expression of naïve and cytotoxic markers peaked at the respective ends of the pseudotime axis, while markers such as GZMK, CMC1 and DUSP2 exhibited peaks during the transitional phase, suggesting a gradual shift in cellular function. Further refinement of this analysis was achieved through BEAM (Branch Expression Analysis Modeling), which identified genes differentially expressed along these differentiation paths. The cell fate toward transitional CD8 T cells was characterized by enrichment of genes such as MALAT1, NEAT1, SYNE2, CMC1, SH2D1A, and GZMK (Figure 5F), suggests a state of cellular preparation for a functional shift. Conversely, the trajectory toward cytotoxic CD8 T cells was characterized by enrichment of genes such as NKG7, GZMB, GNLY, PRF1, and FGFBP2 (Figure 5F), which are well-known for their roles in cell-mediated cytotoxicity. These findings indicate that transitional CD8 T cells may serve as a versatile intermediary, capable of modulating the immune response following anti-PD-1 therapy. This modulation may influence the balance between a regulated immune response and an aggressive cytotoxic attack, potentially mitigating or exacerbating the risks of irAEs.

Peripheral blood samples were collected from patients with melanoma admitted to the second affiliated hospital of Zhejiang University, School of Medicine. This study included patients with resectable stage II–III melanoma, as staged by the American Joint Committee on Cancer (AJCC, 8th edition). Six patients received adjuvant anti-PD-1 therapy, three of whom developed immune-related adverse events (irAEs), and two patients were under observation without treatment. All participants provided written informed consent, and the study protocol was approved by the local institutional ethics committee. Routine blood tests and clinical examinations were performed at regular intervals to monitor for early signs of irAEs. Peripheral blood mononuclear cells (PBMCs) were isolated from whole blood via density gradient centrifugation. For patients experiencing irAEs, PBMCs were collected within one week of irAE onset. For those without irAEs, samples were collected at a comparable time point following the initiation of anti-PD-1 therapy.

Cells were loaded onto the 10X Chromium Single Cell Platform (10X Genomics) at a concentration of 1,000 cells per μl (Single Cell 3′ library and Gel Bead Kit v.3) as described in the manufacturer’s protocol. Generation of gel beads in emulsion (GEMs), barcoding, GEM-RT clean-up, complementary DNA amplification and library construction were all performed as per the manufacturer’s protocol. Qubit was used for library quantification before pooling. The final library pool was sequenced on the Illumina HiSeq instrument using 150-base-pair paired-end reads.

Raw sequencing data in FASTQ format were processed using CellRanger (10x Genomics) to generate digital expression matrices. Low-quality barcodes and cells with fewer than 1,000 UMIs were excluded. Subsequent analyses were performed using the Seurat R package (v4.1.0) (40). Cells with fewer than 600 genes, more than 4,500 genes, or more than 10% mitochondria content were filtered out. Ribosomal genes were not filtered out pre-normalization. Data were normalized and integrated using the SCTransform method. Cell type annotation was performed via reference mapping using the MapQuery function according to Seurat v4 Vignettes (40). Then cells were split according to celltype.l1 for subpopulation analysis. Data were integrated for batch correction before further analysis. RunPCA, FindNeighbors, FindClusters, and RunUMAP were used to calculate the dimension-reduction coordinates for visualization and to perform unsupervised clustering. We used the first 30 principal components (PCs) for downstream analysis. For FindNeighbors and RunUMAP, a dim parameter of 1:30 was used. For FindClusters, a resolution parameter of 0.2 and 0.1 was used in all T cells and CD8 T cells sub-analysis, respectively. Differential expression analysis was conducted with the FindAllMarkers function. Gene Set Enrichment Analysis (GSEA) was performed using the clusterProfiler (v4.2.2) and enrichplot (1.14.2) packages. Pseudotime trajectory analysis was conducted with Monocle2 (v2.20.0). Prior to analysis, cells were downsampled to equal numbers across samples, and genes were ordered based on dispersional and empirical criteria. Dimensionality reduction for trajectory visualization was achieved using the DDRTree method.

Cell population proportions were compared using non-parametric test (Mann-Whitney test) due to limited sample size. Bar charts were generated using GraphPad Prism (La Jolla, CA, USA). A p-value of less than 0.05 was considered statistically significant.