The activity of base excision repair is positively correlated with the infiltration of CD4+ T cells in melanoma

scRNA-seq data analysis

Raw count matrix of GSE222446 was downloaded from GEO (https://www.ncbi.nlm.nih.gov/gds) [24]. The Seurat object was constructed by combining the read count matrix of each patient using the R package ‘Seurat’ (version 5.2.1) [25]. Only genes expressed in at least 1 % of cells were retained for subsequent analysis. Each patient was considered as a batch effect. Normalization was performed via natural-log transformation with a scale factor of 10,000. Genes with the highest variance (n=2,000) were used to calculate principal components (PCs). Harmony integration was applied to obtain the PC reduction space corrected for batch effects. The original Seurat clusters were generated at a resolution of 0.2 using the first 20 batch-corrected PCs. The melanoma signature [26] and ‘AUCell’ analysis (version 1.20.2) [27] were used to assist in identifying the melanoma population. ‘scType’ (Code copyright 2021) [28] was employed to identify cell types using canonical cell type markers from the CellSTAR database (https://cellstar.idrblab.net/) [29].

To further identify melanoma subclusters, a total of 53,112 melanoma cells were initially clustered into 6 subclusters. Hallmark gene sets from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) [30] were used to characterize each subcluster with AUCell analysis. Finally, four melanoma subtypes were identified based on the average AUCell scores in the original subclusters. A total of 16,551 T cells were further integrated into 11 original subclusters at a resolution of 0.6 using the first 20 PCs, and eight subpopulations were finally identified, employing an approach similar to that used for major cluster identification.

Measurement for the relative abundance of melanoma cell subtypes

To elucidate the association between specific melanoma cell subtypes and tumor immune microenvironment features in bulk RNA-seq samples with comprehensive clinical annotation, an algorithm was developed to measure the relative abundance of a given melanoma subtype among all melanoma cells, leveraging scRNA-seq data. The FindAllMarkers function in the ‘Seurat’ package was used to identify marker genes for each of the four melanoma subtypes against all other major cell types and the other melanoma subtypes. Any marker genes overlapping with the melanoma signature from the article were removed [26]. The top 20 marker genes with the highest significance were selected for each melanoma subtype. The relative abundance of each melanoma subtype was calculated as the ratio of the average expression value of the marker genes for that subtype to the average expression value of the melanoma signature from the article [26].

TCGA SKCM cohort and establishment of immune subtypes

Clinical and survival data, along with the bulk RNA-seq Transcripts Per Million (TPM) matrix from The Cancer Genome Atlas Skin Cutaneous Melanoma (TCGA SKCM), were downloaded from UCSC (https://xenabrowser.net/datapages). The log2 (TPM) format was used to infer immune cell abundance using MCPcounter (version 1.2.0), as well as the relative abundance of melanoma cell subtypes. Immune subtyping of the TCGA SKCM cohort was conducted based on the abundance of 10 types of immune cells inferred from MCPcounter, which was normalized through z-transformation across samples. Three immune subtypes were identified using the ExecuteCC function in the R package ‘CancerSubtypes’ (version 1.17.1), with parameters clusterAlg=“km” and distance=“euclidean”. In order to evaluate the reproducibility of the immune subtypes, two other independent melanoma cohorts were used. GSE65904 consisted of 214 melanoma tumor biopsies profiled with Illumina HumanHT-12 V4.0 expression beadchip. GSE190113 was the bulk RNA-seq dataset that included 105 melanoma samples. The expression matrix corrected for batch effect with ComBat of GSE190113, and the normalized expression matrix of GSE65904 obtained through the R package limma (version 3.54.1) were applied to the MCP counter to perform immune subtyping with the same method as TCGA SKCM. Twenty-three canonical pathways from the article were used as a gene set applied to gene set variation analysis (GSVA) [31]. The average value of the normalized transformation of enrichment score within each immune subtype was visualized with a heatmap.

Two gene sets, the GOBP BER pathway and the KEGG BER pathway from MSigDB c5 and c2 gene sets, respectively, were applied to GSVA (version 1.48.3). The abundance of 18 kinds of immune scores was obtained with function consensus TMEA analysis setting parameter cancer=“SKCM”, StatMethod=“ssgsea” by using the R package consensus TME (version 0.0.1.9000). Tumor purity, established by four independent methods for TCGA SKCM, was obtained from the article [32], and the median value of the four tumor purity values was used for correlation analysis.

Massachusetts General Hospital skin pathology (MGSP) cohort analysis

The clinical data and tumor tissue bulk RNA-seq TPM matrix of the melanoma cohort that underwent immunotherapy were obtained from the supplementary materials accompanying the published article [33]. The calculation for the relative abundance of melanoma subtypes and the abundance of 18 types of immune cells followed the same methodology as used for the TCGA SKCM cohort. The immune score and the abundance of CD4+ T cells referred from consensus TME were used to categorize patients into low and high groups based on the 25 % percentile value, respectively. The immune score deduced by consensus TME is the measurement for total immune infiltration in any given tumor tissue sample.

Immunohistochemical (IHC) staining

The tissue microarray comprising 89 melanoma samples and their clinical data was retrieved from Shanghai Outdo Biotech Company in accordance with relevant regulations. The study was conducted in compliance with the Declaration of Helsinki. Antibodies used in the IHC assay: Anti-CD4: Clone EPR6855 (Epitomics/Abcam, Cambridge, UK), diluted at 1:100. Anti-CD8: Clone C8/144B (Dako/Agilent, Santa Clara, CA, USA), diluted at 1:100. Slides were scanned using a Leica SCN400 slide scanner (Leica Biosystems, Nussloch, Germany) and visualized with Digital Image Hub software (v2.1). Quantitative image analyses of CD4 and CD8 were performed using the “Measure Stained Cells Algorithm” (Leica Biosystems), with parameters optimized for chromogen-specific staining intensity, nuclear counterstaining thresholds, and cellular morphology. CD4:>10 % (high), 5–10 % (moderate), <5 % (low); CD8:>20 % (high), 10–20 % (moderate), <10 % (low). The average value from up to four fields of view per patient was calculated to represent overall marker expression. The expression of APE1 was scored on a semi-quantitative scale (0, 1, 2, 3), with each sample evaluated by three experienced pathologists.

This study was conducted in accordance with the guidelines of the Helsinki Declaration and was approved by the Ethics Committee of Army Medical Center of PLA (2022 number 51).

Statistical analysis

The association between the relative abundance of melanoma subtypes and the abundance of immune cells was determined by univariate and multivariate linear regression in the TCGA SKCM cohort by using the relative abundance of each of four melanoma cell types as an independent variable and enrichment score from consensus TME as a dependent variable. The clinical factors that were adjusted for in the multivariate linear regression included T stage (T3 vs. T0-1; T4 vs. T0-2), N stage (N1-3 vs. N0), and tumor biopsy sites (lymph node vs. others), which may impact the composition of the TME. The regression coefficients of the melanoma subtypes or correlation coefficients were visualized with a heatmap without any transformation. Correlation analysis was performed using the Spearman correlation coefficient with the R package ‘corrplot’ (version 0.92). Kruskal-Wallis test was used to compare the percentage of CD4 and CD8 positive staining or the ordinary variables of categorized CD4 and CD8 staining density among APE1 groups with correction for multiple comparisons. Kaplan-Meier curves and the log-rank test were used to evaluate the difference in survival rates between groups. Univariate and multivariate Cox regression were used to identify independent prognostic factors for PFS and OS. All analyses were performed with R version 4.4.3 (R Foundation for Statistical Computing, Vienna, Austria).

CD4+ T cell abundance correlates with optimal PFS and OS in melanoma patients receiving immunotherapy

To assess the impact of CD4+ T cells in melanoma, the MGSP dataset, comprising 121 melanoma samples from patients actually treated with PD-1 inhibitors, was analyzed [23]. The clinical characteristics are summarized in Supplementary Table S1 At the cutoff value of 25 % percentile value for CD4+ T cells, Kaplan-Meier plots showed that the PFS and OS of patients with high CD4+ T cell infiltration were significantly better than those with low CD4+ T cell infiltration (log-rank p<0.001, p=0.012, respectively) (Figure 1A and B). Univariate Cox regression showed that patients with high immune score or with high infiltration of CD4+ T cells have significantly decreased progression risk compared with patients with low immune score or low infiltration of CD4+ T cells (Supplementary Table S2). Moreover, the abundance of CD4+ T cells remained significant in multivariate Cox regression after being adjusted for gender, ECOG, stage, biopsy site, ICB drugs, ICB therapy line, BOR, and immune score (Supplementary Table S2). However, although high infiltration of CD4+ T cells was significantly associated with lower death risk in the univariate Cox regression, this association had borderline significance in the multivariate Cox regression (Supplementary Table S2). These results indicate that CD4+ T cell infiltration may independently influence the long-term benefit from ICB treatment in melanoma patients [23].

Figure 1:

Kaplan-Meier plots for melanoma patient groups from the MGSP cohort. (A) Poor progression-free survival (PFS) of patients with high or low CD4+ T cell infiltration in the MGSP cohort. (B) Overall survival (OS) of patients with high or low CD4+ T cell infiltration in the MGSP cohort. The high and low CD4+ T cell infiltration groups were determined by the 25 % percentile value of the enrichment score of CD4+ T cells from consensus TME in the MSGP cohort.

BER-active melanoma subtype positively correlates with CD4+ T cell infiltration

Melanoma exhibits marked heterogeneity across molecular subtypes, each displaying diverse clinical behaviors and responses to ICBs. To dissect subtype-specific immune interactions and evaluate the impact of DNA damage repair in immune response, we analyzed a single-cell RNA-seq dataset (GSE222446), which included 13 melanoma samples. Using ‘Seurat’ clustering corrected for batch effect, seven major cell populations were identified: melanoma cells, B cells, endothelial cells, epithelioid cells, fibroblasts, monocytes/macrophages, and T cells (Figure 2A). Focusing on T cells, only two CD4+ T cell subpopulations, including regulatory T cells and follicular helper CD4+ T cells, were identified. Then six CD8+ T cell subpopulations were characterized (Figure 2B). Melanoma cells were stratified into four molecular subtypes: BER type, inflammation type, lipid metabolism type, and myelocytomatosis oncogene (MYC) type (Figure 2C). As shown in Figure 2D, distinct biological features were identified in terms of the activity of canonical pathways specific to each subtype. For example, the BER type melanoma is active in DNA damage and repair pathways, including E2F targets, the G2M checkpoint, and the BER pathways. The cell number for each cell type was shown in Supplementary Table S3. The marker genes for these four subtypes of melanoma cells are shown in Supplementary Table S4.

Figure 2:

Identification of major cell types and subpopulations of melanoma cells and T cells in the GSE222446 dataset. (A) UMAP view of seven major cell types. (B) UMAP of eight T cell subsets. (C) UMAP of four melanoma subtypes. (D) heatmap representative of different biological characteristics of four melanoma subtypes. (E) scatter plot showing positive correlation of the proportion of CD4+T cells and BER type melanoma.

Correlationship analysis demonstrated a positive association between proportion of CD4+ T cells within all T cells and the proportion of BER type melanoma (p=0.047, Figure 2E), but no significant association was observed between proportion of CD4+ T cells and the proportion of inflammation type, lipid metabolism type, or MYC type melanoma cells (Supplementary Figure S1A–C). These findings suggest that the BER pathway in highly activated melanoma cells may promote CD4+ T cell infiltration, potentially serving as biomarkers for prognosis and therapeutic response prediction.

Validating the relationship between BER-active melanoma and CD4+ T cell infiltration in TCGA SKCM

To further confirm the findings regarding the positive correlation between CD4+ T cell infiltration and the BER pathway, as well as patient prognosis, melanoma samples from the TCGA database were utilized. The clinical characteristics are summarized in Supplementary Table S5. Firstly, we conduct a correlation analysis in order to verify the feasibility of the method to evaluate the relative abundance of each melanoma subtype in bulk RNA-seq data. The results showed that BER activity and APE1 expression were positively correlated with tumor purity but negatively correlated with immune score, indicating that BER pathway activity and APE1 expression may be affected by tumor purity in bulk RNA-seq data. In contrast, the four subtypes were not correlated with tumor purity, indicating that the measurement of the relative abundance of each melanoma subtype was less affected by tumor purity (Supplementary Figure S2).

The heterogeneity of TME may have an impact on the association of CD4+T cells and the relative abundance of melanoma subtypes. To comprehensively elucidate this association, we classified patients into three subtypes based on their TME:D subtype featured with the lowest immune infiltration, NEF subtype with high infiltration of neutrophil cells, endothelial cells, and fibroblasts, and NKT subtype characterized with high infiltration of NK cells, T cells, and B cells. As expected, the results from GSVA demonstrated that NKT subtype exhibits high activity of type I and type II IFN response, immune checkpoint, NEF subtype displays elevated activity of angiogenesis, pan-fibroblast signature and D subtype is distinguished from these two subtypes by high tumor cell intrinsic features such as cell cycle and DNA damage repair (Supplementary Figure S3A). The reproducibility of three immune subtypes was maintained by two other independent melanoma cohorts profiled with microarray and bulk RNA-seq, respectively (Supplementary Figure S3B and C).

Consistent with the results from the MGSP cohort, Kaplan-Meier analysis demonstrated that melanoma patients with high CD4+ T cell infiltration exhibited better overall survival (OS) compared to those with low infiltration (log-rank p<0.001) (Figure 3A). As to three immune subtypes, the NKT subtype was significantly associated with prolonged OS, as shown in Figure 3B (log-rank p<0.001, NKT vs. D adjusted p<0.001; NKT vs. NEF adjusted p=0.001). Subsequently, we investigated the associations between melanoma subtypes and 19 immune pathways. Univariate linear regression showed that the relative abundance of BER-subtype melanoma cells was significantly positively associated with CD4+ T cell infiltration only in the NKT subtype (Supplementary Figure S4). As shown in Figure 3C, the positive correlation of relative abundance of BER-subtype melanoma cells with CD4+ T cell infiltration remained significant within the NKT subtype in multivariate linear regression (p<0.001). However, in the NEF subtype, CD4+ T cell infiltration was not significantly related to BER-subtype melanoma, and a negative correlation was observed between CD4+ T cell infiltration and BER-subtype melanoma in the D subtype (Figure 3C).

Figure 3:

CD4+ T cells and immune subtypes in melanoma patients from the TCGA SKCM cohort. (A) Kaplan-Meier plot showing different OS between patients with high or low CD4+ T cell infiltration. (B) Kaplan-Meier plot showing different OS among patients with D, NEF, or NKT subtype. (C) heatmap showing the regression coefficients and their significance from multivariate linear regression in each immune subtype. The asterisk ★, ★★, ★★★ labeled in the squares indicate significance level of 0.05–0.01, 0.001–0.01, and less than 0.001, respectively. The groups of high or low CD4+ T cell infiltration were determined by the median value of the enrichment score of CD4+ T cells from consensus TME.

IHC assay identified the positive correlation between APE1 expression and CD4+ T cell infiltration

To further investigate the association between BER activity and tumor-infiltrating CD4+ T cell abundance, a melanoma tissue microarray comprising 89 melanoma patients was used. The baseline clinicopathological characteristics of the cohort are summarized in Supplementary Table S6. IHC was employed to assess the expression levels of APE1, the rate-limiting enzyme in the BER pathway, and the density of CD4+/CD8+ TILs. Representative IHC images for APE1, CD4, and CD8 cells are shown in Figure 4A. The proportions of categorized CD4+ T as an ordinary variable were significantly different among the three APE1 categories revealed by the Kruskal-Wallis test (Supplementary Table S7). Considering the percentage of CD4 or CD8 positive staining as a continuous variable, the percentage of CD4 positive staining significantly increased with the elevation of APE1 staining density. The boxplots were shown in Figure 4B. As illustrated in Figure 4C, APE1 expression in melanoma cells demonstrated a significant positive correlation with categorized CD4+ T cell infiltration (Pearson’s r=0.58, p<0.001), consistent with the results of BER type melanoma observed in GSE222446 and TCGA SCKM NKT subtype. In contrast, no significant correlation was observed between APE1 expression and CD8+ T cell infiltration.

Figure 4:

APE1 expression levels correlate with CD4+T cell infiltration. (A) representative IHC images of APE1 and TILs infiltration of CD4+, CD8+ cells. (B) boxplots showing the difference in the percentage of CD4+ T cells and CD8+T cells among three APE1 staining density groups. (C) correlation plots of CD4+ T or CD8+ T cells infiltration and APE1 expression. The number values in the colored circles represent the spearman correlation coefficients with significance less than 0.05.