Integrated Analysis of the ETS Family in Melanoma Reveals a Regulatory Role of ETV7 in the Immune Microenvironment

The ETS family modulates immune response and drug efficiency to targeted therapies, but their role in melanoma is largely unclear. In this study, the ETS family was systematically analyzed in multiple public data sets. Bioinformatics tools were used to characterize the function of ETV7 in melanoma. A prognostic model was constructed using the LASSO Cox regression method. We found that ETV7 was the only differentially expressed gene with significant prognostic relevance in melanoma. Enrichment analysis of seven independent data sets indicated ETV7 participation in various immune-related pathways. ETV7 particularly showed a strong positive correlation with CD8+ T cell infiltration. The prognostic model based on ETV7 and its hub genes showed a relatively good predictive value in training and testing data sets. Thus, ETV7 can potentially regulate the immune microenvironment in melanoma.


INTRODUCTION
Melanoma, originating from pigment-producing melanocytes, comprises 75% of deaths related to skin cancer (1). In 2019, a total of 96,480 new cases of melanoma of the skin reported were in the United States (2). At its early stage, the disease could be surgically removed, with a relatively good prognosis; however, when it spreads to distant organs, the 5-year survival rate of patients with this disease sharply declines to roughly 10%, as indicated by a recent review (3). Targeted therapies and immune checkpoint inhibitors (ICIs) have revolutionized the treatment of metastatic melanoma. However, numerous problems still need to be addressed. For instance, current evidence suggests that 40 to 65% of patients with advanced-stage metastatic melanoma show minimal or no RECIST response to ICIs at the outset (de novo resistance), and 43% of responders acquire resistance by 3 y (4,5). Moreover, melanoma patients receiving mitogen-activated protein kinase (MAPK)-targeted therapy also develop resistance, which leads to relapse (6). Several underlying mechanisms have been revealed, such as NRAS mutations and BRAF amplification, but 40% of them present unknown resistance mechanisms beyond genetic alterations (6)(7)(8).
Recent studies find that transcriptional factors (TFs) relevant to the MAPK pathways modulate drug efficiency to MAPK inhibitors (9,10). For instance, overexpression of ETV1, ETV4, or ETV5 can sufficiently restore cell proliferation in the presence of trametinib, a MEK inhibitor (9). An early study also observes that overexpression of ETV1 confers resistance to MAPK inhibitors in BRAF-mutant melanoma (11). ETV1, ETV4, and ETV5 are members of the E26 transformation-specific or Etwenty-six (ETS) family, which is one of the largest transcriptional factor families and participates in multiple biological processes such as cellular differentiation, cell cycle control (12), cell migration (13), and cell proliferation (14). The family consists of 28 TFs in humans and binds to similar DNA sequences of MAPK (15). Various studies have revealed that aberrant expression of many ETS family members is strongly associated with tumor initiation, progression, and metastasis in cancer (16)(17)(18). Some ETS inhibitors, such as VPC-18005 and BRD32048, exert antitumor effects in pre-clinical studies (19,20). In skin melanoma, the function of the ETS family is poorly studied. ETS1 can promote aggregation and invasion of melanoma (21,22), and ETV2 is required during tumor angiogenesis (23); meanwhile, the role of the rest of the ETS family members is largely unclear.
In addition, the ETS family has been shown to play a role in immunity (24)(25)(26). ELF1, ETV4, ETV3L, ETS1, and ETS2 can upregulate the expression of leukocyte-associated immunoglobulinlike receptor-1, which inhibits the maturation, differentiation, and activation of immune cells (25,27). ELF4 can inhibit the proliferation of naive CD8+ T cells by increasing the expression of KLF4 (24,28). However, the regulation of immune response in cancer by ETS family members, as well as the mechanism underlying this regulation, has yet to be determined.
In the current study, we systematically assessed the expression profile, prognostic significance, and role of ETS family members in human skin cutaneous melanoma (SKCM) by integrating data from The Cancer Genome Atlas (TCGA) database, the Genotype-Tissue Expression (GTEx) Project, Oncomine database, cBioPortal database, and Gene Expression Omnibus (GEO) database. Our results indicate that ETV7 is markedly downregulated in melanoma, and low ETV7 expression is related to the poor prognosis of SKCM patients. Enrichment analysis and immune profile analysis indicate that ETV7 may regulate the differentiation and activation of T cells. This finding suggests a previously unrecognized involvement of ETV7 in immunity and the potential of the gene to modulate the immune response to the disease.

Data Acquisition and Processing
The workflow of this study is presented in Figure 1. The TCGA_SKCM data set and GEO data sets (GSE65904 and GSE19234) were acquired and processed using the method described in our previous study (29). After data filtering was conducted, 448 tumor samples with survival data in the TCGA_SKCM data set, 210 tumor samples in GSE65904, and 44 samples in GSE19234 were used for further analysis. The single-cell RNA sequencing data (GSE72056) were downloaded and analyzed directly in the current study because the data had been processed and normalized (30). All data used in the present study were acquired from public databases, requiring no further approval by an ethics committee.

Online Database Analysis
The Oncomine database (https://www.oncomine.org/resource/ login.html) was used to validate the transcription level of the genes of interest in skin melanoma by retrieving expression data (log2-transformed) in three cohorts of melanoma vs. normal tissues for statistical comparison, with the following default thresholds: P-value < 1E-4, fold change >2, and gene ranks in the top 10%.
The GEPIA2 (http://gepia2.cancer-pku.cn/) database was also used to compare the transcript level of the genes of interest between melanoma tumors and normal samples from TCGA and GTEx projects.
The cBioPortal for Cancer Genomics (http://www.cbioportal. org/) was used to identify genes that exhibit a strong positive correlation (R > 0.5, q-value < 0.05) with ETV7 in all available melanoma-related data sets with gene expression data.
The human protein atlas (https://www.proteinatlas.org/) was used to analyze the expression of ETV7 across a set of normal human cells or tissues.

Functional Analysis and Enrichment Analysis
The co-expressed genes of ETV7 shared in five independent data sets were constructed into a protein-protein interaction (PPI) network in the STRING database (http://string-db.org). Cytoscape version 3.7.2 was used to visualize these networks. Hub genes were identified as the top 10 nodes, based on the score generated by the cytoHubba plugin in Cytoscape (ranked by degree).
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the clusterProfiler package in R version 3.6.2 (18).
Gene set enrichment analysis (GSEA) was used to investigate pathways enriched in the high-and low-ETV7 groups. C2.cp.kegg.v7.1.symbols.gmt was chosen as the gene set database. The pathways were considered significantly enriched with the following criteria: nominal p-value < 0.05, false discovery rate qvalue < 0.25, and normalized enrichment score > 1.

Immune Profile Analysis
The immune and stromal scores of each sample were estimated using the ESTIMATE algorithm in the 'estimate' package in R version 3.6.2 (31). Immune cell infiltration in each sample of the TCGA_SKCM data set was conducted by preparing and uploading gene expression data into the TIMER2.0 website (http://timer.cistrome.org/) in accordance with the instruction on the website (32).

Construction of the Prognostic Model
The least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted using the glmnet package in R. The analysis generated key gene signatures and their corresponding coefficients, which were used to calculate the novel score, as follows: score = −0.08241781*CCR5 − 0.01208383*IFNG − 0.13002695*TBX21 − 0.04765633*CXCL10 − 0.33444001*CXCR3-0.05630411*CCL5. To facilitate the interpretation of results from different data sets, the risk score was calculated by subtracting the minimum score of the cohort from this score, and dividing the difference by the absolute value of the maximum score of the cohort-that is, namely, risk score = (score − Min)/absolute(Max).

Statistical Analysis
The data collected were analyzed by default as described using web resources. The remaining data were analyzed using R version 3.6.2. The median expression of genes or the median value of the risk score was used as cutoff value in dividing patients into two subgroups. Univariate Cox regression and multivariate Cox regression were conducted using the survminer package in R. The survival analyses were compared using the Kaplan-Meier method with the logrank test. Time-dependent receiver operator characteristic (ROC) analyses and subsequent calculation of the area under the curve (AUC) were performed using the timeROC package in R. Correlation analysis between the gene expression of ETV7 and that of the remaining genes was conducted by R software version 3.6.2 with spearman method. Wilcoxon test was conducted to compare gene expression between groups. Packages in R used for data analysis and graph plotting included ggplot2, ggpubr, limma, vennDiagram, tidyverse, rms, org.Hs.eg.db, dplyr, Rtsne, and plyr. P < 0.05 was considered statistically significant (*, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001).

Identification of Differentially Expressed ETS Family Members in Melanoma
The TCGA_SKCM project has only one normal sample. Thus, we first investigated the expression of all ETS family members in melanoma and normal skin by manipulating the data from the GTEx database. GEPIA2 integrates gene expression data from the TCGA and GTEx projects via a standard processing pipeline (33). Using the default cutoff value of GEPIA2 (|Log2FC| cutoff = 1, and q-value cutoff = 0.01), we found that 13 of the 28 family members were differentially expressed ( In particular, ELK3, ETS1, ETV1, ETV4, ETV5, and SPI1 were significantly upregulated in melanoma (Figures 2A-F), whereas EHF, ELF3, ELF5, ERG, ETS2, ETV7, and SPDEF were significantly downregulated in the tumor ( Figures 2G-M). Further, the expression profiles of ETV1, ETV5, ETV7, ELF3, ETS1, ETS2, ELK3, ELF5, and SPDEF were validated by comparing their transcriptional level between melanoma and normal tissues across three independent cohorts from the Oncomine database ( Figures 2N-V). However, the remaining members were not differentially expressed between tumor and normal tissues (Supplementary Figures 1P-S).

ETV7 as an Independent Prognostic Predictor in Melanoma
We subsequently evaluated the prognostic significance of the above nine differentially expressed ETS family members in TCGA_SKCM. Analysis suggested that EVT7 exhibited prognostic significance in melanoma patients ( Figures 3A-H). Moreover, melanoma patients with high ETV7 showed significantly longer disease-specific survival (DSS) than those with low ETV7 (p < 0.0001, Figure 3I). We further evaluated the prognostic relevance of these nine differentially expressed genes (DEGs) in the GSE65904 data set and found a similar result for ETV7 (p = 0.038, Figure 3J). By contrast, the remaining genes showed no prognostic significance ( Supplementary  Figures 2A-H). In addition, data from the GSE19234 data set showed that melanoma patients with low ETV7 had significantly shorter survival time since metastasis (p = 0.045, Figure 3K).
We then investigated the prognostic relevance of ETV7 in different subgroups of melanoma patients. As shown in Figures  4A, B, low ETV7 expression was correlated with a significantly short OS in melanoma patients younger than 60 y ( Figure 4A

Coexpression and Enrichment Analysis
The aforementioned analysis indicated that ETV7 was significantly downregulated in melanoma and was associated with the prognosis of melanoma patients. To understand the function of ETV7 in melanoma, we first identified the ETV7related genes, which were defined as those exhibiting strong correlation (coefficient >= 0.5, p < 0.05) with ETV7 at the transcriptional level in melanoma. A total of 623 ETV7-related genes were found in the TCGA_SKCM data set (Supplementary Table 1). GO analysis indicated that these genes were enriched in immune-related biological processes and molecular functions such as T cell activation, regulation of T cell activation, regulation of lymphocyte activation, response to interferon-gamma, regulation of leukocyte proliferation, MHC protein binding, cytokine receptor activity, chemokine receptor binding, MHC class II receptor activity, and CCR chemokine receptor binding (Supplementary Figure 3A). KEGG analysis suggested that these genes were enriched in antigen processing and presentation, cell adhesion molecules, cytokine-cytokine receptor interaction, Th1 and Th2 cell differentiation, Th17 cell differentiation, chemokine signaling pathway, natural killer (NK) cell-mediated cytotoxicity, and PD-L1 expression and PD-1 checkpoint pathway in cancer ( Figure 5A, Supplementary Table 2). We further searched the cBioportal database and found that in addition to the TCGA_SKCM data set, four of nine melanoma-related data sets contained gene expression data: MSKCC NEJM, DFCI Science, UCLA Cell, and TGEN Genome Res. Coexpression analysis of ETV7 in these four data sets identified 1 021 ETV7-related genes in the MSKCC NEJM cohort, 981 genes in the DFCI Science cohort, 440 genes in the UCLA Cell cohort, and 667 genes in the TGEN Genome Res cohort. The GO and KEGG analyses of these genes in each cohort consistently indicated that they were enriched in similar immune-related processes or pathways (Figures 5B-E, Supplementary Figures 3B-E). The coexpression and subsequent GO and KEGG analyses of GSE19234 further verified the association between ETV7 and immune-related pathways ( Figure 5F, Supplementary Figure 3F).
In addition, GSEA was applied to investigate the difference in signaling pathways between the low and high ETV7 expression subgroups in the TCGA_SKCM and GSE65904 data sets. The results inferred that various immune-related pathways were enriched in the high ETV7 expression subgroups. These pathways included antigen processing and presentation, Tolllike receptor signaling pathway, T cell receptor signaling pathway, RIG-I like receptor signaling pathway, JAK-STAT signaling pathway, B cell receptor signaling pathway, NK cellmediated cytotoxicity, cytokine-cytokine receptor interaction, and Fc gamma R-mediated phagocytosis (Figures 6A-I, Supplementary Table 3).

Immune Profile
The aforementioned analysis suggests a close link between ETV7 and immune-related pathways. We thus employed the ESTIMATE method to evaluate the overall TME status in the high-and low-ETV7 subgroups of melanoma patients. In the TCGA_SKCM and GSE65904 cohorts, the immune and stromal scores were significantly higher in the high-ETV7 subgroups than in the low-ETV7 subgroups (Figures 7A, B). All phenotypic and functional markers of T cells-CD3E, CD4, CD8B, FOXP3, GZMB, PRF1 and TBX21-were expressed at higher levels in melanoma patients with high ETV7 expression ( Figure 7C). Notably, the inhibitory immune receptors or ligands (CTLA4, LAG3, and PDCD1) and activating immune receptors (CD27, CD40, CD80, ICOS, TNFRSF4, and TNFRSF9) were significantly elevated in patients with high ETV7 expression, indicating a complex immune response in this subgroup ( Figures 7D, E). Among the immune modulators, ENTPD1 was highly expressed in the high-ETV7 subgroup of melanoma patients, whereas NT5E showed no significant difference ( Figure 8F). Moreover, IFNg signatures (CXCL10, CXCL9, IDO1, IFNG, and STAT1) and myeloid lineage phenotypic and functional markers (CD14, CD163, CD33, and CD68) were also significantly increased in melanoma patients with high ETV7 (Figures 7G, H). However, ARG1, a previously reported M2 macrophage marker, was not differently expressed in the two subgroups ( Figure 7H).
Further, multiple immune deconvolution methods-namely, TIMER, QUANTISEQ, CIBERSORT, MCPCOUNTER, XCELL, and EPIC-were manipulated to estimate the abundances of immune cells infiltrating into melanoma (34,35). Analyzing the relationship between ETV7 and the abundance of different immune cells, we observed a strong positive correlation between ETV7 and CD8+ T cells, as determined by TIM ( Figure 8A Previous studies indicate that ETV7 is predominantly expressed in hematopoietic tissues (36,37); thus, the strong positive correlation between ETV7 and CD8+ T cell infiltration To test this hypothesis, we analyzed the GSE72056 data set, a public single-cell RNA sequencing (scRNA-seq) data set of melanoma patients (30). In this data set, the tumor cells were designated by copy number variation analysis, whereas the nonmalignant cells were annotated as six types of cells-T cells, B cells, macrophages, endothelial cells, cancer-associated fibroblasts, and NK cells-based on preferentially or uniquely expressed marker genes (30). The remaining unresolved cells were referred to as undefined cells in the current study. Consistent with the original study by Itay Tirosh et al. (30). non-linear dimensionality reduction (t-distributed stochastic neighbor embedding (t-SNE)) analysis revealed that these eight clusters of cells could be distinguished separately ( Figure 9A). Tumor cells exhibited markedly elevated expression of MLANA, a widely used biomarker of melanoma ( Figure 9B).  Figure  4B) (30). As shown in Figure 9D, ETV7 was not exclusively expressed in T cells but could be detected in several types of cells, including tumor cells (green circle), T cells, and macrophages. ETV7 was also not detected in most CT8+ T cells (labeled as CD8A+ and/or CD8B+, Figures 9D-F), and a negligible fraction of ETV7-positive T cells ( Figures 9D-F, blue circle) were not recognized as CD8+ T cells. Collectively, these results indicate that the strong positive correlation between ETV7 and infiltration of CD8+ T cells cannot be attributed to ETV7 expression in these cells. In addition, we found that skin (green block, Figure 9G) showed the highest ETV7 expression across a set of normal human cells or tissues, including T cells (denoted by a blue block, Figure 9G, Supplementary Figure 4C).

Construction of an ETV7-Related Prognostic Model
As shown in Figure 10A, 98 genes are commonly shared among the ETV7-related genes in the five aforementioned data sets. A PPI network was constructed in the STRING database and visualized in Cytoscape ( Figure 10B). Consistent with the results of GO, KEGG, and GSEA analyses, ClueGo analysis demonstrated that these 98 genes were predominantly involved in immune-related biological processes such as T cell activation, chemokine-mediated signaling pathway, T cell differentiation, and response to interferon-gamma ( Figures 10C, D), as well as in immune-related pathways such as the T cell receptor signaling pathway, cytokine-cytokine receptor interaction, and TNF signaling pathway ( Figures 10E, F) (38). The 10 hub genes  were identified as CCR5, IFNG, TBX21, CXCL10, PRF1, CD2, CXCR3, CXCL9, CCL5, IL15 ( Figure 10B). All melanoma patients exhibiting low expression of these hub genes showed a significantly shorter OS ( Supplementary Figures 5A-J).
To construct an ETV7-related prognostic model for the disease-specific survival of melanoma patients, ETV7 and 10 hub genes were input into a LASSO Cox regression model in the GSE65904 data set (n = 210). Six genes were selected using the partial likelihood deviance method (Figures 11A, B). These genes were CCR5, IFNG, TBX21, CXCL10, CXCR3, and CCL5. The risk score was calculated by inputting the selected signature genes into the aforementioned formula. The median value of the risk score was set as the cutoff value, dividing patients into lowand high-risk subgroups. Prognostic analysis with the Kaplan-Meier method showed that melanoma patients with low risk scores had a significantly longer DSS ( Figure 11C, p = 0.0047) in the GSE65904 cohort. Time-dependent ROC analysis in this cohort demonstrated that the risk score had a favorable predictive value ( Figure 11D, AUC at 1 year = 0.64, AUC at 3 y = 0.69). The prognostic value of the established risk score was validated in the TCGA cohort ( Figure 11E, p < 0.0001). The ROC analysis in this data set showed that the AUC at 1 y reached 0.7 and at 3 y was 0.65 ( Figure 11F). Further, melanoma patients with high risk scores showed significantly shorter survival times since metastasis than those with low risk scores in the GSE19234 cohort (p = 0.021, Figure 11G). The ROC analysis in this data set indicated that the AUC at 1 y reached 0.69 and at 3 y was 0.98 ( Figure 11H).

DISCUSSION
In this systematic analysis of the ETS family in melanoma, we first integrated the data from the TCGA and GTEx databases and three available cohorts of melanoma patients in the Oncomine database; we then found that ELK3, ETS1, ETV1, and ETV5 were significantly upregulated in melanoma, whereas ELF3, ELF5, ETS2, ETV7, and SPDEF were significantly downregulated in the tumor. Prognostic analysis of these DEGs in the TCGA_SKCM and GSE65904 data sets consistently showed that only ETV7 had significant prognostic relevance in melanoma ( Figures 3A, J). Moreover, the gene could serve as an independent prognostic predictor after correcting for other confounding factors (Table 1). Indeed, all patients in various subgroups of melanoma with low ETV7 expression showed a significantly shorter OS than those with high ETV7 expression ( Figure 4). Several studies reported that ETS1 supported cell growth (39), prevented apoptosis (40), and facilitated the invasion of melanoma (22); however, an early study demonstrated that the gene was widely expressed in benign and malignant melanocytes, and immunohistochemical analysis indicated that its expression had no significant association with the clinical outcome (41), supporting the results of the current study. Two studies also observed that  (42,43). Although ETV1 showed no significant prognostic relevance in melanoma in the present study, its importance in a subset of melanomas with ETV1 amplification could not be underestimated and ignored. In addition, no studies have reported on the expression and function of ELK3, ETV5 ELF3, ELF5, ETS2, and SPDEF in melanoma. Although in this study, their expression was not correlated with prognosis, their involvement in the development of certain cases of melanoma could not be excluded, and future research has to be conducted. As a member of the ETS transcription factor family, ETV7 in melanoma has not been previously investigated. Our knowledge of the gene under physiological and pathological conditions is considerably limited. Initially identified at the beginning of this century, ETV7 was found to be highly related to the oncogenic ETV6 (36). However, recent studies have found that the two genes exerted opposite biological effects, and ETV7 can inhibit the transcriptional activities of ETV6 (44,45). ETV7 has earlier been reported to act as a hematopoietic oncoprotein; however, three recent studies found that the gene could suppress the proliferation, migration, and invasion of some solid tumors such as oral squamous cell carcinoma and nasopharyngeal carcinoma (45)(46)(47)(48)(49). To explore the role of ETV7 in melanoma, we analyzed seven independent melanoma-related data sets (TCGA_SKCM, GSE65904, GSE19234, and four data sets from the cBioportal database) by conducting GSEA, GO, and KEGG analyses, which are widely used bioinformatics tools in the functional characterization of specific genes (50)(51)(52)(53). All of these enrichment analyses suggest broad ETV7 involvement in immune-related processes and pathways, including the following: T cell activation; Th1, Th2, and Th17 cell differentiation; PD-L1 expression and PD-1 checkpoint pathway in cancer; and antigen processing and presentation ( Figures 5-7). Melanoma patients with high ETV7 expression had significantly higher immune scores than those with low ETV7, suggesting higher immune cell infiltration in the tumor microenvironment. Indeed, when the infiltration of various immune cells in melanoma was estimated, ETV7 showed a strong positive correlation with the infiltration of CD8+ T cells ( Figure 9). Moreover, melanoma patients with high ETV7 exhibited a significantly higher expression of various T cell markers, such as CD3E, CD4, CD8B, and TBX21 (54); this finding verified the close link between ETV7 and T cell infiltration. With the ETV7 expression in hematopoietic cells considered, the strong correlation between ETV7 and the infiltration of CD8+ T cells may be attributed to ETV7 expression in T cells rather than tumor cells. Thus, we first evaluated ETV7 expression across a set of normal human cells or tissues in the human protein atlas; the skin was found to exhibit the highest ETV7 expression. We then analyzed the public scRNA-seq data set (GSE72056) and found that most of the  CD8+ T cells showed no ETV7 expression, and some ETV7positive T cells were not CD8+ T cells (Figures 9C-F). ETV7 could be detected in some melanoma tumor cells ( Figure 9D). ETV was highly expressed in normal skin and showed markedly decreased expression in melanoma ( Figures 2L, N and 9G); thus, not all melanoma cells might be expected to express ETV7 ( Figure 9D).
Moreover, early studies inferred a positive regulatory role of ETV7 in T and B cells (45,48) promote B cell proliferation by targeting the cell cycle and apoptotic regulators (48). The results obtained from the present study, combined with other studies, highly suggest that ETV7 can have an essential role in T cell differentiation, proliferation, infiltration, and activation in melanoma.
Notably, several members of the ETS family regulate the immune system (55,56). ETS1 regulates the differentiation of several types of immune cells, such as T helper cell subsets and cytotoxic T cells; it also directly controls the expression of cytokine and chemokine genes (57). The mechanism by which ETV7, a member of the ETS family, regulates the immune process and response in melanoma has yet to be determined. However, 10 hub genes were identified in the current study ( Figure 10B). Among them, CD2 is a transmembrane glycoprotein typically known for its participation in the costimulatory pathway of T cell activation (58). TBX21 is essential for naive T lymphocyte development and interferongamma production (59). IL15 is identified as a cytokine that stimulates the proliferation of T lymphocytes (60). PRF1 encodes the protein perforin, which is present in T and NK cells and facilitates the release of granzymes and subsequent cytolysis of target cells (61). CCR5, CXCL10, CXCR3, CXCL9, CCL5, and IFNG also participate in the migration, trafficking, or differentiation of T cells (61)(62)(63)(64)(65). Whether ETV7 binds the promoter of these genes and stimulates their expression remains unknown, but our analysis showed a high correlation between them and ETV7 at the transcriptional level. To clarify this concern, future studies have to be conducted.
We further constructed a prognostic model by using LASSO Cox regression, a broadly selected machine learning algorithm to minimize the risk of overfitting (66). The model had a relatively good predictive value in training and testing cohorts ( Figures  11D, F). All six gene signatures were immune-related and regulated the proliferation, differentiation, migration, or activation of T cells (59,60,62,63,65). In addition, melanoma patients with low expression of these six gene signatures had a significantly shorter OS (Supplementary Figure 5). T cells, particularly CD8+ T cells, demonstrate a predominant role in the anticancer effects of ICIs. Thus, the model can potentially help select patients who can be responsive to ICIs; however, this hypothesis should be tested by prospective analysis in multicenter cohorts.

CONCLUSIONS
A systematic analysis of the ETS family in melanoma identified an essential role of ETV7 in regulating the immune microenvironment of the disease. Further studies are recommended to explore the exact mechanism that allows ETV7 to regulate the proliferation, migration, infiltration, and activation of immune cells, particularly CD8+ T cells.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
HQ and HZ contributed equally to this work. The study was conceived and designed by ZS. Analysis of data was performed by HQ, HZ, and ZS. Bioinformatics analysis in R was conducted by ZS and XZ. Useful discussion and advice were provided by YL, FL, and LS. ZS drafted the article. All authors contributed to the article and approved the submitted version.