9-Gene Signature Correlated With CD8+ T Cell Infiltration Activated by IFN-γ: A Biomarker of Immune Checkpoint Therapy Response in Melanoma

Purpose To identify CD8+ T cell-related factors and the co-expression network in melanoma and illustrate the interactions among CD8+ T cell-related genes in the melanoma tumor microenvironment. Method We obtained melanoma and paracancerous tissue mRNA matrices from TCGA-SKCM and GSE65904. The CIBERSORT algorithm was used to assess CD8+ T cell proportions, and the “estimate” package was used to assess melanoma tumor microenvironment purity. Weighted gene co-expression network analysis was used to identify the most related co-expression modules in TCGA-SKCM and GSE65904. Subsequently, a co-expression network was built based on the joint results in the two cohorts. Subsequently, we identified the core genes of the two most relevant modules of CD8+T lymphocytes according to the module correlation, and constructed the signature using ssGSEA. Later, we compared the signature with the existing classical pathways and gene sets, and confirmed the important prognostic significance of the signature in this paper. Results Nine co-expressed genes were identified as CD8+ T cell-related genes enriched in the cellular response to interferon−gamma process and antigen processing and presentation of peptide antigen. In the low expression level group, inflammation and immune responses were weaker. Single-cell sequencing and immunohistochemistry indicated that these nine genes were highly expressed in CD8+ T cells group. Conclusion We identified nine-gene signature, and the signature is considered as the biomarker for T lymphocyte response and clinical response to immune checkpoint inhibitors for melanoma


INTRODUCTION
Melanoma is one of the deadliest forms of skin cancer. In addition to external causes such as ultraviolet radiation, familial genetic factors cannot be ignored. Patients with melanoma diagnosed early and treated surgically had higher survival rates, and those with advanced metastatic melanoma had lower survival rates (1). In recent years, melanoma treatment has included conventional chemotherapy and radiation therapy and therapies targeting B-raf proto-oncogene, serine/threonine kinase (BRAF), and MAP kinase-ERK kinase. The emergence of immune checkpoint inhibitors has extended survival from metastatic melanoma. Ipilimumab, an anti-cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) antibody, and the programmed cell death protein (PD-1) antibodies nivolumab and pembrolizumab are now widely used for immunotherapy (2).
PD-1, an immune checkpoint protein on T cells, combines with programmed cell death ligand 1 (PD-L1) to inhibit the activity of lymphocytes in the peripheral circulation, thereby inhibiting inflammation and immune response. In the tumor microenvironment, PD-1 has a promoting effect on the immune escape mechanism. Unlike CTLA-4, PD-1 can be induced on the surface of T cells and activated B lymphocytes and NK cells. This fact suggests that blocking PD-1 can rejuvenate disabled immune cells (3). ICI targeting cytotoxic T lymphocyte antigens, or PD-1 receptors, activates T cells and other immune cells, allowing the immune system to attack melanoma (4).
Despite the great success of immunotherapy, there remain many patients unresponsive to or resistant to PD-1 blockers. The reasons for this resistance include lack of CD8 + T cells associated with endogenous oncogenic activation of the WNT-b-catenin pathway in melanoma cells (5), lack of PD-1 receptors, mutation and antigen presentation of tumor antigen, dynamic changes of immune microenvironment, genetic mutations, and epigenetic changes in essential tumor proteins (6).
CD4 + and CD8 + T cells participate in cellular immunity. The former, also known as helper T lymphocytes, assist humoral and cellular immunity. The latter, also known as cytotoxic cells, kill target cells (7). Most tumor cells express antigens that mediate the recognition of CD8 + T cells. Infiltration of CD8 + T cells and high expression of the immune system have important clinical significance for immunotherapy of cancer suppression (8). Therefore, in the present study, we explored the mechanism related to CD8 + T cell infiltration based on the WGCNA method to identify the gene signature.
CIBERSORT packages (9) and weighted gene co-expression network analysis (WGCNA) were used to determine the coexpression factors and related biological functions associated with CD8 + T cell infiltration. The results were validated in two datasets (10).
We identified nine co-expressed genes correlated with CD8 + T cell infiltration in melanoma, possibly involved with the biological process of IFN-g. Next, we measured the correlations between the 9-gene signature expression with angiogenesis, wound healing, and immune responses. Finally, we determine that these factors inhibited tumor proliferation and differentiation, thereby confirming our hypothesis.
The proposed 9-gene signature contains genes related to T cell activation, cytotoxicity, and the regulation of IFN-g downstream, which may indicate that this signature can affect CD8 + T cells through IFN-g activation in melanoma. In addition, by comparing several classical scores previously proposed, we found that 9-gene signature not only has a good predictive ability for the prognosis and survival of melanoma patients, but also can be used as one of the biomarkers to evaluate the response to immunotherapy.

CD8 + T Cell Proportions and Tumor Purity
CIBERSORT is an algorithm for analyzing the proportion of cells in tumor tissues (9). LM22 represents the genetic markers of 22 immune cell subtypes; it was downloaded from the CIBERSORT website portal (https://cibersort.stanford.edu/). Based on this method, we analyzed the proportion of CD8 + T cells, and a P value less than 0.05 was considered statistically significant. Using expression data to estimate stromal cells and immune cells in tumor tissues (ESTIMATE) (16), we evaluated the tumor purity of melanoma based on the ESTIMATE algorithm.

Weighted Gene Co-Expression Network Analysis
We know that the gene expression patterns involved in the same pathway or biological process are similar (17). WGCNA constructs a co-expression network by converting coexpression correlation into connection weight or topological overlap value. Based on this background, we used WGCNA to construct a co-expressed gene network related to the relative content of CD8 + T cells (10). We normalized the TCGA-SKCM data with Log2(exp+1). We set the soft threshold to 5, R square = 0.89 in TCGA-SKCM, R square = 0.88 in GSE65904, and the number of genes in the smallest module to 30.

Protein-Protein Network and Functional Enrichment
After identifying the modules most related to CD8 + T cells, we used Pearson correlation coefficients and module correlations to screen CD8 + T cell-related genes. We performed an intersection analysis in the two cohorts based on module correlations greater than 0.7 and CD8 + T cells correlations greater than 0.4. Cytoscape software was used to generate the protein-protein interaction network of the intersection factor. Kyoto Encyclopedia of Genes and Genomes (https://www.genome.jp/ kegg/) (18) and Gene Ontology (GO) (http://geneontology.org/) (19) were used to determine the biological functions of these intersection factors. These methods were implemented in the database DAVID, v6.8 for annotation, visualization, and integrated discovery (20).

Immune Microenvironment Correlation Analysis
In GSE65904 and TCGA-SKCM, we investigated the correlation between CD8 + T cell-related genes and immune responses in melanoma, then drew the heat map. We selected several classic immune-related metagene sets, including major histocompatibility complex class I (MHC-1), major histocompatibility complex class II (MHC II), lymphocyte specific kinase (LCK), hematopoietic cell kinase (MCK), immunoglobulin G (IgG), signal transduction and activation transcription 1 (STAT1) and interferon (21).

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) can determine the meaning and differences between two biological states using a predefined data set (22). We divided the gene matrix in TCGA into high expression and low expression groups according to the median of the proportion of CD8 + T cells. We then determined the mechanism of CD8 + T cell infiltration-related coexpressed genes.

Timer
The Tumor Immune Estimation Resource (TIMER; https:// cistrome.shinyapps.io/timer/) (23) was used to identify the correlations between CD8 + T cells and other types cancer. Pearson correlation coefficients higher than 0.4 were considered significant.

Single-Cell Sequencing
Single-cell cohort data sets GSE72056 were analyzed to determine the correlation between CD8 + T cell and coexpression genes. GSE72056 contains 4645 single-cell sequencing samples, and the platform is GPL18573 Illumina NextSeq 500 (Homo sapiens). We used the "seurat" package to filter the sequencing data, batch calibration, and data standardization. Then, we used the "t-SNE" package to perform cell cluster analysis. Finally, we used "SingleR" to perform cell subpopulation annotation identification. By measuring the expression level of each cell, we demonstrated correlations between CD8 + T cells and co-expressed genes.
ssGSEA Based on the above methods, we screened a gene signature, then ssGSEA was used to score the gene set. Based on these scores, we explored the prognostic value of gene signature in different melanoma cohorts and in different types of cancers. Meanwhile, in TCGA-SKCM, the prognostic significance of ssGSEA score of the gene signature was compared with the more classical gene sets at present. For instance, CYTOLYTIC score by Rooney et al (24), hallmark angiogenesis gene set and hallmark interferon gamma response gene set from the molecular signature database (25), immune response-related gene sets, CD4 + T cell score and CD8 + T cell score by CIBERSORT.

Immunohistochemical Verification
Archives were registered from recruited melanoma patients, and resection and biopsies of freshly collected melanoma and paratumor tissue were performed. All participants provided written informed consent. The Ethics Committee of Fuzhou No. 1 Hospital Affiliated with Fujian Medical University approved the study. The extracted human tissues were fixed with 4% formaldehyde buffer. Deparaffinized specimens were then sectioned into 4-µm sections. Tissue sections were incubated at 60°C for 2 h before dewaxing, autoclaved at 115°C for 3 min for antigen retrieval in a citric acid buffer (pH 6.0), and quenched for endogenous peroxidase activity with 0.3% H 2 O 2 solution for 15 min. Then, the sections were blocked with normal goat serum for 45 min and incubated with the specific primary antibody against PSMB10 (dilution 1:200) overnight at 4°C. Subsequently, the sections were treated with goat anti-mouse secondary antibody for 30 min at room temperature. Protein expression was visualized using 3,3′-diaminobenzidine. Images were captured using a Nikon Eclipse 80i microscope (Nikon Corporation). The immunohistochemical integral optical density (IOD) was analyzed by ImageJ. Fifteen pairs of IOD with PSMB10 expression levels in melanoma tumors and paracancerous tissues were submitted to GraphPad Prism. The IOD sum and area of each photo were measured, and the average optical density between groups was compared.

Statistical Analysis
Statistical analysis was carried out using GraphPad Prism 8 and R 3.6.3 (https://www.r-project.org/). The subgroups were divided based on the median value. Kaplan-Meier survival analysis was used to generate overall survival curves, and the log-rank test was used to calculate the significance. The nine gene score was evaluated using single-sample gene set enrichment analysis (ssGSEA). The "survival", "ggplot2", "corrplot", "pheatmap", "GSVA" and "limma" packages were built using R version 3.6.3. Differences with P < 0.05 were considered significant.

Identification of CD8 + T Cell Infiltration Modules and Genes
The flow chart of this study is illustrated in Figure 1.  Figure 2A, and the correlation coefficients between each phenotype and co-expression module are shown in Figure 2B. We found that the yellow module had the strongest correlation with CD8 + T cells in the TCGA-SKCM cohort (Cor = 0.69; P= 1e-26 ). Based on these findings, we supplemented the heat map of the correlation between the yellow module and CD8 + T cells and the yellow module ( Figure 2C) (the module correlation coefficient was greater than 0.7; the CD8 + T cell correlation coefficient was greater than 0.4).
We used the dynamic hybrid cutting method to build a hierarchical clustering tree. Each leaf on the tree represents a gene, and each branch represents a co-expression module in GSE65904 ( Figure 2D). A total of 15 co-expression models were generated, and the correlation coefficients between each phenotype and co-expression module were calculated ( Figure 2E). The heat map of the correlation between the green-yellow module and CD8 + T cells ( Figure 2F) (the module correlation coefficient was greater than 0.7; the CD8 + T cell correlation coefficient was greater than 0.4). The genetic significance of the first 20 CD8 + T cells-related genes in the yellow module of TCGA-SKCM is shown in Supplementary  Table 2. The genetic significance of the first 20 CD8 + T cellsrelated genes in the green-yellow module of GSE65904 is shown in Supplementary Table 3. Next, we screened intersection factors between the two modules based on the module FIGURE 1 | Flowchart for identification of co-expressed genes promoting CD8 + T cells infiltration. SKCM-FPKM contains 470 melanoma cancer tissue samples. GSE65904 contains 214 melanoma samples. CIBERSORT algorithm was applied to calculate the infiltration of CD8 + T cells in melanoma tissue samples from two datasets. WGCNA was used to generate a co-expressed gene network, and the two datasets were intermingled to obtain nine essential genes. The main pathways of critical gene enrichment were obtained in GSEA analysis. Survival analysis was used to determine the influence of critical genes on the outcome. Single-cell sequencing and immunohistochemical experiments confirmed that the expression of these genes was high in CD8 + T cells. The comparison between 9-gene signature with angiogenesis geneset, wound healing geneset, immune response geneset and IFN-g genset were analyzed. Follow-up cohorts with immunotherapy were also used to validate the role of critical genes in immunotherapy. The 9-gene signature that related to CD8 + T cell infiltration in other cancers were verified.  correlation greater than 0.7 and the CD8 correlation greater than 0.4. Twenty-two factors were obtained ( Figure 3A). Subsequently, in the 22 genes we identified core genes that could represent TCGA-SKCM and GSE65904 co-expressed modules. These core genes are defined as the correlation of module greater than 0.8 in both TCGA-SKCM and GSE65904 cohorts ( Table 1). Their protein-protein interaction network is shown in Figure 3B. The nine core genes were marked as orange circles. We conducted functional enrichment of these intersection factors and found that the interferon response was the most significant, as shown in Figure 3C.
To calculate the correlations between these and CD8 + T cell infiltration proportions, subgroups were created according to the median of nine gene expression values in the TCGA-SKCM ( Figure 3D) and GSE65904 ( Figure 3E) cohorts. We found higher infiltration proportions in high expression groups (P < 0.05), suggesting that these genes are related to CD8 + T cell infiltration.

Immune Microenvironment Analysis
We chose seven metagenes to analyze the correlation between the eight genes and immune responses, representing various inflammatory and immune responses. We found that CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1, and PSMB10 positively correlated with seven of these clusters in both TCGA-SKCM ( Figure 4) and GSE65904 ( Figure 5). The above results demonstrated that nine core genes were associated with stronger T cell response and immune response. The genes in the groups were uploaded in Supplementary Table 4.

Survival Analysis and GSEA Analysis
To analyze their influence on cancer specific survival, we performed survival analysis.  Figures 6A, B). GSEA analysis showed that the T cell receptor signaling pathway, antigen processing and presentation, chemokine signaling pathway, and nature killer cell-mediated cytotoxicity were related to the high expression group ( Figure 6C). The P-values are displayed in Table 2. We found that these biological pathways were immune-related and were involved in tumor immunity.

Single-Cell Analysis
To determine which cells these genes came from, we conducted a single-cell analysis. After annotating the subsets with "SingleR", we compared the results then found that the expression content of PSMB10, GZMA, GZMH, PRF1, CCL5 was relatively high in the CD8 + T cell subsets, confirming our previous finding in TCGA-SKCM (Figure 7). The expression and distribution of CD8A, CD4, CD3E, and CD14 are shown in Supplementary  Figure 2. The figure shows that these genes are more likely derived from CD8 + T cells than other immune cells.

9-Gene Signature Cancer Prognosis Prediction
In the above study, we identified a 9-gene signature correlated with CD8 + T cells in melanoma. Based on the median grouping of 9-gene   signature ssGSEA score, we figured out that patients with high expression of 9-gene signature in the TCGA and GSE65904 melanoma cohorts had a better prognosis. At the same time, we also demonstrated that patients in the group with high expression of 9-gene signature had better prognosis in multiple cancers of TCGA ( Figure 8). The above results proves that patients with higher expression of the gene signature correlated with CD8 + T cell infiltration activated by IFN-g have a better prognosis and survival.
To further explore the value of this 9-gene signature, we analyzed its correlation with angiogenesis related pathways, wound healing related pathways, IFN-g related pathways, etc. It turns out that the 9-gene signature was positively correlated with immune response ( Table 3) angiogenesis, IFN-g response, IFN-g production, positive regulation of IFN-g production, positive regulation of response to IFN-g, wound healing, wound healing involved in inflammatory response, positive regulation of vascular wound healing, etc ( Figure 9A). Later, through univariate analysis of 9-gene signature with CYTOLYTIC score, hallmark angiogenesis gene set, hallmark interferon gamma response gene set, CD4 + T cell CIBERSORT score and CD8 + T cell CIBERSORT score, we found that our 9-gene signature had more significant prognostic value ( Figure 9B).

Validation of Immunotherapy Cohort
Based on these findings, we believe that the 9-gene signature is related to immunotherapy. We validated this in GSE78220, GSE91061 and GSE93157, which included immunotherapy follow-up data. In GSE91061 cohort, by comparing the 9-gene signature ssGSEA scores of the patients in the complete response (CR)/partial response (PR) group and the patients in the progressive disease (PD) group, it can be seen that the differences were statistically significant ( Figure 9C). Then, we verified the survival predictive ability of 9-gene signature in GSE78220 and GSE93157 cohorts, and found that patients with low expression of 9-gene had a worse prognosis on immunocheckpoint inhibitor therapy (Supplementary Figure 3). This suggests that 9-gene signature may be one of the prognostic indicators for immune checkpoint therapy.

Immunohistochemical Pathological Analysis
To investigate the clinical significance of PSMB10 in melanoma, immunohistochemical analysis of the protein expression levels of PSMB10 in and around carcinomas was performed in the melanoma cohorts of our hospital. One patient was a 58-year-old woman with lesions on the distal left middle finger, and one was a 72-year-old man with lesions on the right plantar skin. On gross observation, it was found that the color of the lesion skin was black and brown, with irregular shape with fuzzy boundary, and a few skin lesions were accompanied by necrosis, erosion, and ulcer.
Immunohistochemical detection showed that PSMB10 was positive in the two malignant melanoma specimens, a positive rate of 100%. PSMB10 had lower staining intensity in cancer tissue, while PSMB10 had higher staining intensity in normal tissue ( Figures 10A-D ). At the same time, we supplemented 15 quantitative immunohistochemical analyses of melanoma and paracancerous tissues, and the results showed that the IOD/Area of PSMB10 in paracancerous samples was higher ( Figure 10E). We then combined the PSMB10 gene expression and CD8 + T lymphocytes to draw the survival curves of multiple groups. Patients with high levels of CD8 + T lymphocytes and PSMB10 had the best prognosis, while patients with low levels of CD8 + T lymphocytes and PSMB10 had the worst prognosis ( Figure 10F).

DISCUSSION
In this research, we selected samples from two data sets, calculated CD8 + T cell infiltration in melanoma tissue samples using the CIBERSORT package, identified co-expressed genes related to CD8 + T cell infiltration using the WGCNA algorithm, and performed preliminary screening based on intersection results. We subsequently selected nine genes that the value of module correlation rank the top in the WGCNA analysis. We performed functional enrichment in a biological process in which these genes were enriched in the IFN-g mediated signaling pathway and the antigen presentation and processing pathway. The obtained 9-gene signature was confirmed to be related to the activation of IFN-g in the correlation analysis of the immune response related gene set and IFN-g related gene set. In GSEA analysis, the pathways enriched in the high expression group of critical genes were associated with tumor immunity. Essential genes play essential roles in CD8 + T cells infiltration in some cancers: CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1, and PSMB10 were identified as promoting factors for CD8 + T cell infiltration, all of which had independent prognostic effects.   Harlin et al. found that CCL5 produced by melanoma cells helped chemokine recruitment of CD8 + effector T cells (26). PRF1 is expressed by a protein called perforin, expressed by CD8 + T cells, and is essential for cell-mediated cytotoxicity and effective control of pathogens (27). Taube et al. demonstrated that CD8A, PRF1, and CCL5 were overexpressed in PD-L1 + melanoma and were involved in the activation of CD8 + T cells (28). GZMA is the most abundant protease in cytotoxic particles of NK cells (29). Inoue et al. examined mRNA levels of immunerelated genes in melanoma patients before and after immunotherapy with nivolumab and found increased CD8 and GZMA expression levels (30). IRF1 is an effective antiviral, antitumor, and immunoregulatory protein. CD8 + T cells in IRF1-deficient melanoma showed increased cytotoxicity, the expression of PD-L1 was upregulated, and tumor growth was more easily restored (31).
LAG3 is the third alternative inhibitory receptor targeted in the tumor microenvironment (32). Frohlich et al. found that LAG3 methylation in melanoma tissues was associated with CD8 + T cells infiltration and IFN-g signaling (33). NKG7 is one of the most highly expressed genes in NK cells and is essential for cytotoxic degranulation of NK cells and CD8 + T cells and the activation and pro-inflammatory response CD4 + T cells (34). Fairfax et al. found that patients who responded to immunotherapy had more CD8 + T cell-related large clones overexpressed gene NKG7 than did nonresponders (35). GBP5 and PSMB10 have been reported in tumors as part of the downstream interferon genes and are also considered to be the coordinators of tumor disease immunity (36,37). Li et al. reported that GZMH is a T cell gene in single-cell analysis of T cells in melanoma (38).
Nine co-expressed genes correlated with CD8 + T cell infiltration were significantly enriched in the IFN-g pathway, suggesting that IFN-g may be closely related to CD8 + T cell infiltration. Dangaj et al. found that T cell infiltration required tumor cell-derived CCL5, and CXCL9 secretion by IFN-g differentiated myeloid cells was amplified; in immunoreactive and immune-responsive tumors, the synergistic effect of tumor-derived CCL5 and IFN-g-induced CXCR3 ligand secreted by bone marrow cells is the key to coordinate T cell infiltration (39). Jia et al. analyzed the distribution of tumor-infiltrating T cells and the expression of PD-L1 in the orthotopic murine glioma model; GBP5, IRF1, as IFN-g-induced genes were positively correlated with PD-L1 scores as a measure of alternative IFN-g levels (40). It was proposed that tumors with high neutrophil burdens are characterized by inadequate T cell responses, with decreased cytotoxic T cell genes such as CD8A, CD8B, GZMA, and GZMB. There was decreased infiltration of CD3 + T cells and CD8 + T cells and decreased expression of IFN-g-related genes (41). Lichtenegger et al. found that blocking LAG-3 led to higher T cell activation and an increase in IFN-g secretion compared with inhibition of other pathways; they concluded that the novel immune response was strongly enhanced by blocking LAG-3 or blocking both LAG3 and PD-1 (42). The promoter region of the PSMB10 gene contains two IFN stimulusresponse elements cross-regulated by IFN; this was confirmed by in vitro mutagenesis (43). Curtsinger et al. stimulated CD8 + T cells in mice using antigenspecific B7-1. They found that they could rapidly produce a small amount of IFN-g, with the production peaking at about 8 hours and decreasing after 24 hours (44). When CD8 + T cells are exposed to mild temperatures, they promoted the production of specific IFN-g, which increased the lethality of tumor target cells (45). Karachaliou et al. treated 21 melanoma patients with pembrolizumab and found that patients with high IFN-g expression had significantly longer progress-free survival than those with low IFN-g expression (46).
In summary, CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1, and PSMB10 are co-expression genes related to CD8 + T cell infiltration. The lack of CD8 + T cells in central tumor areas has become a significant obstacle to immunotherapy In the immunohistochemical results, the PSMB10 expression level in normal tissue was higher than that of tumor tissue. (E) The IOD/area of paracancerous tissue was significantly higher than that of melanoma group. (F) In both groups with low PSMB10 expression and high PSMB10 expression, the group with low CD8 + T cell content had worse prognosis than the group with high CD8 + T cell content. In the low or high CD8 + T cell content groups, the prognosis of the low PSMB10 expression group was worse than that of the high PSMB10 expression group.
for solid tumors, especially melanoma. Therefore, novel therapeutic strategies that promote the accumulation of CD8 + T cells in central tumor regions are urgently needed.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. The data set of this article was downloaded from the open source database TCGA and GEO. TCGA: (http://cancergenome.nih.gov/); GEO: GSE65904, GSE78220, GSE72056, GSE91061 and GSE93157 were downloaded from the GEO database (http://www.ncbi. nlm.nih.gov/geo/).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Medical Ethics Committee of Fuzhou No. 1 Hospital Affiliated with Fujian Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
KY and YW conceived ideas, designed methods, analyzed and visualized the data, and wrote the draft paper. YL contributed analysis tools and analyzed the data. ZY managed the project, supervised and led the experiment, and reviewed and revised the paper. All authors contributed to the article and approved the submitted version.