Identification of a novel tumour microenvironment‐based prognostic biomarker in skin cutaneous melanoma

Abstract Skin cutaneous melanoma (SKCM) is one of the most destructive skin malignancies and has attracted worldwide attention. However, there is a lack of prognostic biomarkers, especially tumour microenvironment (TME)‐based prognostic biomarkers. Therefore, there is an urgent need to investigate the TME in SKCM, as well as to identify efficient biomarkers for the diagnosis and treatment of SKCM patients. A comprehensive analysis was performed using SKCM samples from The Cancer Genome Atlas and normal samples from Genotype‐Tissue Expression. TME scores were calculated using the ESTIMATE algorithm, and differential TME scores and differentially expressed prognostic genes were successively identified. We further identified more reliable prognostic genes via least absolute shrinkage and selection operator regression analysis and constructed a prognostic prediction model to predict overall survival. Receiver operating characteristic analysis was used to evaluate the diagnostic efficacy, and Cox regression analysis was applied to explore the relationship with clinicopathological characteristics. Finally, we identified a novel prognostic biomarker and conducted a functional enrichment analysis. After considering ESTIMATEScore and tumour purity as differential TME scores, we identified 34 differentially expressed prognostic genes. Using least absolute shrinkage and selection operator regression, we identified seven potential prognostic biomarkers (SLC13A5, RBM24, IGHV3OR16‐15, PRSS35, SLC7A10, IGHV1‐69D and IGHV2‐26). Combined with receiver operating characteristic and regression analyses, we determined PRSS35 as a novel TME‐based prognostic biomarker in SKCM, and functional analysis enriched immune‐related cells, functions and signalling pathways. Our study indicated that PRSS35 could act as a potential prognostic biomarker in SKCM by investigating the TME, so as to provide new ideas and insights for the clinical diagnosis and treatment of SKCM.


| BACKG ROU N D
Skin cutaneous melanoma (SKCM) is a malignant transformation of melanocytes derived from neural crest stem cells. 1 Although SKCM accounts for only approximately 5% of all skin tumours, it causes more than 75% of deaths from skin tumours. 2 Recently, the approval of targeted-and immune-based therapeutic agents with substantial activity has brought new hope for the treatment of SKCM. [3][4][5] However, approximately half of the patients treated with these therapies have no sustained response and display no satisfactory improvement in prognosis. 6,7 Moreover, some patients have immune-related adverse events. [8][9][10] Therefore, it is urgent and important to explore biomarkers that can identify the population that is more likely to benefit from these treatments.
Tumour cells constantly interact with the surrounding microenvironment. Increasing evidence has demonstrated that the tumour microenvironment (TME) plays a considerable role in the development of various tumours, 11 including SKCM, 12 and targeting the TME could complement traditional treatment and improve the therapeutic response and clinical outcome for this malignancy. 11,13,14 Infiltrating stromal and immune cells form the major fraction of normal cells in tumour tissues, and not only perturb the tumour signal in molecular studies but also play a vital role in cancer biology. 15,16 Previous studies have shown that melanoma cells deeply interact with the TME and the immune system 17,18 This knowledge has led to the identification of novel therapeutic targets and treatment strategies, including TME-based predictive and prognostic biomarkers. 18 Estimation of stromal and immune cells in malignant tumour tissues using the expression data (ESTIMATE) algorithm is an approach to study tumour-infiltrating immune cells and their interactions with cancer cells to infer the fraction of immune and stromal cells in tumour samples underlying the expression of gene signatures. 16,19 2 | MATERIAL S AND ME THODS

| Data processing
The public RNA-Seq data (fragments per kilobase of transcript per million fragments mapped) and mRNA expression profile data of SKCM from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) were retrospectively analysed, and relevant clinical data (sex, age, stage, tumour-node-metastasis (TNM) classification and survival data) were also obtained. A total of 471 cases of tumour tissue and one case of adjacent tissue were included. Moreover, 812 normal tissue samples were obtained from the Genotype-Tissue Expression (GTEx) database. 20,21 2.2 | TME scores generation ESTIMATE is a tool to predict tumour purity, as well as the presence of infiltrating immune and stromal cells in tumour tissues. 16 The ESTIMATE algorithm is based on single-sample gene set enrichment analysis (ssGSEA) and mainly generates three scores: ImmuneScore (which infers the infiltration of immune cells in tumour tissues), StromalScore (which represents the presence of stromal cells in tumour tissues), and ESTIMATEScore (which captures tumour purity). 16 Here, we generated StromalScore, ImmuneScore and ESTIMATEScore using the ESTIMATE package (version 1.0.13) to estimate the presence of immune and stromal cells in the TME for each tumour tissue. The higher the three TME scores, the larger the presence of the corresponding cells in TME. 22 Tumour purity was negatively correlated with the three scores. Moreover, we introduced tumour purity as a supplement.

| TME scores and clinicopathological characteristics
We generated four TME scores (ImmuneScore, StromalScore, ESTIMATEScore and tumour purity) and analysed the stage and TNM classification data of SKCM to determine the relationship between the ratio of stromal and immune components with clinicopathological characteristics.

| Identification of differentially expressed prognostic genes
We grouped all SKCM samples into high or low score groups according to the median score of TME scores, and then applied the survival package (version 3.2-7) and survminer package (version 0.4.8) to survival analysis. 23 The Kaplan-Meier method with the log rank test was used to plot the survival curve. We defined those that could group the population into different survival statistically as differential TME scores. After the identification of differential TME scores, we used the limma package (version 3.12) to screen differentially expressed genes in differential TME scores. We then used a Venn diagram to display the overlapping genes between differential TME scores and used these genes as differentially expressed prognostic genes.

| Univariate Cox regression analysis
We performed univariate Cox regression analysis on differentially expressed prognostic genes to assess their correlation with SKCM prognosis and used the forestplot package (version 1.10.1) for visualization.

| Development of a prognostic prediction model
Least absolute shrinkage and selection operator (LASSO) is a regression-based method that allows numerous covariates in the model and thus could regulate those covariates that might influence the overall regression. 24,25 Here, we further obtained more reliable prognostic genes by LASSO to prevent overfitting of highly correlated genes from further filtering. The Kaplan-Meier method with a log rank test was used to plot the survival curve of a single filtered gene and their combination. We used the survivalROC package (version 1.0.3) and pROC package (version 1.17.0.1) to draw receiver operating characteristic (ROC) curves and quantified the area under the ROC curves (AUCs) to evaluate the diagnostic performance of these genes alone and in combination. To verify the clinical significance of the differentially expressed prognostic gene-based risk model, we used univariate and multivariate Cox regression analyses to determine the clinical relevance of the model. We then explored the possibility of using these genes as prognostic factors to independently affect the survival prognosis of SKCM. Finally, a prognostic prediction model was constructed to predict the overall survival (OS).

| Identification of a novel prognostic biomarker and functional enrichment analysis
Based on previously established AUCs, we identified a novel prognostic biomarker. The limma package was used to evaluate the expression of the novel prognostic biomarker in SKCM tissues and normal tissues. Moreover, we also conducted clinicopathological characteristics analysis.
CIBERSORT is an algorithm that is widely used to characterize the cell composition of complex tissues through biomarker expression. 26 The LM22 signature, a special genetic marker that contains 547 genes, is usually used to distinguish 22 immune cell subtypes downloaded from CIBERSORT. In this study, the CIBERSORT package was used and the LM22 signature algorithm was used to calculate the infiltration

| TME scores and clinicopathological characteristics
We generated four TME scores (ImmuneScore, StromalScore, ESTIMATEScore and tumour purity) and analysed the stage and TNM classification data of SKCM to determine the relationship between different TME cells and clinicopathological characteristics.
Unfortunately, the results were not statistically significant ( Figure   S2), but the overall trends were relatively consistent.

| Survival analysis of TME scores
We analysed the correlation between the four TME scores and survival of SKCM patients. We found that there was no significant difference between StromalScore and SKCM OS (p = 0.084, Figure 1A) or ImmuneScore (p = 0.069, Figure 1B). In contrast, the ESTIMATEScore was positively correlated with OS (p = 0.037, Figure 1C), and tumour purity was negatively correlated with prognosis (p = 0.036, Figure 1D). Therefore, we considered the ESTIMATEScore and tumour purity as differential TME scores. Our results indicate that different TME cells might help predict the outcome of SKCM patients.

| Screening of differentially expressed prognostic genes
We identified the differentially expressed genes in differential TME scores through a comparison analysis of high-and lowdifferential TME score groups. There were 35 and 37 differentially expressed genes in ESTIMATEScore and tumour purity, respectively ( Figure 2A,B). Using the intersection, we identified 34 differentially expressed prognostic genes ( Figure 2C).

| Development of a prognostic prediction model
LASSO was used to generate the key differentially expressed prognostic genes. In LASSO, as logλ changed, the corresponding coefficients of certain differentially expressed prognostic genes were reduced to zero, indicating that their effects on the model could be omitted because they were shrinking parameters ( Figure 4A). Following crossvalidation, a total of seven differentially expressed prognostic genes (SLC13A5, RBM24, IGHV3OR16-15, PRSS35, SLC7A10, IGHV1-69D and IGHV2-26) achieved the minimum partial likelihood deviance ( Figure 4B) and thus were used to develop the risk model. The clinical relevance heat map of the risk model showed the influence of the model on the clinical variables related to SKCM ( Figure 4C).
As IGHV3OR16-15, IGHV1-69D and IGHV2-26 are pseudogenes, we only included SLC13A5, RBM24, PRSS35 and SLC7A10 in further analyses. The survival curves showed that whether the four differentially expressed prognostic genes were separated or combined together, they were able to effectively predict the outcome of SKCM patients ( Figure 4D, Figure S3). The combination of four differentially expressed prognostic genes had a favourable diagnostic performance (AUC = 0.582, Figure 4E). Furthermore, we calculated the AUCs of each differentially expressed prognostic biomarker and concluded that PRSS35 had the best diagnostic performance (AUC = 0.721, Figure 4F); therefore, PRSS35 was considered as a novel prognostic biomarker.
We compiled data from 346 SKCM patients with complete clinical information to perform univariate and Cox multivariate regression analyses. In the univariate analysis, several factors, such as age, stage, TNM classification and the risk score obtained by the constructed risk model, can affect the prognosis of SKCM (Table 1). In the multivariate Cox regression analysis, age, T classification, N classification and risk score could still affect the prognosis of SKCM (Table 1). This shows that the risk model we constructed can be utilized independent of other clinical traits and can aid in identifying prognostic biomarkers. Therefore, we determined the characteristics of the prognostic prediction model in SKCM patients ( Figure 5A) and conducted a nomogram to predict the 1-, 3-and 5-year survival of SKCM patients ( Figure 5B).

| Identification of a novel prognostic biomarker and functional enrichment analysis
After identification of PRSS35 as a novel prognostic gene, we validated that PRSS35 was highly expressed in SKCM tissues compared F I G U R E 1 Correlation of TME scores with the survival of SKCM. A, ImmuneScore. B, StromalScore. C, ESTIMATEScore. D, Tumour purity. SKCM, Skin cutaneous melanoma; TME, tumour microenvironment   231  150  100  71  52  41  27  226  149  114  7 9  5 7  4 9  4 2   229  148  100  7 4  5 6  4 5  3 3  228  151  114  76  53  45  36   228  147  99  72  53  42  29  229  152  115  78  56  48  40   230  152  115  78  56  48  40  227  147  99  72  53  42      to normal tissues (p < 0.0001, Figure 6A). Moreover, the exploration of clinical clinicopathological characteristics indicated that PRSS35 was associated with age (p = 0.0105) and T classification   drug tolerance or resistance, toxicity and the considerable expense of these methods. 7,28 Moreover, there is still a lack of robust prognostic biomarkers to guide clinical decision-making. Consequently, we attempted to identify a novel prognostic biomarker based on TME to contribute to SKCM therapy. In this study, we sought to develop a TME-related prognostic biomarker and to classify TNM in SKCM patients through a combined analysis of multiple databases. We identified that PRSS35 is involved in the TME of SKCM. More importantly, we further concluded that PRSS35 might reflect the TME status of SKCM and act as a novel TME-based prognostic biomarker in SKCM.

| DISCUSS ION
Tumour microenvironment has a considerable role in tumorigenesis and development. Transcriptome analysis of SKCM from TCGA showed that different immune components from the TME facilitate the prognosis of SKCM patients. Our results highlight the significance of exploring the interaction between SKCM cells and TME.
The National Comprehensive Cancer Network recommends antiprogrammed cell death protein 1 monotherapy (pembrolizumab, nivolumab) as first-line therapy for metastatic or unresectable SKCM and combination targeted therapy (dabrafenib/trametinib, vemurafenib/cobimetinib, encorafenib/binimetinib) as first-line therapy for the BRAF V600-activating mutation SKCM. 29 Despite the promising efficacy of immune checkpoint inhibitors, approximately half of the patients either still have no durable response or suffer from immune-related adverse events. 6-10 Therefore, the universality of these regimens and the identification of patients who are more likely to have survival benefits from these treatments are unknown, and it is necessary to explore TME-based prognostic biomarkers in SKCM.
In our study, we started with the transcriptome analysis of SKCM from TCGA and found that the ESTIMATEScore and tumour purity were correlated with survival. We then identified seven potential prognostic markers and narrowed it down to PRSS35, which showed the most promise as a novel prognostic biomarker.
Serine proteases play an important role in cancer progression and metastasis. 30,31 As a member of the serine protease family, PRSS35 may also contribute to the aetiology of several cancers. Our results suggest that PRSS35 is upregulated in SKCM tissues, which is consistent with the results obtained from studies on ovarian cancer. 32 Unfortunately, PRSS35 has not been reported to be immune-related, except in aortic stenosis and aortic insufficiency. 33 In other words, we report for the first time that PRSS35 upregulation is significantly associated with SKCM prognosis.
Furthermore, we uncovered the relationship between PRSS35 and TME. CIBERSORT indicated that immune-related cells, such as CD8 + T cells, neutrophils, activated mast cells and M0 macrophages. It is well known that these are immune-related cells. 34 The GSEA results suggested    There are some limitations to this study, although the established novel TME-based prognostic biomarker in SKCM is powerful. First, the transcriptome data used in the present study were obtained from TCGA and GTEx databases. Although we conducted homogenization, the source bias could not be ignored and might have affected the extrapolation of our results. Second, we were not able to obtain satisfactory results in the relationship between TME scores/PRSS35 and clinicopathological characteristics due to the limited clinical information and clinical sample size included in our study. Third, only one SKCM data set was applied in the present study; however, more data sets with richer clinical information and a larger sample size were collected to verify our findings. Finally, this was an in silico analysis; therefore, molecular biology experiments, as well as prospective, well-designed, multicentre studies are required to validate the findings.

B cells naive B cells memory Plasma cells T cells CD8 T cells CD4 naive T cells CD4 memory resting T cells CD4 memory activated T cells follicular helper T cells regulatory (Tregs) T cells gamma delta NK cells resting NK cells activated
In conclusion, we obtained TME scores using ESTIMATE and screened the differentially expressed prognostic genes after the identification of differential TME scores, and finally identified PRSS35 as a novel prognostic biomarker. LASSO, Cox regression F I G U R E 8 Functional enrichment analysis and functional enrichment analyses were also conducted using CIBERSORT and ssGSEA. By investigating the TME to obtain new insights for the clinical diagnosis and treatment of SKCM, our study determined that PRSS35 might act as a potential prognostic biomarker for this malignancy.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.