A total of 25 SKCM samples with scRNA-seq data, previously published by Li et al. (19), was downloaded for screening Tregs-associated marker genes. Data from The Cancer Genome Atlas (TCGA) bulk-seq and clinical data for 458 patients with SKCM were obtained to identify survival-related genes and create a signature. To confirm the predictive power of the developed signature, separate groups of patients were obtained from the GSE65904 datasets (n=209). Data from single-cell RNA sequencing of melanoma patients treated with carfilzomib or ICI were gathered from three separate datasets (GSE161801, GSE120575, and GSE189125) to assess the effectiveness of Treg.Sig in predicting the response to immunotherapy. The research made use of publicly accessible datasets that had received ethics approval from the original studies. The study was conducted according to the Helsinki declaration.

Analysis of single-cell RNA sequencing data was conducted with the R software tool “Seurat” (version 4.0.5) (20). Merged individual data and excluded low-quality cells based on gene detection criteria (nFeature_RNA >= 200 & nCount_RNA >= 500). Then, the “harmony” de-batching method was performed on 30 cases of SKCM scRNA-seq data. Following that, the scRNA-seq data that had been standardized underwent UMAP analysis for dimension reduction and were grouped using the “Seurat” tool in R software, with the parameters scale_factor=10000, nfeatures=3000, npcs=50, dims=30, and resolution = 0.2. Cell types were annotated based on markers using the SingleR package (version 1.0.0), and the expression patterns of each marker gene within clusters were analyzed with the “DotPlot” function in Seurat. Subsequently, potential marker genes for the Tregs cell cluster were pinpointed utilizing the FindAllMarkers feature.

To analyze bulk data, we used weighted co-expression network analysis (WGCNA) on the differentially expressed genes (DEGs) to find modules that are most relevant to Tregs. Initially, the correlations between genes and genes were computed in order to create a similarity matrix. Next, the optimal soft-thresholding power β was chosen using the pickSoftThreshold function within the “WGCNA” software package. By applying soft-thresholding, the similarity matrix was transformed to create a scale-free topology. Third, the adjacency matrix was transformed into a topological overlap matrix by similarity, and the corresponding dissimilarity was also calculated.

Finally, co-expression gene modules were determined by utilizing the R software tool “Dynamic Tree Cut” with a deepSplit value of 2, a minimum module size of 30, and a maximum block size of 20,000. Modules that were very alike were combined if the module eigengene height in the clustering was less than 0.25. As a result, the modules were connected to clinical characteristics in order to pinpoint the module with genes that were most pertinent to the Tregs.

A novel machine learning-guided computational framework based on transcriptome data from scRNA-seq and bulk-seq was developed for identifying Treg.Sig as follows ( Figure 1 ):

For performance evaluation:

To be clinically useful, a nomogram was created using the multivariate Cox regression coefficients for age, sex, cancer stage, and Treg.Sig. The nomogram’s predictive performance was validated by calculating the AUCs, examining the consistency between predicted 1-, 3-, and 5-year overall survival (OS) probabilities and actual observations with calibration curves, and measuring the clinical utility of the nomogram with decision curve analysis (DCA).

To further estimate the prognostic performance of Treg.Sig, we compared Treg.Sig with other melanoma-specific signatures. 51 signatures were included ( Supplementary Table S2 ), consisting of lncRNA, miRNA, and mRNA signatures, for comparison of the 1-year, 3-year, and 5-year AUC in external datasets GSE65904.

64 types of immune infiltrating cells were calculated. The XCELL algorithm was used to quantify the infiltrating level of immune cells. Furthermore, the potential immunotherapy value of Treg.Sig was assessed using the immune proportion score (IPS), which measures the percentage of PD-L1 expression on immune cells, and the TIDE score, which indicates the expression signature of two main mechanisms of tumor immune evasion: T cell dysfunction and T cell exclusion. These two biomarkers could forecast the response to immune checkpoint inhibitors in individuals with cancer (21). The IPS of each individual with SKCM were acquired from The Cancer Immunome Atlas (TCIA) website (https://tcia.at/home), while the TIDE score was retrieved from the TIDE web platform (http://tide.dfci.harvard.edu).

Kyoko encyclopedia of genes and genomes (KEGG) and Gene ontology (GO) were used to analyze the immune enrichment of Treg.Sig marker gene with the help of R software packages: “ClusterProfiler” (version 4.0.5), “org.Hs.eg.db” (version 3.13.0), “ggplot2” (version 3.3.5), and “enrichplot” (version 1.12.3). The R software package “clusterProfiler” was utilized to conduct gene set enrichment analysis (GSEA) in order to explore the connection between Treg.Sig and immune response.

Currently, advanced SKCM patients are primarily treated with chemotherapy and ICI therapy. Our study selected four commonly used chemotherapeutic agents of SKCM, including AP_24534, Pazopanib, Paclitaxel, and Tivozanib. The half-maximal inhibitory concentration (IC50) chemotherapy drugs was estimated by utilizing the R software package “pRRophetic algorithm” and “pRRophetic” with data from the Cancer Genome Project (CGP) cell lines (22), in order to forecast the response to chemotherapy in groups categorized as high- and low-Treg.Sig.

Furthermore, the ICI response was predicted using SKCM patients from two GEO datasets (GSE120575 and GSE189125), with both scRNA-seq data and ICI treatment results. Meanwhile, the GSE161801 cohort comprised scRNA-seq data from pre-treatment melanomas receiving carfilzomib therapy. These data was also used to determine the potential value of Treg.Sig in predicting therapeutic responsiveness. The AddModuleScore function in the Seurat package was utilized to compute Treg.Sig scores in scRNA-seq data.

Candidate differentially expressed genes (DEGs) were submitted to the STRING online database for analysis of protein-protein interactions (PPI), and the resulting network was visualized and hub genes were identified using Cytoscape software. Data from the HPA database (https://www.proteinatlas.org/) was used to analyze the protein expression of the hub gene STAT1 in SKCM tissues. To investigate the potential immunotherapy response of STAT1, six immune cell types were evaluated: activated dendritic cells (aDC), B cells, CD8 T cells, cytotoxic cells, macrophage cell, and T helper cells; along with three immune checkpoint molecules: PD-1, PD-L1, and PD-L2.

Moreover, the frequency of STAT1 mutations in melanoma was determined by analyzing the TCGA PanCancer Atlas studies (n=448) on cBioPortal (https://www.cbioportal.org/). In addition, in order to investigate the connection between STAT1 mutation and response to ICI treatment, data from clinical and tumor whole exome sequencing (WES) of melanoma patients from two cohorts treated with ICIs, including 38 samples from UCLA (Cell 2016) (23) and 64 samples from MSKCC (NEJM 2014) (24), were obtained for analysis. Every individual in the cohort treated with ICIs has received treatment with antibodies that target PD-1 and CTLA-4. In our research, melanoma patients with any nonsynonymous mutations in STAT1, such as missense, translation start site, nonstop, splice site, frameshift, and nonsense variants, were categorized under STAT1 mutations (STAT1 Mut). Conversely, patients without any STAT1 mutations were classified as having the STAT1 wild-type (STAT1 Wt).

The Servicebio Technology Institute provided the B16-F10 melanoma cell line. These cells were cultured in RPMI-1640 supplemented with 10% FBS, 1% penicillin-streptomycin. All cells were grown in a 5% CO2 humid atmosphere at 37°C. For N1 neutrophils polarization, mouse primary neutrophils were treated with 1ug/ml LPS and 20 ng/ml IFN-gamma. For N2 neutrophils polarization, mouse primary neutrophils were isolated from femur and tibia, then purified by using MACS-based neutrophil isolation kit (Milteyni Biotec) following the manufacturer’s instructions, and then co-cultured with B16-F10 melanoma cells. For STAT1 inhibition, mouse primary neutrophils were cultured with 10uM Fludarabine (Flu, STAT1 inhibitor) or DMSO (vehicle).

To implant subcutaneous tumors, B16-F10 melanoma cells were injected into the sides of 8- to 10-week-old mice using a saline solution containing 5×105 cells in 150µl. The size of the tumor was assessed weekly using a digital caliper and determined using the formula: volume = ½ (length × width × width) in millimeters cubed. Mice were sacrificed at end points of the study to analyze tumor weight, following the method outlined in a previous study (25). Mice received 200 µg of anti-PD-1 antibodies or IgG control per mouse through intraperitoneal injection weekly for specified weeks starting on day 7 post subcutaneous implantation of tumor cells. Some mice were injected intraperitoneally with Fludarabine (100 mg/kg per day) to inhibit the phosphorylation of STAT1. To perform adoptive transfer of N2 neutrophils, the N2 phenotype was acquired as outlined in the “Cell culture and treatment” section. The recipient mice were then given intra-tumor injections of N2 neutrophils (5 × 10⁶) that had been suspended in ice-cold PBS at a volume of 200 µl. Mice were sacrificed on the 28th day after implantation, and tumors were harvested for further experiments.

To assess cell proliferation, the Cell Counting Kit-8 (CCK-8; Vazyme, Nanjing, China) method was utilized. B16-F10 melanoma cells were plated and then co-cultured with N2 neutrophils. Following this, each well received 10ul of CCK-8 solution and was incubated for 2 hours at 37°C, shielded from light. The viability of the cells was determined by measuring the absorbance at 450 nm with an enzyme-labeled meter (A33978, Thermo, USA), monitoring the changes at 24, 48, 72, and 96-hour intervals.

The transwell assay was conducted using a two-chamber invasion assay. In the upper chamber, 5×104 B16-F10 melanoma cells were seeded. Additionally, a bottom chamber with 5×104 N2 neutrophils was available. Following incubation, the cells beneath the filter were treated with 4% PFA, dyed with 0.1% crystal violet solution (Beyotime, Shanghai, China), and observed under a light microscope for counting.

RNA extraction utilized TRIzol reagent (Invitrogen, Carlsbad, California), and the NanoDrop 2000 quantified the concentration of total RNA. cDNA synthesis followed, employing the ReverTra Ace qPCR RT Master Mix with gDNA remover (FSQ-301; Toyobo). The SYBR Green Kit (Vazyme, Nanjing, China) facilitated the quantitative PCR, using GAPDH as the internal standard. Primer designs were courtesy of Tsingke Biotech (Beijing, China), with sequences detailed in Supplementary Table S3 .

Cell lysis was performed using RIPA buffer on ice, with protein concentrations determined via the BCA assay (Beyotime, Shanghai, China). Proteins from lysates underwent separation on 10% SDS-PAGE, followed by PVDF membrane transfer. Blocking occurred in 5% non-fat milk for an hour at ambient temperature before overnight primary antibody incubation at 4°C. Following triple washes with TBST, membranes were exposed to secondary antibodies for an hour at room temperature. Detection utilized an ECL kit.

Tissue sections that were fixed in formalin and embedded in paraffin were processed by removing wax, rehydrating, and retrieving antigens. In order to inhibit natural peroxidase activity, the sections were exposed to 3% H2O2 and then blocked with 3% BSA for a duration of 1 hour. Tissue samples were stained with anti-mouse antibodies targeting iNOS, Arg-1, Ly6G, and p-STAT1 overnight at 4°C, followed by secondary antibodies for one hour at ambient temperature, and DAPI to tag nuclei. An inverted fluorescence microscope captured the staining results, measured by ImageJ software.

Cell suspensions from cancer samples (n=3) were washed with PBS and permeabilized using the eBioscience™ Foxp3/Transcription Factor Staining Buffer Set for flow cytometry analysis. Proceed with the steps outlined in the manufacturer’s guidelines (catalog. 00-5523-00). Antibodies that had been conjugated were introduced and left to incubate for one hour at ambient temperature. The antibodies used included Alexa Fluor™ 700-CD11b (eBioscience, 56-0112-80, 1:100), FITC-Ly-6g (eBioscience, 11-5931-82, 1:400), PE-iNOS (eBioscience, 12-5920-80, 1:300), and APC-Arginase (eBioscience, 17-3697-80, 1:100). A flow cytometry recorded cell fluorescence, with FlowJo software analysis and appropriate control antibodies ensuring data accuracy.

Pearson correlation analysis assessed continuous variables, while group differences were evaluated using two-tailed t-tests or one-way ANOVA. Cox hazards regression model identified independent prognostic factors, with significance set at p<0.05. Data analysis and visualization were conducted with R software version 4.1.0 (http://www.R-project.org).

Following the preprocessing of scRNA-seq data using strict quality control measures as indicated ( Supplementary Figure S1 ), 64,071 infiltrating immune cells of high quality were extracted from 30 individuals diagnosed with melanoma, encompassing patients at different stages of the disease and with varied treatment backgrounds. Subsequently, UMAP was performed using Seurat to reduce the dimensionality and identified 12 cell clusters ( Supplementary Figure S2A ), which was subsequently identified as specific cell types through annotation using the SingleR package in R. T-cells, myeloid cells, B-cells, fibroblasts, and NK cells were identified as the main cell types in the study ( Supplementary Figure S2B ). Furthermore, the accuracy of the automated Treg cell cluster annotation was further confirmed by manually validating the Tregs cluster using Treg cell markers (FOXP3, IL2RA, TIGIT, IKZF2, BATF, TNFRSF4, and TNFRSF9) (26) ( Supplementary Figures S2C–F ). Then, 1041 genes related to Tregs were obtained using the Findmarkers function in Seurat, and were considered as candidate genes for further analysis.

WGCNA is a useful bioinformatics technique for analyzing gene correlation patterns in microarray or RNA sequencing data. It detects modules containing genes that are closely related, establishes connections between modules and external sample characteristics, and pinpoints potential biomarkers or targets for treatment. The soft-thresholding power value was a crucial parameter that impacted the autonomy and average connectivity level of co-expression modules. The optimal soft-thresholding exponent was determined based on a scale-free topology fit index of 9 ( Supplementary Figure S3A ). A power value of 6 was selected, which clustered all genes into 24 modules. The hierarchical cluster tree displayed by Supplementary Figure S3B illustrated 24 modules of co-expressed genes analyzed using WGCNA, where genes with similar expression profiles were clustered together. Genes within the identical color block exhibit comparable expression patterns and could potentially share functional relationships. Analyzing the clinical data revealed a strong correlation between the blue module and Tregs, with a Pearson correlation coefficient of 0.47 and a p-value of 6e-27 ( Supplementary Figure S3C ). Furthermore, a strong relationship was observed between gene co-expression in the blue module and two clinical features: OS state (Pearson Correlation = -0.18, P=2e-04) and OS times (Pearson Correlation = 0.13, P=0.007) ( Supplementary Figure S3C ). This indicates that, among the 24 modules, genes in the blue module, consisting of 398 genes, are highly associated with Tregs and tumor recurrence. Further analysis was conducted on the relationship between 398 module membership and gene significance for Tregs, showing a Pearson correlation coefficient of 0.4 (P=1e-16) ( Supplementary Figure S3D ). In conclusion, 398 candidate genes of Tregs were obtained through the bulk sequence.

Finally, 130 intersected candidate genes were identified from scRNA-seq and bulk-seq ( Supplementary Figure S4A ), and were further used to screen out the most valuable Treg marker genes.

Univariate Cox regression identified the prognostic value of 130 potential genes for feature selection. A total of 114 Tregs genes were found to have a significant correlation with OS (P<0.05) in Supplementary Table S1 . Next, the most prognostic genes were further screened out using six feature selection algorithms, all of which were widely used in the machine learning literature. Ten hub genes, including “GBP2”, “RAP1A”, “BTG1”, “IL27RA”, “STAT1”, “TXNDC11”, “GBP4”, “ARID5A”, “CLEC2B”, and “SEL1L3”, were identified for the construction of the prognostic model ( Supplementary Figure S4B ).

Before modeling, we measured the reliability of above 10 hub genes by GO and KEGG enrichment analyses. As illustrated in Supplementary Figure S4C , these Tregs marker genes were associated with immune characteristics, including positive regulation of T-helper 1 type immune response, regulation of T-helper 17 type immune response, and IL-17 production. Supplementary Figure S4D displays the top 10 KEGG pathways that are enriched, primarily related to immune functions, including the NOD-like receptor signaling pathway, Th17 cell differentiation, JAK-STAT signaling pathway, chemokine signaling pathway, Kaposi sarcoma-associated herpesvirus infection, and inflammatory bowel disease. Furthermore, unsupervised machine learning (consensus clustering) indicated that for cluster k values of 2 or 3, analyses of OS rates showed different survival outcomes among the subtypes ( Supplementary Figures S4E, F ). These enrichment terms and unsupervised machine learning confirmed that a combination of machine learning and scRNA- and bulk- sequence data was reliable for Tregs marker screening.

Conducting a multivariate Cox analysis yielded the coefficients associated with the development of Treg.Sig. The detailed formula was as follows: Treg.Sig = (-0.0119 ∗expression value of GBP2) + (-0.0839 ∗expression value of RAP1A) + (-0.0690 ∗expression value of BTG1) + (-0.0829 ∗expression value of IL27RA) + (-0.0294 ∗expression value of STAT1) + (-0.1577 ∗expression value of TXNDC11) + (-0.1128 ∗expression value of GBP4) + (-0.2128 ∗expression value of ARID5A) + (-0.1169 ∗expression value of CLEC2B) + (-0.0492 ∗expression value of SEL1L3).

Then, the Treg.Sig score for each patient was calculated using the above formula. Patients were categorized into high-risk (n=229) and low-risk groups (n=229) based on the median cutoff point of 0.03. Supplementary Figure S5A displays the distribution of risk scores and survival status for each patient. The heatmap exhibited detailed expression level of the enrolled ten genes ( Supplementary Figure S5B ). In the TCGA-SKCM dataset, Kaplan-Meier survival analysis indicated that patients at high risk had notably poorer OS (P<0.001) ( Supplementary Figure S5C ). The prognostic accuracy of Treg.Sig was confirmed by time-sensitive ROC curves showing AUC values of 0.777, 0.762, and 0.721 for 1, 3, and 5 years of OS, respectively (by DeLong’s test; Supplementary Figure S5C ). The performance of Treg.Sig was externally assessed in GSE65904 cohort to validate its robustness. Kaplan-Meier analysis showed that patients with a low Treg.Sig score had notably improved OS (P=0.0032) ( Supplementary Figure S5D ). The validation cohort also demonstrated strong performance in the ROC curves of the risk score ( Supplementary Figure S5D ). The Treg.Sig was identified as a significant factor in the multivariate Cox regression model, with a hazard ratio of 0.421 (95%CI 0.315–0.526, P<0.001), independent of age, gender, and clinical stage ( Supplementary Figure S5E ). Furthermore, Treg.Sig worked better than these factors in comparative ROC curves ( Supplementary Figure S5F ).

To assist clinicians in predicting the OS of SKCM patients, a nomogram was developed that integrated Treg.Sig with age, gender, and tumor stage for practical use. The new nomogram evaluated the chances of survival at 1, 3, and 5 years ( Supplementary Figure S6A ). The ROC curves, calibration, and DCA plots were drawn to evaluate the reliability of the nomogram. In the internal validation cohort, the AUCs in the ROC curve for predicting OS at 1, 3, and 5 years were 0.790, 0.809, and 0.815, as shown in Supplementary Figure S6B . In the GSE65904 external validation cohort, the values were 0.761, 0.756, and 0.738, respectively (by DeLong’s test; Supplementary Figure S6C ). Furthermore, the calibration plots showed a line segment near the 45 degree, suggesting a strong agreement between the prediction and observation in the internal validation cohort ( Supplementary Figure S6D ) and GSE65904 external validation cohort ( Supplementary Figure S6E ). And the DCA showed that Treg.Sig would offer a robust net benefit in internal validation cohort ( Supplementary Figure S6F ) and GSE65904 external validation cohort ( Supplementary Figure S6G ).

Over the last couple of years, there has been an increase in the development of various gene-expression signatures that have different predictive abilities. The predictive power of Treg.Sig was examined and compared to 51 signatures ( Supplementary Table S2 ), encompassing both lncRNA- and mRNA-derived signatures, to assess AUC at 1-, 3-, and 5-year OS. The 51 signatures showed connections to various biological characteristics, such as immune cell infiltration, cell death processes, cellular transformation, hypoxia, genetic modifications, N6-methyladenosine, etc. The Treg.Sig exhibited superior performance compared to all other signatures in predicting 1-year and 3-year OS in the GSE65904 dataset ( Figure 2 ).

The Treg.Sig was further validated in several SKCM subgroups stratified by different clinical characteristics, included genders, OS, and tumor stages. In particular, the elevated risk score correlated with higher clinical stage and OS status in SKCM ( Supplementary Figures S7A–C ). Furthermore, the OS was notably poorer in the high-risk cohorts compared to the low-risk cohorts, either in the early clinical stage (P=2.3×10−4) ( Supplementary Figure S7D ) or advanced clinical stage (P<0.0001) ( Supplementary Figure S7E ) SKCM subgroups. Clark’s level of invasion significantly predicts tumor invasion in melanoma (27). Likewise, we found that Treg.Sig exhibited a robustly predictive power, whether in Clark 1, 2, and 3 subgroups (P=0.0046) ( Supplementary Figure S7F ), or Clark 4 and 5 subgroups (P=0.0008) ( Supplementary Figure S7G ).

The T lymphocytes were crucial in the immune response against the tumor. Consequently, we examined how Treg.Sig is associated with immune cell infiltration and immunotherapy response in patients with SKCM. Initially, we examined how the signature is connected to 64 different types of immune cells that infiltrate. xCELL algorithms demonstrated that SKCM patients with high-risk scores had a higher proportion infiltration of mesenchymal stem cells, osteoblast, and NKT cells. Despite this, their composition included fewer aDC cells, B cells, CD4+ memory T cells, CD4+ naive T cells, CD8+ Tcm, CD8+ Tem, cDC cells, iDC cells, M2 macrophages, memory B−cells, monocytes, pDC cell, Th2 cells, and Tregs, suggesting an inflamed yet somewhat immunosuppressive environment ( Supplementary Figure S8A ).

In melanoma cancer, T cells that have infiltrated the tumor exhibit anti-tumor immune responses. Therefore, utilizing T cell-focused immunotherapy approaches could offer a novel opportunity to enhance cancer treatment (28). Thus, we conducted additional analysis to explore the connection between Treg.Sig and immunotherapy by examining the association between Treg.Sig and established immunotherapy markers. Initially, we utilized the TIDE online algorithm to forecast the likelihood of response to ICI therapy. While low-Treg.Sig patients exhibited higher TIDE scores (Chi-square test, P=0.0038; Supplementary Figure S8B ), this apparent paradox may reflect compensatory immune evasion mechanisms in tumors with baseline low immunosuppression. Considering the TIDE score’s predictive performance is context-dependent, we included the immune proportion score (IPS) in our analysis to provide a more comprehensive assessment, as it is a more precise biomarker for the protein expression of PD-L1. Predictably, patients with low risk exhibited notably elevated levels of PD-L1 protein expression compared to high-risk patients in the ctla4_neg_pd1_neg group (P=0.011) ( Supplementary Figure S8C ), ctla4_neg_pd1_pos group (P<0.0001) ( Supplementary Figure S8D ), ctla4_pos_pd1_neg group (P<0.0001) ( Supplementary Figure S8E ), and ctla4_pos_pd1_pos group (P<0.0001) ( Supplementary Figure S8F ). Overall, these findings suggested that Treg.Sig correlated with the infiltration of immune cells and the immune status of tumors in the tumor microenvironment.

Due to the immune-related characteristics of the Treg.Sig, we then tended to dig into its underlying mechanism of action. Initially, we conducted a correlation analysis to identify the genes that are highly associated with the Treg.Sig. As illustrated in the Supplementary Figures S9A and B , a total of 438 genes were found to have a positive relationship, while 28 genes had a negative relationship, including immune checkpoints such as CTLA-4, LAG-3, TIGIT, CD8A, CD163, and more. Subsequently, GSEA analysis was conducted on these genes using the R software package “clusterProfiler”, indicating their involvement in immune-related processes like allograft rejection, autoimmune thyroid disease, leishmaniasis, systemic lupus erythematosus, and type I diabetes mellitus ( Supplementary Figure S9C ). Taken together, our data indicated that Treg.Sig participates in the immune related response. Thus, studies to detect further potential value of Treg.Sig are necessary.

Subsequently, we examined the feasibility of utilizing Treg.Sig to inform systemic treatments. Initially, the pRRophetic algorithm was utilized to estimate IC50 values in order to forecast the varying responses to chemotherapy in Treg-associated high- and low-risk groups. According to the Cancer Genome Project (CGP) database, a decreased Treg.Sig risk score was associated with increased sensitivity to four different chemotherapy drugs (AP_24534, Pazopanib, Paclitaxel, Tivozanib) based on the Wilcoxon test results (all P< 0.0001) ( Supplementary Figures S10A–D ). Therefore, low Treg.Sig score was correlated with high chemotherapy sensitivity.

We then confirmed the effectiveness in various immunotherapy datasets, such as GSE161801 (carfilzomib), GSE120575 (anti PD1 or/and CTLA4), and GSE189125 (anti PD-1). Melanoma patients exhibited lower Treg.Sig score after new proteasome inhibitors therapy–carfilzomib (Wilcoxon test, P=1.2x10-7) in the GSE161801 scRNA-Seq cohort ( Supplementary Figures S10E, F ). Similarly, individuals who responded to anti PD1 or/and CTLA4 treatment (GSE120575) exhibited significantly lower Treg.Sig scores than those who did not respond (Wilcoxon test, P<0.0001), suggesting a negative correlation between Treg.Sig score and ICI effectiveness ( Supplementary Figures S10G, H ). In the GSE189125 anti PD-1 cohort, the durable clinical benefit (DCB) group showed significantly lower risk scores than the no durable benefit (NDB) group (Wilcoxon test, P=0.048), suggesting that patients with lower risk were more responsive to immunotherapy than those with higher risk ( Supplementary Figures S10I, J ). The findings indicated that individuals with low-Treg.Sig scores have a higher chance of responding well to immunotherapy, suggesting that Treg.Sig could serve as a valuable biomarker for identifying SKCM patients who could benefit from immunotherapy.

Among the 10 Treg.Sig genes, the crucial one STAT1 was predicted through PPI network, which was selected for further analysis ( Figure 3A ). The protein expression levels of STAT1 in SKCM were observed by analyzing the immunohistochemical staining images in the HPA database. The levels of STAT1 expression were found to be elevated in SKCM tissue compared to normal skin tissue, as shown in Figure 3B . Notably, patients harboring higher STAT1 status had better survival than those with a lower STAT1 status ( Figure 3C ). In order to delve deeper into the underlying mechanism, we examined the connection between STAT1 and tumor-infiltrating cells. Our analysis revealed a positive correlation between STAT1 and the infiltration of immune cells that target tumors, such as aDC, B cells, CD8 T cells, cytotoxic cells, and macrophage cells ( Figures 3D, E ). Adding further complexity, an analysis was conducted on immune checkpoint proteins such as PD-1, PD-L1, and PD-L2. Likewise, there was a strong connection found between STAT1 and the primary immune checkpoint markers (PD-1 R=0.650, P<0.001; PD-L1 R=0.829, P<0.001; and PD-L2 R=0.763, P<0.001) ( Figures 3F, G ). These results indicated the genetic and expression alteration landscape of STAT1 in melanoma, suggesting dysregulating STAT1 participates in tumor immune microenvironment in melanoma contexts.

Subsequently, we explored weather STAT1 mutation was correlated with ICI treatment response. Initially, we confirmed the mutation rate of STAT in different types of cancers. Melanoma ranked 1st with a mutation frequency of 24.32% among all 30 cancers, as depicted in ( Figure 4A ), followed by mature B-Cell neoplasms and Endometrial Carcinoma. Subgroup analyses revealed that the occurrence rate of STAT1 mutations in melanoma was 3.87% ( Supplementary Figure S11A ). Various clinical features of individuals in the TCGA cohort, such as age, sex, survival rate, and tumor mutation burden, were analyzed and are depicted in Figure 4B . The Lollipop plot illustrates a total of 15 missense mutations (green) and 1 truncating (black) mutation in 750 amino acids long STAT1 protein ( Figures 4C, D ). Importantly, there was a correlation between STAT1 mutation and response to ICI treatment, showing a higher response rate in patients with STAT1-Wt compared to those with STAT1-Mut ( Figures 4E, F ). In the meantime, the STAT1-Wt group exhibited higher rates of complete response (CR) and partial response (PR), while the STAT1-Mut group had a higher rate of progression disease (PD) as shown in Figure 4G .

The anticancer role of STAT1 is thought to be essential in the process of tumorigenesis, where it fosters the disruption of transformed cells by the immune cells (29). Thus, we assessed the effect of STAT1 blockage on anti-PD-1 efficacy in the context of B16-F10 melanoma-bearing mice, thereby elucidating whether STAT1 is associated with immunotherapy resistance. We first performed a molecular docking analysis to assess the affinity of Flu for its target, STAT1. The findings indicated that Flu attaches to STAT1 by forming observable hydrogen bonds and robust electrostatic interactions, resulting in a binding energy of -5.894 kcal/mol, suggesting a very stable attachment ( Figure 5A ). Significantly, we noticed increased suppression of tumor growth when using anti-PD-1 therapy ( Figure 5B ). Nevertheless, tumors in the Flu treatment group eventually acquired adaptive resistance to the anti-PD-1 treatment, reaching a primary tumor size comparable to the controls by the fourth week ( Figure 5B ). We isolated neutrophils from harvested tumor tissue by using MACS neutrophil isolation kit, and analyzed whether cell phenotype of neutrophils is affected after mice treated with Flu. As illustrated in Figures 5C, D , STAT1 phosphorylation of neutrophils in B16-F10 tumor tissues was blocked in mice treated with Flu. Strikingly, blockade of STAT1 by Flu reversed N1 polarization-induced elevation of iNOS, IL-1, and TNF-α mRNA levels in neutrophils ( Figure 5E ). Furthermore, in order to confirm these findings, we utilized siRNA technology to knock down the expression of STAT1. The results suggested that STAT1 siRNA significantly down-regulated the expression of STAT1 ( Figure 5F ), accompanied by immunofluorescence showing the expression of N1-tagged iNOS was also reduced ( Figure 5G ).

The above results which showing patients/mice with STAT1 mutation/inhibition harboring lower response rate during ICI treatment promote us to explore the underlying mechanisms of resistance. We suspected that malfunction of STAT1 may increase the susceptibility of neutrophils switching to a protumor suntype in melanoma tumor microenvironment (TME). We first validated whether STAT1 inhibition increase the infiltration of neutrophils. As shown in Figures 6A, B , during anti-PD-1 therapy, treatment of Flu in B16-F10 bearing mice did not alter the proportion of CD11b+Ly-6G+ neutrophils, indicating that STAT1 inhibition did not affect neutrophils infiltration into tumor tissue. As expect, treatment of Flu decreased the phosphorylation of STAT1 in neutrophils in B16-F10 melanoma tissue ( Figure 6C ). Surprisingly, flow cytometry analysis showed that treatment of Flu in B16-F10 melanoma bearing mice significantly promoted the phenotype switch of neutrophils toward a N2 subtype (Arg-1+) ( Figures 6D, E ). These results indicated that melanoma TME could skew more neutrophis to N2 subtype when STAT1 was inhibited. To further validated this, we used cell-cell coculture system to explore weather B16-F10 melanoma cells could drive more neutrophils toward N2 subtype in vitro. We isolated neutrophils from bone marrow and detected the expression of N2 marker after co-culture. The results showed neutrophils expressed Arg-1 after cocultured with B16-F10 melanoma cells, and this trend was intensified while STAT1 was inhibited by Flu ( Figure 6F ). Further investigations were conducted into the tissue-specific chemokine and cytokine profile influencing neutrophil recruitment and polarization. It was discerned that within Flu treated B16-F10 melanoma sites, there were no alteration in the levels of CCL2, CCL3, and CCL5, which are pivotal in mediating neutrophil recruitment ( Figure 6G ). Nevertheless, there were noticeable increased in the concentration of IL-4 and IL-13, cytokines crucial for the polarization towards an N2 neutrophil phenotype, whereas the levels of IFN-gamma, a cytokine essential for N1 neutrophil programming, were decreases ( Figure 6H ).

In order to study the immediate impact of these N2-neutrophils on tumor progression in vivo, we conducted the transfer of N2-neutrophils (Ly-6G+Arg-1+), which were extracted from the bone marrow of healthy mice without tumors and then co-cultured with B16-F10 melanoma cells ( Figure 7A ). We employed fluorescent tracer technology to track intra-tumor injected N2 neutrophils. As shown in Figure 7B , PKH67-labeled N2 neutrophils remained in the tumor tissue 1 day post injection. Surprisingly, in the context of anti-PD-1 treatment, administering Flu to mice with B16-F10 melanoma notably reduced their OS compared to the control group ( Figure 7C ). Furthermore, the decline in OS was exacerbated by the introduction of N2-neutrophils through adoptive transfer ( Figure 7C ). Flow cytometry analysis was used to examine the expression of pertinent N1 and N2-neutrophil markers in B16-F10 melanoma allografts. In our study, we observed that both Flu and adoptive transfer of N2-neutrophils led to a rise in N2-neutrophil (Arg-1+) levels in B16-F10 melanoma allografts when compared to the DMSO control group, which predominantly showed N1- neutrophil infiltration (iNOS) ( Figure 7D ). To figure out the underlying mechanisms that STAT1 regulated melanoma resistance, we analyzed the correlation between STAT1 and epithelial-mesenchymal transition (EMT). As an important mechanism by which cancer cells obtain highly invasive phenotype, EMT is a process acquiring mesenchymal features in epithelial cells. It was demonstrated that immune cells within the tumor microenvironment play a crucial role in influencing this plasticity of tumor cells (30). Correlation analyses were conducted to examine the impact of STAT1 on these processes, revealing a negative correlation between STAT1 and EMT markers ( Figure 7E ). To verify this result, we examined the expression of epithelial (E-cadherin) marker by immunohistochemical staining. As shown in Figure 7F , both Flu and adoptive transfer of N2-neutrophils resulted in elevated levels of E-cadherin expression.

Based on the aforementioned in vivo findings, we conducted in vitro studies to examine how N2 neutrophils impact melanoma EMT, invasion, and proliferation. In order to achieve this goal, melanoma B16-F10 cells were cultured with N2 neutrophils ( Figure 8A ). As shown in Figures 8B , C , B16-F10 cells co-cultured with N2 neutrophils showed higher epithelial markers (Ecadherin, EpCam, Krt18) expression but lower mesenchymal markers (Vimentin, N-cadherin, ZEB1). Additionally, N2 neutrophils increased the invasion of B16-F10 cells as shown in transwell assays ( Figures 8D, E ). Likewise, in CCK-8 experiments, it was noted that the presence of N2 neutrophils led to a notable enhancement in the proliferation ability of B16-F10 melanoma cells as opposed to the control group ( Figure 8F ). Collectively, these findings indicate that the absence of STAT1 is associated with the pro-tumor activity of N2 neutrophils and resistance to melanoma immune checkpoint inhibitors.