Identification of CircRNA signature associated with tumor immune infiltration to predict therapeutic efficacy of immunotherapy – Nature Communications

CircRNA profiling in two independent melanoma patient cohorts undergoing ICB treatment

To characterize the significant roles that circRNAs may play in differentiating between response and non-response in patients undergoing ICB treatment, we retrieved two independent total RNA sequencing datasets for individuals treated with single-agent anti-PD-1 and/or combined anti-CTLA-4 and anti-PD-1. Cohort 1 consisted of 88 patients treated with anti-PD-1 monotherapy (n = 47, including pre-treatment [PRE] samples, n = 38, and early during therapy [EDT] samples, n = 9) or combined anti-CTLA-4 and anti-PD-1 therapy (patients, n = 41, including PRE, n = 32 and EDT, n = 9; Supplementary Data 1), while Cohort 2 consisted of 69 melanoma patients undergoing anti-PD-1 therapy (nivolumab or pembrolizumab) (Supplementary Data 2). The unmapped reads per patient ranged from 1,646,487 to 28,592,124 reads in Cohort 1 and from 1,044,568 to 7,607,887 reads in Cohort 2 (Supplementary Fig. 1a, b), and these data were used to quantify circRNA expression. To reliably identify circRNAs, we combined four well-established circRNA-detection tools with user-friendly computational algorithms²¹, including CIRI2⁵⁶, find_circ⁵⁷, CircExplorer2⁵⁸, and CircRNA_finder⁵⁹, to quantify back splice-spanning reads (Fig. 1a, b; See Methods). There was no significant correlation between the number of unmapped reads and the number of circRNAs detected by any of these circRNA-detection tools (Supplementary Fig. 1c, d), suggesting that the detectable number of circRNAs is not dependent on the number of unmapped reads. Different circRNA-detection tools identified varying numbers of circRNAs⁶⁰ (Fig. 1a, b), so we kept those circRNAs identified by at least two tools with ≥2 back-splicing reads. We ultimately identified 89,204 total circRNAs from the 88 samples in Cohort 1 (Fig. 1c, Supplementary Data 3), and 43,911 circRNAs from the 69 samples in Cohort 2 (Fig. 1d, Supplementary Data 4). The detectable number of circRNAs ranged from 1440 to 16,653 in each patient from Cohort 1 (Supplementary Fig. 1e), and from 39 to 11,652 in each patient from Cohort 2 (Supplementary Fig. 1f). To minimize potentially spurious events, we only considered circRNAs identified in more than 20% of the total samples in each cohort. A total of 5350 circRNAs were retained in Cohort 1, ranging from 774 to 4525 in each patient (Fig. 1e), while a total of 3654 circRNAs were retained in cohort 2, ranging from 20 to 3013 in each patient (Fig. 1f).

In patient Cohort 1, the 5350 identified circRNAs originated from 2678 host genes, with most of these host genes (1515/2678 = 56.6%) having at least one circRNA, while 16 host genes were associated with more than 10 circRNAs (Fig. 1g). In patient Cohort 2, the 3654 identified circRNAs originated from 2110 host genes, and 1333 host genes were associated with one circRNA, while 5 host genes were associated with more than 10 circRNAs (Fig. 1h). We found that circRNAs detected in this study were ubiquitously located across whole genomic regions (Supplementary Fig. 1g, h), and most circRNAs were back-spliced from exonic regions (Supplementary Fig. 1i, j). We performed overlap analyses among the two cohorts in this study, circAtlas⁶¹, and the MiOncoCirc^62,63 database. A significant overlap of detectable circRNAs was observed between Cohort 1 and Cohort 2, with 90.4% (3293/3644) of the circRNAs in Cohort 2 having been detected in Cohort 1 (Fisher test, p < 2.2e-16; Fig. 1i). We further found that 96.6% of the circRNAs (5168/5350) in Cohort 1 and 97.3% (3544/3654) in Cohort 2 were identified in MiOncoCirc database, suggesting circRNAs in these two cohorts have also been identified in various tumor tissues. 98.5% of the circRNAs (5271/5350) in Cohort 1 and 98.8% (3612/3654) in Cohort 2 were identified in the circAtlas database⁶¹. Furthermore, 95.4% of circRNAs (5106/5350) in Cohort 1 and 96% (3508/3654) in Cohort 2 were identified in both the circAtlas and MiOncoCirc databases, suggesting the conserved expression of these circRNAs in both normal and tumor samples.

Identification of circRNAs associated with immune responses through circRNA-miRNA-mRNA regulatory axes

To better understand the molecular regulation of tumor responses to immunotherapy, we used a linear mixed effects (LME) model to identify circRNAs that were differentially expressed between the response (n = 54) group, defined as partial/complete response (PR/CR) or stable disease (SD) with progression-free survival (PFS) >3 months and the non-response (n = 16) group, defined as progressive disease (PD) or SD with PFS < = 3 months before treatment (PRE), with the expression of these differentially expressed circRNAs then being further examining early during therapy (EDT) (response, n = 13; non-response, n = 5). In the PRE biopsies, we found 227 upregulated and 23 downregulated circRNAs (P < 0.05 and |log₂ (fold change) | ≥ 0.5) in non-responder patients compared to responder patients. In EDT biopsies, we detected 1547 upregulated and 10 downregulated circRNAs after ICB treatment (Fig. 2a, b). More upregulated circRNAs and fewer downregulated circRNAs were found in the non-responder group at both the PRE and EDT time points, suggesting a potential association between the overexpression of circRNAs and resistance to immunotherapy. CircRNAs were relatively highly abundant in the non-responder group, whereas no significant differences in the expression of their host genes were observed (Supplementary Fig. 2), suggesting the independent roles of circRNAs in the response to immunotherapy.

Volcano plot of upregulated (red) and downregulated (blue) circRNAs between non-response and response melanoma samples in pre-treatment (PRE; a) and early during therapy (EDT; b) samples, respectively (P < 0.05, |log2 (fold change)| ≥0.5; See Methods). c The overlap of upregulated circRNAs identified in PRE (a) and EDT (b) samples (81 shared circRNAs with 80 in MiOncoCirc and 74 host genes). d Strategy of ICB associated circRNA-miRNA-mRNA axes (upper panel, see method), and the network of all circRNA-miRNA-mRNA interaction pairs (bottom panel). “+” represent the positive correlation between the expression of circRNAs and mRNAs. e Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment and GO biological process analysis among the mRNAs identified in ICB-associated circRNA-miRNA-mRNA axes (P < 0.05). The statistical analysis was performed by a two-sided linear mixed-effects model (LME) in a, b and two-sided Fisher’s test in e. PRE pre-treatment, EDT early during therapy. Source data are provided as a Source data Fig. 2.

To further explore the functional effects of circRNAs in immunotherapy, we have assumed that circRNA expression at PRE time point may be elevated in order to be elevated at EDT time point, and thus focused on 81 upregulated circRNAs shared between non-responder and responder patients at the PRE and EDT timepoints for subsequent analysis (Supplementary Table 3), which were derived from 74 host genes (Fig. 2c). We focused on the circRNA-miRNA-mRNA interactions in which circRNA expression was positively correlated with mRNA expression to identify potential functional roles for these upregulated circRNAs in immunotherapy responses. We constructed circRNA-miRNA-mRNA regulatory networks mediated by common miRNAs that bind to the upregulated circRNAs and mRNAs and are correlated with the expression of those circRNAs. We used the miRanda algorithm⁶⁴ to predict high-confidence binding sites, and selected the top 30 miRNAs for each of 81 upregulated circRNAs. A total of 773 miRNAs and 2429 interactions were detected for the circRNA-miRNA axis. Next, we identified the potential targets of miRNA-mRNA interactions using the Tarbase databases⁶⁵, which incorporates experimentally supported interactions, and the TargetScan database. Taken together, these analyses enabled us to identify 184,587 predicted circRNA-miRNA-mRNA interactions based on the miRNA target sites shared by circRNAs and mRNAs. We filtered low-confidence associations based on multiple criteria including expression correlations (see Methods), after which 8449 (81 circRNAs − 183 miRNAs − 2046 mRNAs) ICB response-associated interactions were retained for further study (Fig. 2d; Supplementary Data 6).

Pathway analyses based on these 2046 mRNAs suggested their significant enrichment in cancer signaling pathways, including the Hippo, p53, mTOR, and AMPK signaling pathways (Fig. 2e. These mRNA targets were also enriched in several biological processes, including cell cycle checkpoint, cellular response to hypoxia, autophagy, and regulation of Wnt signaling pathways. For example, targeting Wnt/β-Catenin signaling may reverse the resistance to immunotherapy by altering antigen presentation⁶⁶. Our analysis revealed the unexpected upregulation of circRNAs in non-responder patients, suggesting that certain circRNAs may mediate tumor resistance to immunotherapy through the alteration of cancer signaling pathways.

Identification of an ICBcircSig prognostic model to predict immunotherapeutic efficacy

CircRNAs may serve as potential biomarkers in cancer⁶⁷, but it is unclear whether they can serve as biomarkers capable of predicting patient responses to cancer immunotherapy. To identify prognosis-related cirRNAs, we performed univariate Cox regression analyses examining the relationship between patient PFS and the expression levels of 227 upregulated circRNAs in PRE biopsies from Cohort 1. We found that high levels of expression for 25 circRNAs were significantly associated with poorer PFS (log-rank test, FDR < 0.05, and Cox FDR < 0.05). To identify the optimal circRNAs for use as prognostic biomarkers, we employed a LASSO Cox regression model analysis⁶⁸ using the expression profiles for these 25 circRNAs and associated clinical information, and ultimately selected four circRNAs with non‐zero regression coefficients (Supplementary Fig. 3a, b). We considered these four circRNAs as variables in a subsequent multivariate Cox regression analysis and found both circTMTC3 and circFAM117B to be significant predictors of patient survival (Supplementary Fig. 3c). Specifically, the expression of these two circRNAs was associated with worse PFS (Fig. 3a, b), and patients that exhibited CR/PR or SD following anti-PD-1 therapy exhibited significantly lower expression of circTMTC3 and circFAM117B as compared to patients with PD (Fig. 3c, d). We further collected 11 and 4 pre-treatment tumor biopsies from patients who did or did not respond to anti-PD1 treatment (Supplementary Table 1). Expression levels of circTMTC3 and circFAM117B were detected in these samples by qRT-PCR. We found that both circRNAs were significantly upregulated in the non-responder group relative to the responder group (Fig. 3e, g). Sanger sequencing further confirmed that these PCR products spanned the circular junction for both circTMTC3 and circFAM117B (Fig. 3f, h).

Kaplan–Meier survival curves show expression of circTMTC3 (a) and circFAM117B (b) in the ICBcircSig associated with PFS. The expression of circTMTC3 (c) and circFAM117B (d) expression were significantly higher in PD groups (n = 16) than CR/PR or SD groups (n = 54). e–h The qRT-PCR and Sanger sequencing validate expression and junction of circRNAs. The expression of circTMTC3 (e) and circFAM117B (g) were detected by qRT-PCR in pre-treatment tumor biopsies from patients show response (n = 11) and non-response (n = 4) to anti-PD-1 therapy. The schematic diagram show generation of circTMTC3 (f, left panel) and circFAM117B (h, left panel). Sanger sequencing confirmed that PCR products spanned the circular junction of predicted circTMTC3 (f, right panel) and circFAM117B (h, right panel) in the patient with non-response to anti-PD-1 therapy. The magenta dash line indicates the circular junction site. i ROC curves quantifying ICBcircSig score response prediction AUC in cohort 1. j Kaplan–Meier survival curves of PFS between high- and low-risk patients stratified by ICBcircSig score using the optimal cutoff. k Time-dependent ROC curve at 12 and 24-months of PFS for the ICBcircSig score. l Forest plot for the HRs of multivariate Cox model of the ICBcircSig score and clinicopathological variables. Black vertical lines indicate the 95% confidence interval (CI). m Boxplot of ICBcircSig score among NR (n = 16) and R (n = 54) groups. n Kaplan–Meier survival curves of OS between high- and low-risk patients stratified by ICBcircSig score. A log-rank test was used in a, b, j, and n. A two-sided Wilcoxon rank-sum test was used in c–e, g, and m. The box in c–e, g, and m showed the median ± 1 quartile, with the whiskers extending from the hinge to the smallest or largest value within 1.5× IQR from the box boundaries. PFS progressive free survival, OS overall survival, CR/PR complete response/partial response, SD stable disease, PD progressive disease, ROC receiver operating characteristic curve, HR hazard ratio, R responder, NR non-responder. Source data are provided as a Source data Fig. 3a–e, g, i–k, m, n.

We further constructed an ICB-related circRNA signature (ICBcircSig) score by weighting the expression values of circRNAs in ICBcircSig based on their established multivariate Cox regression coefficient values. When assessing the clinical relevance of these scores, we found that the ICBcircSig score was able to distinguish between responder and non-responder patients with an AUC of 0.8 (Fig. 3i). We then determined that patients with a high ICBcircSig score had worse PFS compared with those patients with a low ICBcircSig score (log-rank test, p < 0.001, Fig. 3j). The 12 and 24-month progression rates in the high ICBcircSig score group were 100% and 100%, respectively, which were significantly higher than the rates (27% and 30%, respectively) in the low ICBcircSig score group. The respective AUCs of the time-dependent ROC curves for this ICBcircSig score were 0.76 and 0.75 for 12- and 24-month PFS (Fig. 3k). The predictive ability of the ICBcircSig score was validated through the random sampling of 90% of these samples 10 times, yielding mean time-dependent ROC curve AUC values for the ICBcircSig score of 0.757 and 0.754 for 12- and 24-month PFS (Supplementary Fig. 3d, e), respectively.

To further examine whether the ICBcircSig scores could serve as an independent prognostic factor in melanoma patients, we performed multivariate Cox regression analyses to adjust for the potential confounding effects of other conventional clinical factors, including age, gender, ICB treatment, CD274 (PD-L1), and PDCD1. The ICBcircSig score (hazard ratio [HR] = 2.975, 95% confidence interval [CI] 1.723–5.138, p < 0.001) remained an independent prognostic risk factor for PFS even after adjusting for these confounding factors (Fig. 3l). We further investigated the association between the ICBcircSig score and patient response to ICB treatment and demonstrated that responders exhibited significantly lower ICBcircSig scores as compared to non-responders (R vs NR, p = 1.8 × 10⁻⁴; Fig. 3m). In terms of overall survival, a higher ICBcircSig score was also associated with a shorter survival duration (log-rank test, p = 0.028, Fig. 3n).

Validation of the performance of ICBcircSig prognostic model in an independent patient cohort

We further assessed the performance of the ICBcircSig score in an independent patient cohort. CircTMTC3 and circFAM117B were associated with worse PFS and tended to be enriched in the patients that exhibited a PD response to ICB treatment in Cohort 2 (Supplementary Fig. 4a–d). We further calculated the ICBcircSig scores for each sample, and found ICBcircSig scores to be higher in the non-responder group (Fig. 4a), consistent with our observation in Cohort 1. ICBcircSig scores were also able to distinguish between responders and non-responders (AUC = 0.66; Fig. 4b). Survival analyses revealed that patients with high ICBcircSig scores exhibited worse PFS as compared to patients with low ICBcircSig scores (Fig. 4c), and the respective AUCs of the time-dependent ROC curves for the ICBcircSig score were 0.69 and 0.65 for 12- and 24-month PFS (Fig. 4d). The predictive ability of the ICBcircSig score was further validated via the random sampling of 90% of these samples 10 times, yielding respective mean AUC values for ICBcircSig scores of 0.678 and 0.645 for 12- and 24-month PFS (Supplementary Fig. 4e, f). To assess whether the ICBcircSig score is an independent predictor for PFS, we considered stage, gender, CD274 expression, PDCD1 expression, TMB, lactate dehydrogenase (LDH), and ICBcircSig scores as variables in a multivariate Cox analysis, and we ultimately found that ICBcircSig score (HR = 1.32, 95% CI: 1.0721–1.630, p = 0.009) was an independent predictor significantly associated with worse PFS (Fig. 4e). We further found that higher ICBcircSig scores were associated with worse overall survival (Fig. 4f), and these ICBcircSig score model exhibited strong prognostic performance, with respective AUC values for the ICBcircSig score of 0.71 and 0.67 for 12- and 24-month OS (Fig. 4g). ICBcircSig score remained an independent predictor of melanoma patient OS even after adjusting for confounding factors (Fig. 4h).

a Boxplot of ICBcircSig score between response (n = 38) and non-response (n = 30) groups. b Receiver Operating Characteristic (ROC) curves quantifying ICBcircSig score prediction AUC. c Kaplan–Meier survival curves of PFS between high- and low-ICBcircSig score patients stratified by the optimal cutoff in cohort 2. d Time-dependent ROC curve at 12 and 24-months of PFS for the ICBcircSig score. e Forest plot for the HRs of multivariate Cox model of the ICBcircSig score and clinicopathological variables. f Kaplan–Meier survival curves of OS between high- and low-risk patients stratified by ICBcircSig score using the optimal cutoff in validation data. g Time-dependent ROC curve at 12 and 24-months of OS for the ICBcircSig score. h Forest plot for the HRs of multivariate Cox model of the ICBcircSig score for OS and clinicopathological variables. I, j Boxplot of expression of circTMTC3/circFAM117B (i) ICBcircSig score (j) distribution between response (n = 14) and non-response (n = 7) groups for in-house cohort 3. A log-rank test was used in c, f and i. k ROC curves quantifying ICBcircSig score response prediction AUC in in-house cohort 3. l Kaplan–Meier survival curves of PFS between high- and low-risk patients stratified by ICBcircSig score using the optimal cutoff in in-house cohort 3. m Time-dependent ROC curve at 12 and 24-months of PFS for the ICBcircSig score in in-house cohort 3. A two-sided Wilcoxon rank-sum test was used in a, i, and j. The boxes in a, i and j indicate the median ± 1 quartile, with the whiskers extending from the hinge to the smallest or largest value within 1.5× IQR from the box boundaries. Black vertical lines in e and h indicate the 95% confidence interval (CI). PFS progressive free survival, OS overall survival, ROC receiver operating characteristic curve, HR hazard ratio, AUC Area Under the ROC Curve. Source data are provided as a Source data Fig. 4a–d, f, g, i–l.

In addition, we validated the ICBcircSig score using our in-house cohort (Cohort 3) consisting of 24 patients with melanoma undergoing anti-PD-1 treatment (Supplementary Table 2, Supplementary Data 7). Consistently, circTMTC3, circFAM117B, and ICBcircSig score values were significantly enriched in non-responders (Fig. 4I, j). In this cohort, the ICBcircSig score exhibited an AUC of 0.85 when predicting patient responses to ICB treatment (Fig. 4k). A higher ICBcircSig score was associated with worse PFS (Fig. 4l), and the AUCs of time-dependent ROC curves for the ICBcircSig score at 12 and 18-months were 0.82 and 0.76, respectively (Fig. 4m). Taken together, our results provide proof-of-concept evidence that this ICBcircSig score may accurately predict patient responses to immunotherapy.

The ICBcircSig score outperforms other transcriptome-based signatures

Previous studies have explored multiple transcriptome-based signatures to predict patient response to ICB treatment^{35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51}. We next compared the performance of the ICBcircSig score with 20 previously reported signatures (Supplementary Table 3) to assess the predictive ability of these tools in our patient cohorts. The ROC classification curves for non-responders and responders based on the ICBcircSig score yielded AUC values of 0.80, 0.66, and 0.85 for Cohort 1, Cohort 2, and our in-house cohort (Fig. 5a), respectively. The mean AUCs of time-dependent ROC curves for this ICBcircSig score when predicting the PFS of patients in these three respective cohorts were 0.771, 0.676, and 0.825, respectively, with these values being higher than those for other extant predictors (Fig. 5b). We applied a univariate Cox regression model to assess the relationship between each signature and PFS, and found that 12, 3, and 1 of these 21 signatures were significantly associated with the PFS of patients in Cohort 1, Cohort 2, and our in-house cohort, respectively, with ICBcircSig score achieving the highest HR in all patient cohorts (HR = 4.01, 95% CI: 2.38–6.75, p < 0.05 in Cohort 1, HR = 1.35, 95% CI: 1.11–1.66, p < 0.05 in Cohort 2, HR = 1.35, 95% CI: 1.22–2.22, p < 0.05 in the in-house cohort; Fig. 5c). We further compared the ability of transcriptomic biomarkers to differentiate between ICB responders and non-responders by performing the Wilcoxon Rank Sum test (Fig. 5d, left panel) and univariate logistic regression analyses (Fig. 5d, right panel). The results indicated that ICBcircSig scores were both significantly able to distinguish patients with different responses to ICB treatment and were superior to other predictive tools. Finally, we evaluated the association between each signature and overall survival based on AUC values, and found that ICBcircSig scores also exhibited the highest mean AUC of 0.66 (Supplementary Fig. 5). Taken together, these findings demonstrate that this ICBcircSig score outperforms existing transcriptome-based signatures when predicting responses to ICB treatment.

a Performance of ICBcircSig score and other published signatures based on ROC curves quantifying ICBcircSig score response prediction AUC in each cohort. b AUC of time-dependent ROC of 6, 12, and 18-months for each cohort. c HR for FPS in each cohort. The statistical difference is performed by cox regression model. Red dashed line indicated p value 0.05. d P value of two-sided Wilcoxon rank-sum test (left panel) and logistic regression test (right panel) between response versus non-response, red dashed line indicated p value 0.05. Blue point means 6 months, red point means 12 months and green point means 18-months in b (n = 3). The boxes in b indicate the median ± 1 quartile, with the whiskers extending from the hinge to the smallest or largest value within 1.5× IQR from the box boundaries. ROC receiver operating characteristic curve, HR hazard ratio, AUC Area Under the ROC Curve. Source data are provided as a Source data Fig. 5.

The association between ICBcircSig score and cancer hallmarks and immune features

To explore the functional implication of the ICBcircSig score in immunotherapy, we performed an enrichment analysis of 50 hallmark genes^69,70 between high ICBcircSig score and low ICBcircSig score groups based on a gene set enrichment analysis (GSEA) approach. The low ICBcircSig score group was enriched for genes associated with the immune response, including the IFN-α response (normalized enrichment score [NES] = −2.58, p = 0.0024), IFN-γ response (NES = −2.72, p = 0.0027), inflammatory response (NES = −2.03, p = 0.0027), and IL6/JAK/STAT3 signaling pathway (NES = −1.88, p = 0.0023) in patient Cohort 1 (Fig. 6a). Consistently, immune response-related cancer hallmark pathways such as the IFN-α response (NES = −2.13, p = 0.029) and IL6/JAK/STAT3 signaling pathway (NES = −1.95, p = 0.023) were also enriched in the low ICBcircSig score group in Cohort 2 (Supplementary Fig. 6a), suggesting that patients with low ICBcircSig scores may exhibit a more activated immune microenvironment. To further examine whether the ICBcircSig scores were associated with immune infiltration, we conducted GSEA and single sample GSEA (ssGSEA) analyses of the infiltration of 22 immune subpopulations⁷¹ in tumors from patients with low or high ICBcircSig scores. We observed significant differences in immune infiltration when comparing these two ICBcircSig score-based groups (Fig. 6b and Supplementary Fig. 6b) in both tested cohorts. Specifically, 17 and 18 out of the 22 analyzed immune subpopulations exhibited higher levels of predicted infiltration in the low ICBcircSig score group in Cohort 1 and Cohort 2, respectively. These results suggest that patients with low ICBcircSig scores exhibit characteristically high levels of tumor immune infiltration that may explain their higher response rates to immunotherapy. In addition, we assessed the association between ICBcircSig score and cytolytic activity (CYT), a proxy used to reflect the ability of T cells to kill cancer cells⁵⁰, and GEP, which corresponds to a T cell-inflamed environment⁴⁹, and we observed that the high ICBcircSig score group had a lower CYT score/GEP score than the low ICBcircSig score group in both tested cohorts (Fig. 6c, d and Supplementary Fig. 6c, d), suggesting their less robust cytolytic activity and impaired ability to kill tumor cells.

a Volcano plots for the enrichment of hallmarks for cohort 1 samples with high and lowICBcircSig score based on the Normalized Enrichment Score (NES) from the GSEA. b The enrichment of immune cell for samples with high and low ICBcircSig score based on the Normalized Enrichment Score (NES) from the GSEA. c Cytotoxic T cell score in the high- (n = 11) and low- (n = 59) risk groups stratified by the ICBcircSig score. d GEP score in the high- and low-risk groups stratified by the ICBcircSig score. e–g The expression of circTMTC3 (e), miR-142-5p (f), and PD-L1 (g) in SK-MEL-28 melanoma cell line with circTMTC3 overexpression (oe) or control group (n = 3 in each group). h–j The expression of circFAM117B (h), miR-142-5p (i), and PD-L1 (j) in SK-MEL-28 melanoma cell line with circFAM117B overexpression (oe) or control group. k–m SK-MEL-28 melanoma cells with or without overexpression of circTMTC3 or circFAM117B co-cultured with activated T cells for 48 h were subjected to crystal violet staining (k). SK-MEL-28-to-T cell ratio, 1:3. l, m Statistical analysis comparing the ability of T cell killing in (k), n = 3 in each group. Two-sided Wilcoxon rank-sum test was used in c, d. The boxes in c, d indicate the median ± 1 quartile, with the whiskers extending from the hinge to the smallest or largest value within 1.5× IQR from the box boundaries. The unpaired two-sided Student T-test was used in e–j, one-way ANOVA and Tukey’s multiple comparison test was used in l, m. The error bars in e–j and l, m indicate the mean ± s.d. *p < 0.05, **p < 0.01 and ***p < 0.001. Source data are provided as a Source data Fig. 6a–d.

To validate the functional role of circRNA signatures is based on the circRNA-miRNA-mRNA regulatory axe. We generated the circTMTC3 overexpression (OE) cell lines and circFAM117B OE cell lines in SK-MEL-28 (Fig. 6e, h). According to miranda’s prediction27, hsa-miR-142-5p can interact with circTMTC3 and circFAM117B, we observed hsa-miR-142-5p significantly downregulated in OE cell lines of circTMTC3 or circFAM117B (Fig. 6f, i). Then, we found the expression of PD-L1 is significantly upregulated in both OE cell lines by RT-PCR analysis (Fig. 6g, j). These suggest that ICBcircSig signature circTMTC3 and circFAM117B can regulate the hsa-miR-142-5p/PD-L1 pathway in melanoma cell line. To further investigate the functional role of circTMTC3 and circFAM117B in immunosuppression, we performed an in vitro T cell cytotoxicity-mediated tumor killing assay based on SK-MEL-28 melanoma cells overexpressing circTMTC3 or circFAM117B. The overexpression of circTMTC3 or circFAM117B significantly reduced the CD8⁺ T cell cytotoxicity and the ability of these T cells to eliminate tumor cells (Fig. 6k–m). Together, we thus identified a strong correlation between ICBcircSig score and a series of immune signatures that reflect the complex TME, highlighting the prognostic value of this ICBcircSig score and the roles played by circTMTC3 and circFAM117B in immunosuppression.