Differences in and verification of genetic alterations in chemotherapy and immunotherapy for metastatic melanoma

Gene expression data acquisition and screening

Data were collected according to the rules of public data use on official websites, and all normalized data were collected from The Cancer Genome Atlas data center (TCGA, https://portal.gdc.cancer.gov). Patients were screened based on the administration of chemotherapy (e.g., dacarbazine, cisplatin, etc.,) and immunotherapy (e.g., IL-2, interferon, etc.,). The exclusion criteria included samples with BRAF, NRAS, NF1, and KIT gene mutations and patients who were lost to follow-up.

Differential gene analysis and screening

After completing the grouping of samples (chemotherapy and immunotherapy), the DESeq2 package (version 1.24) of R software (version 4.02, R Foundation for Statistical Computing, Vienna, Austria) was used to analyze the differences in gene expression between the two groups. The filter conditions for significant differences in gene expression were |log2FC| >1.5 and P < 0.05. Then, the STRING protein–protein interaction database was used to screen the key differential genes [9], and the filtering conditions in the R package STRINGdb (version 2.30) were set to score >500 and connectivity >2.

Finally, the key screening method, Lasso regression and univariate Cox prognostic model were used. Due to the collinear relationship between the expression of certain genes, the constructed prognostic model exhibited overfitting. Lasso’s penalty coefficient (λ) was used to filter out the collinearity factors, and only the most representative factors were retained for the establishment of the prognosis model. The degree of Lasso regression complexity adjustment was controlled by the parameter λ. A greater value of λ corresponds to a greater penalty to obtain a model with fewer variables from a complex model. In this study, the parameter λ value was set to 1, and the cutoff value was set to 0.6 in the R package glmnet (version 4.02).

Analysis of the significance of signature genes

This study focuses on the biological significance of the final screening of key genes and the difference in responsiveness to chemotherapy and immunotherapy. In addition, the mutation status and methylation changes of the signature genes were also influencing factors of chemotherapy tolerance. According to the clinical information, the chemotherapy and immunotherapy patients were matched based on the barcodes of patient samples in TCGA database. Therefore, this study also analyzed data on the mutation status, methylation changes and site information of signature genes, and all the data were downloaded from the TCGA database. R packages maftools [10] (version 2.0) and Gviz [11] (version 1.28) were used to compare gene mutation differences and map methylation site information, respectively.

Verification of signature genes

According to the standards of the Ethics Review Committee of Qingdao Municipal Hospital and the conditions of informed consent provided by the patients, skin tissue samples from twenty patients with metastatic MM who underwent chemotherapy and immunotherapy and ten normal patients were collected. According instructions of the western blotting techniques [12], relative protein concentration of signature genes in tissues was detected. All antibody reagents were purchased from Abcam (Abcam.com). The SYBR Green labeling method was used to detect the gene expression level of signature genes by real-time qPCR [13], and the gene amplification primer sequence was obtained from PrimerBank [14]. ABI StepOne Plus was used to calculate the relative expression differences of signature genes. The primers used in this study are listed in Supplementary Table 1.

Quantification and statistical analysis

The western blot results of the signature genes were quantitatively analyzed by Quantity One (version 4.6.2, Bio-Rad Laboratories, Inc.,). The function of this tool can convert the gray value of the western blot image into data for comparative analysis. The real-time qPCR results were converted using the 2^−ΔΔct method to calculate the amplification fold changes, and Student’s t-test was used to test whether the differences between the two groups were statistically significant (P < 0.05).

Chemotherapy achieved a better overall survival expectancy

According to the screening rules, the gene expression data of 46 patients with metastatic melanoma treated with chemotherapy and 47 patients with immunotherapy were ultimately used in this study. Table 1 shows the differences in the treatment period and response to chemotherapy or immunotherapy. The results showed that there was a significant difference in the response rate between the chemotherapy and immunotherapy patients. In addition, taking the final survival status (death or survival) as the outcome indicator, the cumulative overall survival and disease-specific survival events of patients in the chemotherapy group were significantly higher than those in the immunotherapy group (Figure 1A and 1B), whereas significant differences were not observed in progression-free survival (Figure 1C). The results also showed that there was no difference in overall survival between male and female patients (Figure 1D).

Figure 1. Cumulative number of survived events in patients with metastatic melanoma undergoing chemotherapy and immunotherapy. (A–C) The Log-rank test for overall survival, disease specific survival, and progression-free interval survival curves of the chemotherapy and immunotherapy groups. (D) Survival curves of overall survival of male and female patients.

Anti-tumor properties of signature genes

Figure 2A shows the screening flowchart of signature genes, and based on |log2FC| > 1.5 and P < 0.05, 309 DEGs were used to construct the PPI protein interaction network. Subsequently, the genes in the marginal region were filtered out according to score >500 and degree >2 of the protein molecule. In this result, 50 genes (Supplementary Table 2) were used to construct a generalized linear Lasso-Cox regression model to further remove genes with collinear characteristics. Finally, seven signature genes were obtained: CCKBR, KCNJ11, NMU, MMP13, ITGA10, IGFBP1 and CEACAM5.

Figure 2. Flowchart of screening signature genes and its Kaplan-Meier survival curve and receiver operating characteristic curve (ROC). (A) Flowchart of screening signature genes. Including volcano plot of differentially expressed genes (DEGs), DEGs’ PPI network and Lasso-COX fitting model. (B) DEGs immune cell infiltration analysis, Pearson correlation test of Anti-tumor and Pro-tumor suppression based on immune cell type enrichment scores. (C) Log-rank test for survival curves of signature genes. (D) ROC curve, area under curve (AUC), sensitivity, specificity and accuracy values of signature genes.

According to The Cancer Immunome Atlas (TCIA) [15], a tumor immune cell infiltration analysis was performed, and the enrichment scores of these DEGs showed significant tumor suppressive effects in both the chemotherapy group and the immunotherapy group (Figure 2B). A comparison of the immune cell infiltration scores between the two groups showed that the immature dendritic cells (DCs) were significantly higher in the chemotherapy group than the immunotherapy group (P < 0.01), while the central memory CD8 T cells and type 2 T helper cells were significantly higher in the immunotherapy group than the chemotherapy group (P < 0.05) (Supplementary Figure 1A). The Pearson correlation test between immunoinfiltrating cells is shown in the attached Supplementary Figure 1B. In addition, these seven signature genes were divided into high and low groups according to their median expression levels. Significant anti-tumor effects were observed based on these groupings, and the high expression group of CCKBR and KCNJ11 showed better anti-tumor effects (Supplementary Figure 1C).

Functional enrichment in the PPI network

The functional enrichment analyses of the 50 genes screened by the PPI network included gene ontology terms (cell composition, molecular function and biological process), pathways (KEGG and Reactome) and protein domains (InterPro and Ffam database). (Supplementary Figure 2A). The biological functions and effects of these key genes are shown in Supplementary Figure 2. The cell composition category mainly includes cell membrane synapses, keratin filaments, collagen trimers, cell junctions and ion channel complexes (Supplementary Figure 2B). The molecular function category mainly includes ion channel activity and receptor binding, such as cytokine activity and receptor binding function, chemokine activity and receptor binding function, collagen and hormone binding function, etc., (Supplementary Figure 2C). The enrichment results of the integrated protein domains were mainly conserved sites of CXC chemokines, chemokine IL-8-like superfamily, hemopexin superfamily and intermediate filament (Supplementary Figure 2E). The signaling pathways mainly included IL-17, TNF and estrogen signaling, ECM, cytokine-cytokine and neuroactive-ligand receptor interactions, and nicotine and morphine addiction (Supplementary Figure 2F).

Differences in the survival of signature genes

As shown in Figure 2C, patients with high expression of the signature genes CCKBR, KCNJ11, ITGA10, and IGFBP1 had higher overall survival prognoses while patients with low expression of NMU, MMP13, and CEACAM5 had higher overall survival. This finding implies that patients with longer overall survival periods are positively correlated with anti-tumor scores for certain signature genes with high expression (CCKBR, KCNJ11, ITGA10, and IGFBP1) and other signature genes with low expression (NMU, MMP13 and CEACAM5). However, based on the anti-tumor properties described above, this property was significant only in the signature genes CCKBR and KCNJ11 (Supplementary Figure 1).

Then, receiver operating characteristic curves (ROCs) of these signature genes were drawn, and the area under the curve (AUC) was used to indicate the predictive value of the model. As shown in Figure 2D, the AUC values of the CCKBR, KCNJ11, NMU and MMP13 genes were greater than 0.7 and the AUC values of the other signature genes were between 0.7 and 0.6. In addition, the sensitivity, specificity and accuracy of the seven characteristic genes were also calculated to evaluate the predictive value of the signature genes in the ROC model.

Differences in mutations, methylation and treatment response

The signature genes CEACAM5, ITGA10 and MMP13 had the highest proportion of mutations among all 467 patients with cutaneous melanoma of the skin (Figure 3A and 3B). Although the frequency of mutations in the immunotherapy group was higher than that in the chemotherapy group, no significant difference was found between the two groups because the total sample of the study was only 93 cases (Figure 3C). A comparison of the protein domains of the missense mutations showed that the protein domains of the immunotherapy group mainly occurred in the box-labeled region while those of the chemotherapy group mainly occurred outside the box-labeled region.

Figure 3. Genetic variation and relative expression differences of signature genes between chemotherapy and immunotherapy. (A and B) Mutation frequency and type of signature genes. (C) Mutations and protein domain differences of signature genes between chemotherapy group and immunotherapy group. (D) Differences in relative expression levels (log2(CPM)) of signature genes between the two groups. (E) Comparison of complete response (CR), partial response (PR) and progressive disease (PD) differences in signature genes between the two groups. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001, NS: no significance.

Treatment responses include complete response, partial response, and disease progression, which are important indicators for judging the effectiveness of tumor treatment. Based on significant differences in signature genes between the two treatment groups (Figure 3D), they were grouped according to the treatment response, and significant differences in disease progression were observed for all genes except for the MMP13 gene and differences in complete response were mainly observed for the CCKBR, KCNJ11, ITGA10, IGFBP1 and CEACAM5 genes (Figure 3E). No significant differences were found for partial responses.

In addition, a group comparison was performed for methylation site changes of these signature genes. Significant differences were not found between the chemotherapy group and the immunotherapy group or between the treatment response groups. The comparative results are shown in Supplementary Figure 3.

Expression changes of signature genes in metastatic melanoma tissues

Using basic public data, this study verified the expression changes of signature genes in patients with metastatic melanoma after chemotherapy and immunotherapy. The gene expression level and protein concentration were expressed using real-time qPCR and western blotting, respectively, with reference to the relative expression level of the internal control. As shown in Figure 4A–4B, the comparison of the relative protein concentration of characteristic gene expression showed that the protein concentrations of CCKBR, ITGA10 and IGFBP1 in melanoma tumor tissue in the postoperative chemotherapy group were significantly higher than those in the immunotherapy group while the protein concentrations of KCNJ11, NMU, MMP13 and CEACAM5 in the chemotherapy group were significantly lower than those in the immunotherapy group. In addition, the relative differences in gene expression levels (Figure 4C) were consistent with the results of previous data analysis (Figure 3D). It is worth noting that the difference in protein concentration expressed by the KCNJ11 gene is opposite to the difference in gene expression level. Additionally, the hematoxylin-eosin (HE) staining results showed that the skin tissue structure of tumor patients was disordered and deeply stained compared with that of normal individuals (Figure 4D).

Figure 4. Validation of signature genes in patients with metastatic cutaneous melanoma undergoing chemotherapy and immunotherapy. (A and B) Western blot results of signature genes and t-test comparison results of gray values between the two groups. (C) T-test comparison results of the relative expression levels of signature genes in real-time qPCR between the two groups. (D) Hematoxylin-eosin (HE) staining: normal skin tissue, tumor tissues of metastatic melanoma patients receiving chemotherapy and immunotherapy. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001.