Identification of Co-Expressed Genes Promoting CD8+ T Cell Infiltration in Melanoma Tumor Microenvironments


 Purpose: To identify CD8+ T cell-related factors and the co-expression network in melanoma, to illustrate the interactions among CD8+ T cell-related genes in the melanoma tumor microenvironment.Method: We obtained melanoma and para-carcinoma 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 common results in the two cohorts. Results: Nine co-expressed genes (CCL5, GAP5, GZMA, GZMH, IRF1, LAG3, PRF1, NKG7, PSMB10) were identified as CD8+ T cell-related genes that promoted infiltration of CD8+ T cells, which were 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, and the tumor purity was higher, suggesting that these genes were immune protective factors that improved the prognosis of melanoma. Besides this, co-expression factors negatively correlated with angiogenesis factors and positively correlated with angiogenesis inhibitors. Conclusion: These nine co-expressed genes promote high levels of infiltration of CD8+ T cells by the cellular response to the interferon−gamma process. The mechanism might provide new pathways for treatment of patients who are insensitive to PD-1 immunotherapy due to low degrees of CD8+ T cell infiltration.


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 a higher survival rate; however, those with advanced metastatic melanoma had lower survival rates [1] . In recent years, the treatment of melanoma has included conventional chemotherapy and radiation therapy, as well as therapies targeting B-Raf proto-oncogene, serine/threonine kinase (BRAF) and MAP kinase-ERK kinase (MEK). However, many patients with advanced melanoma cannot bene t from it. The emergence of immune checkpoint inhibitors (ICI) 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 in immunotherapy [2] . PD-1, an immune checkpoint protein on T cells, combines with programmed cell death ligand 1 (PD-L1) to inhibit activity of lymphocytes in the peripheral circulation, thereby inhibiting in ammation 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 to express not only on the surface of T cells, but also on activated B lymphocytes and NK cells. Therefore, blocking PD-1 can rejuvenate disabled immune cells [3] .
Immune checkpoint inhibitors targeting cytotoxic T lymphocyte antigens, or PD-1 receptors, activate T cells and other immune cells, allowing the body's own immune system to attack melanoma [4] .
Despite the great success of immunotherapy, there are still many patients who are unresponsive to or resistant to PD-1 blockers, the reasons for it including: lack of T cells. Which is associated with endogenous oncogenic activation of WNT -β-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 key proteins of tumors [6] .
Among T lymphocytes, CD4 + T cells and CD8 + T cells play major roles in cellular immunity. CD4 + T cells, also known as helper T lymphocytes, assist humoral and cellular immunity. CD8 + T cells, also known as cytotoxic cells, kill target cells [7] . Most tumor cells express antigens that mediate the recognition of CD8 + T cells. In ltration of CD8 + T cells and high expression of the immune system have important clinical signi cance for immunotherapy of cancer suppression [8] . Therefore, in this paper, we explored the mechanism related to CD8 + T cells in ltration based on the generative credibility method to try to solve the PD-1 insensitivity caused by the low CD8 + T cells in ltration.
CIBERSORT packages [9] and weighted gene co-expression network analysis (WGCNA) were used to explore the co-expression factors and related biological functions associated with CD8 + T cells in ltration, and the results were veri ed in two different datasets [10] .
In this paper, we identi ed nine co-expressed genes that promote CD8 + T cells production in melanoma.
These nine genes are involved in the biological process of IFN-γ, suggesting that the biological process of IFN-γ is positively correlated with the in ltration level of CD8 + T cells. Next, we investigated the correlation between their expression and angiogenic factors, angiogenic inhibitors, tumor purity, in ammation and immune response, then veri ed that they inhibited tumor proliferation and differentiation, thereby con rming our hypothesis.
These results may help to improve the e cacy of immunotherapy in patients who are refractory to therapy due in part to low CD8 + T cells in ltration. This paper aims to determine whether promoting T lymphocyte in ltration may reduce drug resistance rate of PD-1 treatment.

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 P value less than 0.05 was considered statistically signi cant. Using expression data to estimate stromal cells and immune cells in tumor tissues (ESTIMATE) [12] , we evaluated the tumor purity of melanoma based on the ESTIMATE algorithm.
Weighted gene co-expression network analysis (WGCNA) We know that the gene expression patterns involved in the same pathway or biological process are similar [13] . WGCNA constructs a co-expression network by converting co-expression 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 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. To date, we obtained the co-expression network related to the content of CD8 + T lymphocytes [14] .

Protein-protein network and function enrichment
After obtaining the modules most related to CD8 + T cells, we used Pearson correlation coe cients 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 (KEGG) (https://www.genome.jp/kegg/) [15] and Gene Ontology (GO) (http://geneontology.org/) [16] 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 [17] .
Gene set enrichment analysis (GSEA) Gene Set Enrichment Analysis (GSEA) can determine the meaning and difference between two biological states in a prede ned data set [24] . In this article, we divide 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 in ltration-related co-expressed genes.

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

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 signi cance. The "survival," "ggplot2," "corrplot," "pheatmap," and "limma," packages were built using R version 3.6.3.
Differences with P < 0.05 were considered signi cant.

Identi cation of CD8 + T cell in ltration modules and genes
The ow chart of this study is illustrated in Fig. 1.
We obtained 202 TCGA-SKCM and 143 GSE65904 samples with complete clinical information and proportion of immune cell in ltration assessment after selection by P < 0.05. The ratio of 22 immune cells is shown in Fig. 2A, C. The tumor purity heat map is shown in Fig. 2B, D. We performed cluster analyses in TCGA-SKCM and GSE65904. The cluster heat map of 178 samples in TCGA-SKCM is shown in Fig. 3A, and the correlation coe cients between each phenotype and co-expression module are shown in Fig. 3B.
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 ndings, we supplemented the heat map of the correlation between the factors in the yellow module and CD8 + T cells and the yellow module ( Fig. 3C) (the module correlation coe cient was greater than 0.7; the CD8 + T cell correlation coe cient was greater than 0.4).
In the GSE65904 cohort, we followed the same steps to obtain the green-yellow co-expression module (Fig. 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 (Fig. 4A). A total of 15 coexpression models were generated, and the correlation coe cients between each phenotype and co-Page 6/28 expression module were calculated (Fig. 4B). The green-yellow module was selected as the most relevant module in GSE69504. The genetic signi cance of the rst 20 CD8 + T cells-related genes in the yellow module of TCGA-SKCM is shown in Table 1. The genetic signi cance of the rst 20 CD8 + T cells-related genes in the green-yellow module of GSE65904 is shown in Table 2. Next, we screened intersection factors between the two modules based on the module correlation greater than 0.7 and the CD8 + T cells correlation greater than 0.4. Later, by taking the intersection analysis of the genes screened at the two modules, we obtained 22 co-expression factors with the strongest correlation with CD8 + T cells (Fig. 5A).  Their protein-protein interaction network is shown in Fig. 5B. According to the degree of connection of the protein-protein network, we determined the nine most core factors. We conducted functional enrichment of these intersection factors, and found that the interferon response was the most signi cant, as shown in Fig. 5C.

Clinical analysis of CD8 + T cell in ltration-related genes
To calculate the correlations between these and CD8 + T cell in ltration proportions, subgroups were created according to the median of nine gene expression values in the TCGA-BLCA (Fig. 6A) and GSE65904 (Fig. 6B) cohorts. We found higher in ltration proportions in high expression groups (P < 0.05), suggesting that these genes promote CD8 + T cell in ltration and improves prognosis. Then, the various stages and statuses were applied to determine the prognosis level. For these nine genes, expression levels in 5-year mortality (Fig. 6C,D) groups were signi cantly lower (p < 0.05).

Immune microenvironment analysis
The biological functions of angiogenic factors and angiogenesis inhibitors are listed in Table 3. A negative correlation was found between several genes (CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1 and PSMB10) and angiogenic factors in both TCGA-SKCM and GSE65904 (Fig. 7A, D); positive correlations were found for angiogenic inhibitors (Fig. 7B, E), suggesting that these genes might modulate vascular changes in the tumor microenvironment. With the increase of CD8 + T cells in ltrationrelated gene expression, there was a decreasing trend of tumor purity (Fig. 7C, F). These ndings suggest that these genes might in uence bladder tumor purity and local microenvironment component proportions. To analyze the correlation between the eight genes and immune responses, we choose seven metagenes representing various types of in ammatory and immune responses. We found that CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1 and PSMB10 positively associated with seven of these clusters in both TCGA-SKCM (Fig. 8) and GSE65904 (Fig. 9).   (Fig. 10A, 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 (Fig. 10C). The P-values are displayed Table 4. We found that these biological pathways were immune-related and were involved in tumor immunity that protects against tumor in ltration.
Immune check point and HPA CD8 + T cells in SKCM-FPKM and GSE65904 were positively correlated with CD274 encoding PD-L1 and PDCD1 encoding PD-1 (Fig. 11A-B). In the Human Protein Atlas (HPA), the level of PSMB10 encoded by melanoma tissues was lower than normal tissues (Fig. 11C-F). These genes showed clinical phenotypic correlation at mRNA and protein levels.

Discussion
We focused on reducing drug resistance to PD-1 inhibitors due to lack of CD8 + T cells in ltration. This paper focused on the search for CD8 + T cells in ltration-related factors and explored in ltration mechanisms so as to provide new concepts for improving immunotherapy.
We selected samples from two data sets, calculated CD8 + T cell in ltration in melanoma tissue samples using the CIBERSORT package, identi ed co-expressed genes that promoted CD8 + T cells in ltration using the WGCNA algorithm, and performed preliminary screening based on intersection results.
We subsequently selected nine genes in the survival analysis that were most signi cantly associated with reduced mortality as key genes promoting CD8 + T cells in ltration. We performed functional enrichment in a biological process (BP) in which these genes were enriched in the IFN-γ mediated signaling pathway as well as in the antigen presentation and processing pathway. Through correlation analysis with angiogenic factors, angiogenic inhibitors and tumor purity, we veri ed that these nine genes inhibit tumor progression in the tumor microenvironment.
In GSEA analysis, the biopathways enriched in the high expression group of key genes were associated with tumor immunity. Key genes play important roles in CD8 + T cells in ltration in some cancers: CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1 and PSMB10 were identi ed as promoting factors for CD8 + T cell in ltration, all of which had independent prognostic effects.
CCL5 is a chemokine expressed by T lymphocytes, macrophages, platelets, synovium broblasts, renal tubular epithelial cells and tumor cells [26] . Harlin et al. found that CCL5 produced by melanoma cells helped chemokine recruitment of CD8 + effector T cells, and lack of key chemokines may limit CD8 + T cell migration, thereby limiting the effectiveness of anti-tumor immunity [27] . PRF1 is expressed by a protein called perforin, which is expressed by CD8 + T cells and NK cells and is essential for cell-mediated cytotoxicity and effective control of pathogens [28] . Taube et al. demonstrated that CD8A, PRF1 and CCL5 were overexpressed in PD-L1 + melanoma by QRT-PCR analysis and were involved in the activation of CD8 + T cells [29] . GZMA is the most abundant protease in cytotoxic particles of NK cells [30] . Inoue et al. examined mRNA levels of immune-related genes in melanoma patients before and after immunotherapy with nivolumab and found increased CD8 and GZMA expression levels [31] . IRF1 is an effective antiviral, antitumor and immunoregulatory protein. Many studies reported that IRF1 is closely related to tumor inhibition. CD8 + T cells in IRF1-de cient melanoma showed increased cytotoxicity, the expression of PD-L1 was up-regulated, and tumor growth was more easily restored [32] . LAG3 is the third alternative inhibitory receptor targeted in tumor microenvironment, which has attracted much attention and research in clinical trials [33] . Frohlich et al. found that LAG3 methylation in melanoma tissues was associated with CD8 + T cells in ltration and IFN-γ signaling, and they believed that LAG3 could be used as a substitute biomarker for cytotoxic anti-tumor response [34] . NKG7 is one of the most expressed genes in NK cells and is essential for cytotoxic degranulation of NK cells and CD8 + T cells as well as the activation and proin ammatory response of CD4 + T cells [35] . Through mapping T cell receptor cloning, Fairfax et al. found that patients who responded to immunotherapy had more CD8 + T cell-related large clones overexpressed gene NKG7 than did non-responders [36] . Although the relationship between GBP5, GZMH, PSMB10 and CD8 + T cells has not been clearly demonstrated in melanoma, some of them have been veri ed in other types of cancer. GBP5 and other immunomodulatory genes were signi cantly up-regulated in myeloid carcinoma of the colon, promoting the in ltration of CD8 + T cells in myeloid carcinoma [37] .
Nine co-expressed genes promoting CD8 + T cell in ltration were signi cantly enriched in the IFN-γ pathway, suggesting that IFN-γ may be closely related to CD8 + T cell in ltration. Dangaj et al. found that T cell in ltration required tumor cell-derived CCL5, and CXCL9 secretion by IFN-γ differentiated myeloid cells was ampli ed; in immunoreactive and immunoresponsive tumors, the synergistic effect of tumorderived CCL5 and IFN-γ-induced CXCR3 ligand secreted by bone marrow cells is the key to coordinate T cell in ltration [38] . with PD-L1 scores as a measure of alternative IFN-γ levels [39] . It was proposed that tumors with high neutrophil burdens are characterized by T cell response dullness, characterized by decreased expression of cytotoxic T cell genes such as CD8A, CD8B, GZMA and GZMB. There was decreased in ltration of CD3 + T cells and CD8 + T cells, and decreased expression of IFN-γ-related genes [40] . Lichtenegger et al.
found that blocking LAG-3 led to higher T cell activation and an increase in IFN-γ secretion compared with inhibition of other pathways; they concluded that the novel immune response was strongly enhanced by blocking LAG-3 or blocking both LAG 3 and PD-1 [41] . The promoter region of PSMB10 gene contains two IFN stimulus response elements, interregulated by IFN; this was con rmed by in vitro mutagenesis [42] .
Curtsinger et al. stimulated CD8 + T cells in mice using antigen-speci c B7-1, and found that they could rapidly produce a small amount of IFN-γ, with the production peaking at about 8 hours and decreasing after 24 hours [43] . When CD8 + T cells are exposed to mild temperatures, they promoted the production of speci c IFN-γ, which increased the lethality of tumor target cells [44] . Karachaliou et al. treated 21 melanoma patients with pembrolizumab and found that patients with high IFN-γ expression had signi cantly longer progress-free survival than those with low IFN-γ expression [45] .
A limitation of this paper is that we only included a total of 800 cases in two cohorts in two databases.
More samples are needed to validate the scienti c accuracy of the results. Although this paper is a comprehensive data analysis, it still needs to be further veri ed using in vitro experiments.
In summary, CCL5, GBP5, GZMA, GZMH, IRF1, LAG3, NKG7, PRF1 and PSMB10 are co-expression promotors of CD8 + T cell in ltration. Among them, CCL5, GZMA, IRF1, LAG3, NKG7 and PRF1 are closely related to CD8 + T cell in ltration in melanoma tissues. While there is no conclusive evidence that CD8 + T cell in ltration in melanoma tissues is closely related to GBP5, GZMH, PSMB10, it is worth exploring. The lack of CD8 + T cells in central tumor areas has become a major obstacle to immunotherapy for solid tumors, especially melanoma. Therefore, novel therapeutic strategies that promote the accumulation of CD8 + T cells in central tumor regions are urgently needed. We highlighted the complexity of immune checkpoint regulation in tumor microenvironments and identi ed several factors that may contribute to synergistic immunosuppression. These factors may be involved in the tumor escape mechanism of anti-PD-1/PD-L1 treatment; therefore, they can be considered for co-targeting in future immunotherapy.

Declarations
Ethics approval and consent to participate The data are from open databases and do not require ethical review.

Consent for publication
All authors agree to be published Availability of data and material The data set of this article was downloaded from the open source database TCGA and GEO.

Competing interests
The authors declare there are no competing interests.

Funding
No.

Authors' contributions
Kexin Yan conceived and designed experiments, downloaded and analyzed data, prepared charts and wrote drafts of the paper. Yutao Wang conceived and designed experiments, contributed analysis tools, analyzed the data and wrote drafts of the paper. Ting Xiao conceived and designed experiments, reviewed and revised this article. All authors approved the nal draft. CIBERSORT algorithm was applied to calculate the in ltration of CD8+ T cells in melanoma tissue samples from two datasets, WGCNA was used to generate a co-expressed gene network, and the results of the two datasets were intermingled to obtain 9 key genes. The correlation between key genes and tumor purity, angiogenic factors, angiogenic inhibitors and in ammatory factors was analyzed to verify the correlation. The main pathways of key gene enrichment were obtained in GSEA analysis. Survival analysis was used to determine the in uence of key genes on prognosis. To verify if the key genes also promote CD8+ T cells in ltration in other cancers.    In SKCM-FPKM, CD8+ T cells in ltration promotes the correlation between genes and in ammatory and immune responses. These reactions were induced mainly by hematopoietic cell kinases, immunoglobulin G, interferon, lymphocyte speci c kinase, major histocompatibility complex class I, major histocompatibility complex class II and activator of transcription 1.

Figure 9
In GSE65904, CD8+ T cells in ltration promotes the correlation between genes and in ammatory and immune responses. These reactions were induced mainly by hematopoietic cell kinases, immunoglobulin G, interferon, lymphocyte speci c kinase, major histocompatibility complex class I, major histocompatibility complex class II and activator of transcription 1.