Low expression of endoplasmic reticulum stress-related gene SERP1 is associated with poor prognosis and immune infiltration in skin cutaneous melanoma

SERP1 expression in pan-cancers and SKCM patients

The association between SERP1 expression and clinical characteristics in SKCM patients were list in Table 1. By analyzing The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) database, we obtained the expression of SERP1 RNA in pan-cancer. As Figure 1A revealed, compared to normal tissues, there was significantly different expression of SERP1 mRNA in 33 incorporated cancers except KICH, KIRC and THCA. Because there were only tumor-related samples without normal tissue samples (MESO and UVM) or too few normal tissue samples (SARC having 2 samples) in the database for these three cancer types, these three cancers could not be compared to normal tissue in the pan-cancer comparison and therefore no results are shown. As the target of our study, SERP1 expression was lower in SKCM tumors than in normal tissues (p = 0.002, Figure 1B). Then SKCM patients were divided into the high SERP1 expression group and low SERP1 expression group based on the median SERP1 expression. We compared the RNA expression between these two groups. There were 111 RNAs that met the established selected threshold and were recognized as differentially expressed genes (DEGs) (Adjust P-value < 0.05 and absolute log-fold change > 3) (Figure 1C), of which 30 were upregulated (logFC is positive) and 81 were downregulated (logFC is negative). Top 15 up-regulated and down-regulated DEGs were illustrated by heatmaps and included in the table (Figure 1D, 1E and Supplementary Table 2). We also analyzed RNAs correlated with SERP1 expression in SKCM patients and showed the top 50 positively and negatively associated RNAs as heatmaps (Figure 2A, 2B and Supplementary Table 3).

Protein-protein interaction (PPI) and gene–gene interaction (GGI) network analyses of SERP1 in SKCM patients

The PPI network analysis was conducted to explore the potential interactions of SERP1 protein. As Figure 2C presents, the network with interaction nodes and edges was built via STRING, whose top 10 proteins are listed in Supplementary Table 4. The GGI network also showed that the functions of the differentially expressed SERP1 and its associated genes (such as SERP2, XBP1, DNAJB9, SEC23B, ARF4, SLC33A1, BET1, COPB1, SSR1, SRPRA, TMED10, SSR2, PGM3, DNAJB1, HSPA5, C6orf120, TGDS, MGAT2, GOLGA5, TSC22D2) were primarily related to response to topologically incorrect protein, cellular response to unfolded protein, response to unfolded protein, cellular response to topologically incorrect protein, endoplasmic reticulum unfolded protein response, response to endoplasmic reticulum stress and regulation of endoplasmic reticulum unfolded protein response (Figure 2D).

Predicted functions and pathways of the DEGs between high- and low-SERP1 expression in SKCM patients

The functions of the DEGs between high- and low- SERP1 expression groups were predicted by analyzing Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) via R software and Cytoscape software. The different RNA functional of DEGs included three categories: the biological process (BP), molecular function (MF), and cellular component (CC). The top three GO terms of BP, CC and MF functional groups found via R software are present as Figure 3A.

The BP ontology included keratinization, epidermis development and keratinocyte differentiation. In the CC ontology, most genes were connected with cornified envelope, intermediate filament and keratin filament. In addition, the structural constituent of epidermis, structural constituent of cytoskeleton with serine-type endopeptidase inhibitor activity made up the majority of MF ontology. Furthermore, the KEGG pathway analysis showed that the DEGs were most closely correlated to Estrogen signaling pathway, GnRH signaling pathway and Inflammatory mediator regulation of TRP channels. Detailed results of GO and KEGG are placed in Supplementary Table 5.

On the other hand, we also performed GO and KEGG analysis of DEGs using ClueGO and the results were presets as Figure 3B. Results followed cut-off threshold P-value < 0.01 were just enriched in BP, CC and KEGG (Figure 3C). The BP ontology included keratinization (20%), antimicrobial humoral response (18%), regulation of epidermis development (12%), establishment of skin barrier (8%), cornification (6%) and hemidesmosome assembly (2%). The CC ontology included integrator complex (10%), intermediate filament (16%), and desmosome (2%). The KEGG pathway was Inflammatory mediator regulation of TRP channels (6%). The top 10 nodes ranked by MCC of cytoHubba were SPRR1B, CDSN, RPTN, IVL, SPRR2G, LOR, SPRR2E, EVPL, PI3 and TGM1 (Figure 3D) and modules with MCODE score = 30 including LCE1C, SPRR2A, LCE3E, LCE3D, LCE2D, LCE2C, LCE1E, SPRR1A, SPRR2D, SPRR2B, SPRR2E, SPRR2G, LCE1B, LCE3A, LCE1A, LCE3C, LCE5A, RPTN, LCE1D, SPRR1B, PI3, LOR, IVL, CASP14, LCE6A, CDSN, TGM1, LCE2B, LCE1F, EVPL were present in Figure 3E.

Gene set enrichment analyses (GSEA) identifies DEGs between high- and low-SERP1 expression related signaling pathways

We additionally conducted GSEA to identify signaling pathways that were differentially activated between high- and low-SERP1 expression groups in SKCM. The results indicated the Top 15 pathways were associated with Reactome GPCR ligand binding, Reactome G alpha i signaling events, Reactome class a1 rhodopsin like receptors, Reactome leishmania infection, Naba core matrisome, Kegg neuroactive ligand receptor interaction, Kegg cytokine-cytokine receptor interaction, Wp GPCRs class a rhodopsin-like, Reactome anti-inflammatory response favouring leishmania parasite infection, Reactome G alpha q signaling events, Kegg chemokine signaling pathway, Reactome peptide ligand binding receptors, Reactome cell surface interactions at the vascular wall, Reactome FC epsilon receptor (FCERI) signaling and Reactome immunoregulatory interactions between a lymphoid and a nonlymphoid cell (Figure 4A–4O). We further list the Top 50 pathways in Supplementary Table 6.

Genetic alteration and protein localization of SERP1 in SKCM patients

We analyzed the genetic mutations of SERP1 expression in SKCM patients via the cBioPortal online tool. Based on TCGA, SERP1 mutated lowly in SKCM via pan-cancers analysis (Figure 5A). Figure 5B showed the mutation rate of SERP1 genes was 1.1%. Figure 5C indicated the major form alterations of the SERP1 genes in SKCM is amplification. Figure 5D presented there was an overall higher amount of amplification occurrence of SERP1 mRNA expression in SKCM patients. We further explored the correlation between genetic mutations and the prognosis of SKCM patients. The statistically significant result was showed in Figure 5E that SERP1 genetic mutations decrease overall survival (OS) of SKCM patients (p = 0.0402). However, it is noteworthy that the expression of SERP1 protein did not show significant differences between tumor tissues and normal tissues in the Human Protein Atlas (Figure 5F).

Correlations between SERP1 and prognosis in SKCM patients

In order to explore the prognosis value of SERP1 in SKCM patients, we performed Kaplan-Meier analysis to evaluate its value on prediction of SKCM on clinical outcomes. The result found that OS (hazard ratio (HR): 0.62, 95% confidence interval (CI): 0.48–0.82, p = 0.001 (Figure 6A), disease specific survival (DSS) (HR: 0.65, 95% CI 0.49-0.87, p = 0.004) (Figure 6B) and progress free interval (PFI) (HR: 0.77, 95% CI 0.61-0.96, p = 0.023) (Figure 6C) for high SERP1 groups were all statistically better than those for the low SERP1 groups. In addition, in order to comprehensive understand the multidimensional prospect of the correlation of SERP1 expression to SKCM patient’s survival, subgroup Kaplan-Meier analysis of OS was also performed according to different clinical variables. As Figure 6D–6O revealed, low SERP1 expression was significantly associated with worse OS in SKCM patients both of gender, female (p = 0.009) or male (p = 0.006), race white (p = 0.001), age ≤ 60 (p = 0.005), T stage is T2-4 (p = 0.001), N stage is N1-3 (p = 0.008), M stage is 0 (p = 0.006), Pathologic Stage is II-IV (p = 0.007), without radiation therapy (p = 0.001), Tumor tissue site is on the trunk (p = 0.003), having melanoma ulceration (p = 0.003) and Melanoma Clark Level Stage is II-V (p = 0.012). However, the high or low SERP1 expression on other cases did not show statistical significance on OS (Supplementary Figure 1A–1L).

Univariate Cox analysis showed that SERP1 expression was an independent risk factor for OS (HR: 0.591, 95% CI 0.405-0.861, p = 0.006) (Figure 7A and Supplementary Table 7) and DSS (HR: 0.584, 95% CI 0.395-0.864, p = 0.007) (Figure 7B and Supplementary Table 8). However, SERP1 expression did not have significant predictive power for PFI (HR: 0.769, 95% CI 0.555-1.067, p = 0.116) (Supplementary Table 9). Also available as independent risk factors through Univariate Cox analysis were N stage (HR: 2.945, 95% CI 1.887-4.597, p < 0.001), Melanoma ulceration (HR: 1.645, 95% CI 1.062-2.546, p = 0.026) and Breslow depth (HR: 2.028, 95% CI 1.312-3.136, p = 0.001) for OS and N stage (HR: 2.893, 95% CI 1.813-4.616, p < 0.001), Melanoma ulceration (HR: 1.709, 95% CI 1.087-2.687, p = 0.020) and Breslow depth (HR: 1.719, 95% CI 1.096-2.695, p = 0.018) for DSS, respectively.

The nomograms predicting OS (Figure 8A) and DSS (Figure 8B) was constructed based on significant risk factors identified from in the univariate Cox analyses. The score of the corresponding risk factors should be first identified by the point scale at the top of the nomogram. All scores were then summed and the corresponding 1-year, 3-year, and 5-year survival probability was obtained corresponding to the scales at the bottom of the nomogram. The time-dependent ROC curves of OS and DSS are shown in Figure 8C, 8D, respectively. The 1, 3, and 5-year AUCs of OS are 0.404 (95% CI 0.308-0.500), 0.388 (95% CI 0.324-0.453) and 0.398 (95% CI 0.334-0.463), respectively. The 1, 3, and 5-year AUCs of DSS are 0.432 (95% CI 0.324-0.539), 0.401 (95% CI 0.333-0.469) and 0.408 (95% CI 0.341-0.475), respectively. These suggesting that low SERP1 has a good predictive efficacy and diagnostic accuracy for survival of SKCM patients. Details of the ROC information are presented in Supplementary Table 10.

Calibration curves showed a strong consistency between the possibilities obtained by the nomogram and the real results of 1-year and good consistency of 3 and 5-year for OS (Figure 8E) and DSS (Figure 8F). Furthermore, the C-index values were 0.686 (95% CI 0.662 to 0.710) and 0.664 (95% CI 0.638 to 0.690) for OS and DSS, respectively. These revealed the good credibility of the nomogram.

Correlation between SERP1 expression and clinical variables in SKCM patients

The Figure 9A–9E revealed the correlation between the expression of SERP1 and the different clinical variables of SKCM patients. SERP1 expression level significantly related to the T stage (p = 0.037), Pathologic stage (p = 0.009), Radiation therapy (p = 0.013), Breslow depth (p < 0.001) and Melanoma ulceration (p = 0.043). The higher SERP1 expression was associated with lower T stage and Pathologic stage, shallower Breslow depth and fewer ulceration. Moreover, the expression of SERP1 was lower in patients with radiation therapy experience.

Correlation of SERP1 and immune cell infiltration in SKCM

The results found that T helper cells (r = 0.510, p < 0.001), Tcm (r = 0.340, p < 0.001), Tgd (r = 0.330, p < 0.001), Th2 cells (r = 0.230, p < 0.001), Macrophages (r = 0.200, p < 0.001), Th1 cells (r = 0.190, p < 0.001), B cells (r = 0.180, p < 0.001), aDC (r = 0.140, p = 0.002), T cells (r = 0.140, p = 0.003), CD8 T cells (r = 0.130, p = 0.004), Eosinophils (r = 0.130, p = 0.005) showed positive association with SERP1 expression. The NK CD56 bright cells (r = -0.210, p < 0.001), NK cells (r = -0.180, p < 0.001), Mast cells (r = -0.170, p < 0.001), pDC (r = -0.130, p = 0.005), Th17 (r = -0.097, p = 0.036) were negatively correlated with SERP1 (Figures 10A, 11). However, the TFH, Cytotoxic cells, iDC, NK CD56dim cells, Neutrophils, Tem, DC, Treg showed no significant correlation with SERP1 (Supplementary Figure 2).

We further evaluated the infiltration levels of 24 immune cells above in high- or low-SERP1 expression group (Figure 10B). The result indicated in high SERP1 expression group, the T cells (p = 0.018), aDC (p = 0.009), B cells (p < 0.001), Eosinophils (p = 0.027), Macrophages (p < 0.001), T helper cells (p < 0.001), Tcm (p < 0.001), Tgd (p < 0.001), Th1 cell (p < 0.001), Th2 cells (p < 0.001) infiltration levels were significantly higher than those in low SERP1 expression group. But the NK CD56 bright cells (p < 0.001), NK cells (p = 0.016), Mast cells (p = 0.003) infiltration level was significant lower in high SERP1 expression group than low SERP1 expression group.

Moreover, the association of SERP1 expression with markers of mainly relative immune infiltration cells and immune checkpoints of SKCM were explored. The former mainly included.

CD8+ T cell, monocyte, tumor-associated macrophage (TAM), M1 macrophage, M2 macrophage, neutrophils, DC, Th1, Th2, Tfh, Th17 and Treg. The latter mainly included PD-1, PD-L1, CTLA4, LAG3, TIM-3, GZMB, TIGIT and BTLA. We found that SERP1 expression is significantly positive correlated with the most (36 of 41) immune cell markers in SKCM (P < 0.05) (Table 2). Although none of the correlations were high (r < 0.5), but it suggests that SERP1 is extensively involved in regulating immune cell function in SKCM. It is worth noting that, SERP1 expression was positively correlated with numerous immune checkpoint gene markers, including PD-1, PD-L1, PD-L2, CTLA4, LAG3, TIM-3, GZMB, TIGIT (Table 3). These results suggest that SERP1 expression may be associated with the immunotherapeutic effect of SKCM.

Gene expression and clinical data

We obtained the mRNA expression and clinical data from TCGA (https://cancergenome.nih.gov). The data for the corresponding 813 normal tissue samples were obtained from the GTEx (https://gtexportal.org/). Samples with “0” values for gene expression or with inadequate prognostic information were excluded.

RNA-sequencing data in Fragments Per Kilobase per Million format (FPKM) were converted and normalized by the Toil [51] process as transcripts per million (PTM) reads and log2 transformed for further analysis. The clinical characteristics of the 472 patients included in the study are summarized in Supplementary Table 1. Location and qualitative data of SERP1 protein in SKCM and normal tissues were evaluated by the Human Protein Atlas (https://www.proteinatlas.org/) [52].

DEGs analysis between SKCM patients with high and low SERP1 expression

We performed a “DESeq2” analysis [53] in R to identify DEGs between SERP1-high and SERP1-low SKCM patients identified by unpaired Student’s t-test. Thresholds were set as an adjusted P < 0.05 and absolute log-fold change > 3. Identified genes were analyzed and presented as volcano plots. The top 15 up- or down-regulated genes were presented as heat maps. All data visualization was achieved using the “ggplot2” package in R. We set a relatively high threshold (log-fold change > 3) with the aim of selecting the DEGs with the greatest degree of change associated with SERP1-high and SERP1-low expression. The top 15 DEGs were selected to further screen the genes most associated with SERP1-high and SERP1-low expression.

SERP1 correlation analysis in patients with SKCM

We performed correlation analysis between SERP1 and other RNAs in patients with SKCM using TCGA data via the “stat” package in R. Pearson correlation coefficients were calculated and the top 50 genes most positively and negatively associated with SERP1 were filtered to construct heat maps using the “ggplot2” package in R.

PPI network and GGI network analysis

To collect and integrate potential protein interactions with SERP1, we searched the STRING database (https://string-db.org/) [54] and conducted a PPI network analysis. A confidence score > 0.7 was set as the significance threshold. GeneMANIA (http://www.genemania.org) [55] uses extensive genomics and proteomics data to discover functionally similar genes. The database is used to generate hypotheses about gene function, analyze gene lists, and prioritize genes for functional analysis. We searched for SERP1 in this database to predict the gene–gene interaction (GGI) network.

Functional enrichment analysis of DEGs between patients with SKCM with high- and low-SERP1 expression

The “clusterProfiler” [56] and “org.Hs.eg.db” packages of R and the “ClueGO” app of Cytoscape (v3.8.2) were used to conduct GO function and KEGG pathway enrichment analyses for statistically significant DEGs. A p-value < 0.01 was set as the cut-off threshold for GO and KEGG pathway enrichment analyses. The results were presented as a bar plot via the “ggplot2” package in R. Additionally, we also imported the DEGs into the STRING database to build a PPI network map and completed the GO and KEGG signaling pathway enrichment analyses using the ClueGO plugin of Cytoscape. The MCODE and cytoHubba plugins were used to identify key modules. The top 10 nodes ranked by MCC of cytoHubba and modules with MCODE score = 30 were presented.

Gene set enrichment analysis

We performed GSEA via the “clusterProfiler” [56] package to determine the biological pathway differences between high- and low-SERP1 groups. Pathways with a false discovery rate (FDR) < 0.25 and an adjusted p-value < 0.05 were considered to be remarkably changed. Gene set permutation was performed 1,000 times for each analysis.

cBioPortal

We searched and downloaded mutation information and corresponding clinical data from cBioPortal (https://www.cbioportal.org/) [57], a comprehensive database that provides visual and multidimensional cancer genomics data.

Kaplan–Meier analysis

We conducted Kaplan–Meier analysis via the “survival” package to compare the OS, DSS, and PFI rates between the high- and low-SERP1 gene expression groups. Subgroup Kaplan–Meier analysis of OS was also performed based on different clinical variables. The p-value was determined by Cox regression analysis and the survival curves were visualized via the “survminer” package.

Correlation between SERP1 expression and immune infiltration in SKCM

We used the gene set variation analysis package [58] to investigate the correlation between the SERP1 expression and tumor-infiltrating immune cells in patients with SKCM. Twenty-four types of immune cells were included in the analysis: T cells, aDC (activated DC), B cells, CD8 T cells, cytotoxic cells, DC, eosinophils, iDC (immature DC), macrophages, mast cells, neutrophils, NK CD56 bright cells, NK CD56 dim cells, NK cells, pDC (plasmacytoid DC), T helper cells, Tcm (central memory T cells), Tem (effector memory T cells), Tfh (T follicular helper cells), Tgd (gamma delta T cells), Th1 cells, Th17 cells, Th2 cells, and Treg (regulatory T cells). The association of SERP1 with immune cell markers as well as immune checkpoints were also explored. Spearman’s correlation was used to evaluate the correlation of gene expression with p-values < 0.05 considered statistically significant.

Statistical analysis

R software was used to perform statistical analyses in this study (version 3.6.3). The “ggplot2” package was used to present SERP1 gene expression as dot plots in patients with pan-cancer and SKCM. The median method of gene expression was selected for cutoff values. SERP1 expression in patients with SKCM with different clinical characteristics was analyzed using the Chi-squared test, Fisher test, and Wilcoxon rank sum test depending on the situation. The “survival” package was used to analyze the effect of SERP1 on survival with other clinical characteristics in SKCM patients via univariate and multivariate Cox regression. The predictive nomogram of 1-, 3-, and 5-year OS, DSS and PFI with SERP1 expression for SKCM patients and other relevant clinical parameters was constructed via a stepwise Cox regression model through the “rms” and “survival” packages. The ROC curve, C-index and calibration curve were used to verify the reliability of nomograms. ROC curves were analyzed and the areas under the ROC curve (AUC) were calculated via the timeROC (version 0.4) package of R software. The calibration curves were analyzed and plotted via rms (version 6.2-0) and survival (version 3.2-10) package of R software. The number of samples per group for our calibration curves is set to 40, the number of repetitions is 200, and the boot method is used.