A lipid metabolism-related gene model reveals the prognosis and immune microenvironment of cutaneous melanoma

Major strategies of this study

A workflow demonstrating the strategies employed in this study is shown in Figure 1.

Figure 1:

The workflow diagram utilized in the present research.

Data gathering and preprocessing

Human SKCM tumors gene expression data (n=471) were extracted from The Cancer Genome Atlas (TCGA) database. Normal tissue sample data (n=812) were obtained from the GTEx repository (http://commonfund.nih.gov/GTEx/). The above data were merged with the shared genes expression levels. The combined expression spectrum was normalized by the normalize between arrays algorithm function in the ‘limma’ package of R3.6.1 Version 3.34.7 (https://bioconductor.org/packages/release/bioc/html/limma.html) [13]. The above data were used as the training set in the present research.

Additionally, we downloaded the GSE65904 dataset [14] from the NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/) database [15], [16], [17]. The above information was incorporated into the GPL10558 platform (Illumina HumanHT-12 V4.0 expression beadchip). This dataset contained information on 214 SKCM samples and the survival data of 210 patients. This dataset served as the validation dataset in this analysis.

Determination of differentially expressed genes related to lipid metabolism

LMGs were extracted from MSigDB (http://www.gsea-msigdb.org/gsea/msigdb/cards/KEGG_FATTY_ACID_METABOLISM.html) in the Gene Set Enrichment Analysis (GSEA) database [18] (http://software.broadinstitute.org/gsea/downloads.jsp). Standardized TCGA-GTEX expression data were used to determine the relevant genes expression levels in the samples.

Differentially expressed genes (DEGs) connected with lipid metabolism between SKCM and normal samples were determined using the ‘limma’ package of R3.6.1 Version 3.34.7 [14], ensuring the false discovery rate (FDR) lower than 0.05. A heatmap was used to display significant DEGs.

Generation of genes related to overall survival and lipid metabolism

Version 2.41-1 of the ‘survival’ package in R3.6.1 language [19] (http://bioconductor.org/packages/survivalr/) was adopted to filtrate the DEGs related to lipid metabolism in TCGA SKCM samples with clinical prognostic source. Then, survival-related LMGs with p <0.05 were generated by univariate Cox hazards regression analysis.

Establishment and validation of risk model

The LASSO algorithm in the R3.6.1 language ‘lars’ package version 1.2 [20] (https://cran.r-project.org/web/packages/lars/index.html) was applied to screen the optimal combinations of LMGs among those significantly associated prognosis and perform survival regression analysis. Then, the risk score (RS) model was established using the LASSO prognostic coefficient of each component in the optimized LMG composition and target gene expression level, and the RS calculation equation was presented below:

(Coef_genes represents the prognostic coefficient of the target gene under LASSO analysis, and Exp_genes represents the corresponding gene expression level.)

RS values were calculated for the TCGA training set and the GSE65904 validation dataset, and the samples were divided into high (RS score≥RS median value) and low (RS score<RS median value) risk sample groups. The Kaplan‒Meier curve method in the ‘survival’ survival package in R 3.6.1 version 2.41-1 [19] was used to assess the interrelation between the actual survival prognosis and the grouping conditions.

The clinical factors of samples in the TCGA training dataset in different risk groups were analyzed according to the high- and low-risk sample groups. Furthermore, the clinical information was compared utilizing the accurate intergroup Fisher test in R3.6.1.

Determination of differentially expressed genes related to risk groups and functional analyses

The TCGA training set samples were grouped into high- and low-risk groups in accordance with the previously determined RSs. Significant DEGs between the two groups were filtered using the ‘limma’ package in R3.6.1 Version 3.34.7, and the standards were FDR <0.05 and |log₂FC| >0.5. The significant DEGs were then analyzed by Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses to evaluate the genes with FDR values <0.05 using the Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) [20, 21].

Immune analysis

CIBERSORT (https://cibersort.stanford.edu/index.php) [22] was employed to evaluate diverse immune cell ratios in TCGA SKCM samples, and the comparison of immune cell distributions between the two risk groups was conducted. In addition, the correlation between the RS characteristic values according to the RS model and immune cells was carried out using the COR function in R3.6.1. The TIDE (Tumor Immune Dysfunction and Exclusion) database (http://tide.dfci.harvard.edu/) [23] was utilized to predict the response of each sample to immunotherapy, denoted by the TIDE score, and the TIDE score of each patient was determined. Then, a statistical analysis was conducted using a t-test to assess the variance in TIDE scores between the RS groups. Furthermore, the ‘estimate’ package in R3.6.1 (http://127.0.0.1:29606/library/estimate/html/estimateScore.html) [24] facilitated the calculation the stromal, immune, and ESTIMATE scores and tumor purity. Finally, the Kruskal–Wallis test in R3.6.1 was utilized to contrast the risk group estimated scores.

Analysis of drug sensitivity

Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) [25] data were combined with patient data, and ‘pRRophetic’ in R3.6.1 [26] (http://127.0.0.1:22402/doc/html/Search?objects=1&port=22402) was utilized to assess each patient’s susceptibility to diverse chemical medicines. IC₅₀ values, representing drug susceptibility, were compared between the risk groups for diverse drugs using the Kruskal–Wallis test.

HPA immunohistochemical analysis

The Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) [27, 28] was utilized to obtain immunohistochemical information on the proteins encoded by the screened genes in melanoma and normal tissue.

Cell lines and culture

Primary normal human epidermal melanocytes were separated from human prepuce tissues (human tissues were obtained from the patients of Hospital for Skin Diseases, Institute of Dermatology, Chinese Academy of Medical Sciences and Peking Union Medical College) refer to our laboratory protocols. This study was approved by the ethics committee of Hospital for Skin Diseases, Institute of Dermatology, Chinese Academy of Medical Sciences and Peking Union Medical College (2023017). Specifically, fresh foreskin tissue was soaked in iodophor for 5 min and washed with phosphorate buffered solution (PBS, 0.01 mol/L) 3 times. The subcutaneous tissue was separated, and the remaining tissue was cut into strips and added to a centrifuge tube containing dermal–epidermal separation solution (catalog number: 4942078001; Roche, Switzerland) at 4 °C overnight. The epidermis was isolated from the dermis with sterile ophthalmic tweezers, and the epidermis tissue was collected and cut into pieces. 0.25 % Trypsin-EDTA (catalog number: 25200056; Gibco, Canada) was applied to digest the tissues at 37 °C for 5 min. Complete medium (catalog number: 2201; Sciencell, United States) added with 10 % fetal bovine serum was added to neutralize the trypsin. Then, 100 µm and 70 µm filters were used to remove the excess ingredients, and the cell suspension was centrifuged (Centrifuge model: KDC-40; Zhongke Zhongjia company, Anhui, China) at 250×g at room temperature for 4 min. Finally, the cell precipitate was collected and resuspended in a melanocyte medium. The M14, A375, A2058, and SK-MEL-28 melanoma cell lines were restored at Hospital for skin Diseases, Institute of Dermatology, Chinese Academy of Medical Sciences and Peking Union Medical College (the M14, A375 and SK-MEL-28 cell lines were obtained from human skin tissues, and A2058 was obtained from human lymph node, all the four cell lines have correct STR identification and have no mycoplasma contamination). All melanoma cells were cultivated in Dulbecco’s Modified Eagle Medium (DMEM) (catalog number: C11995500BT; Gibco, China) added with 110 mg/L L-glutamine and 4,500 mg/L D-glucose combined with 1 % penicillin and streptomycin and 10 % fetal bovine serum. The above cells were incubated in cell incubators.

Real-time PCR analysis

RNAiso Plus (catalog number: 9108; Takara, Dalian, China) was applied to isolate total RNA according to its specification from primary human normal melanocytes (MCs) and melanoma cell lines (M14, SK-MEL-28, A2058, and A375) at passages 5–10. HiScript II Q RT SuperMix (catalog number: R22301; Vazyme Biotech, Nanjing, China) was utilized to synthesize cDNA from these RNAs. The obtained cDNA was subsequently served as a template for PCR magnification on a LightCycler^® 480 Instrument II (Roche, Basel, Switzerland) using ChamQ SYBR qPCR Master Mix (catalog number: Q33102; Vazyme Biotech, Nanjing, China). The 2^−ΔΔCt method was utilized for the calculation of gene expression levels which were then uniformization against GAPDH. The primers used are enumerated in Supplementary Table 1.

Western blotting

Total protein from MC, M14, A375, A2058, and SK-MEL-28 cells was obtained using RIPA lysis buffer. After centrifugation, the supernatant was measured with a BCA kit (catalog number: KGB2101; Keygen Biotech, Nanjing, China) for soluble protein content. In accordance with the protein concentration, the final protein solution was calculated by the principle of equal mass and equal volume, and RIPA lysate was used to supplement the lack of volume. Then 1/4 volume of 5× protein loading buffer (catalog number: P0285; Beyotime, Shanghai, China) was added to each specimen. All the final mixtures were bathed in boiling water for 6 min and centrifuged at low speed for further experiments. Protein mixtures were separated by SDS‒PAGE (catalog number: SLE015; Smart-lifesciences, Changzhou, China) and transferred to a PVDF membrane. Then, the membranes were blocked in EveryBlot Blocking Buffer (catalog number: 12010020; Bio-Rad, USA) for 5 min at room temperature and then incubated with primary antibodies (dilution ratio, 1:1,000; Proteintech, Wuhan, China; antibody codes: ADH4 16474-1-AP, ALDH7A1 10368-1-AP, HADH 19828-1-AP, and HADHA 10758-1-AP.) for 15 h at 4 °C. After being washed in TBST for 30 min, these membranes were immersed in the corresponding secondary diluting antibodies (HRP Goat Anti-rabbit IgG, dilution ratio, 1:2,000; MedChemExpress, Shanghai, China) for an appropriate time at room temperature. A chemiluminescence kit (catalog number: P0018; Beyotime, Shanghai, China) was applied for protein detection.

Statistical analysis

Each set of experiments was conducted in triplicate at least. GraphPad Prism9 and Image J version 1.8.0 software were used to get the result of the experiment analysis and visualization. Student’s t-test was used to analyze the differences between the two groups. “*”, “**” and “***” respectively showed p <0.05, p <0.01, and p <0.001. “ns” showed no obvious difference (p>0.05).

Data preprocessing and identification of differentially expressed genes

A total of 1,283 samples (812 normal and 471 cutaneous melanoma samples) and 60,370 genes were obtained. A total of 42 LMGs were obtained, all of which met the criteria of significant differential expression in the above two groups. The DEGs involved in lipid metabolism are shown in a heatmap (Figure 2A) and volcano map (Figure 2B). Univariate Cox regression analysis of the 42 genes is shown in a forest plot (Figure 2C).

Figure 2:

Identification of differentially expressed genes and these associated with prognosis. (A) Heatmap showed the expression of LMGs in TCGA SKCM tumor and GTEX normal group. (B) DEG volcano plot; red dots represent upregulated genes; blue dots represent downregulated genes. (C) Forest plot showed the univariate Cox regression of the LMGs. The red, blue and black squares respectively represented genes associated with risk, protection and non-significant prognosis. (D) KM curve with significant prognostic significance. Blue and red curves respectively represented low and high expression levels of genes. (E) Expression levels of prognostic significantly related LMGs in TCGA SKCM tumors and GTEX normal control samples, *** indicates p<0.001.

Identification of differentially expressed genes related to prognosis

A total of 10 DEGs with significant prognostic correlations were obtained. Based on the median gene expression level, the samples were classified into high- and low-expression groups, and the KM curve of 10 DEGs was constructed (Figure 2D). The expression levels in tumor and normal samples are displayed (Figure 2E).

The downregulation of ACAA2, ALDH7A1, CPT2, HADHA, HADHB, ADH7, and HADH expression and the upregulation of ACSL5, ACSL4, and ADH4 expression were speculated to be related to improved overall survival. In addition, the expression of genes such as ACAA2, ADH4, ADH7, and ALDH7A1 was predicted to be greater in normal tissue than in tumor tissue.

Establishment of lipid metabolism – related Gene Risk Model

The LASSO algorithm was applied to obtain the optimal gene combination for the 10 genes significantly correlated with prognosis (Figure 3). The following four optimized genes were obtained: ADH4, ALDH7A1, HADH, and HADHA. The RS computational equation was calculated according to the LASSO regression coefficients and four genes as follows:

Figure 3:

Establishment of the lipid metabolism-related Gene Risk Model. (A) LASSO parameter diagram. (B) The location where the horizontal and vertical black lines cross was the location where the optimal parameters (−4.113, 0.247), and the corresponding gene was the optimal gene combination.

Validation of the Gene Risk Model

The RS and survival time distributions of patients in the TCGA and GSE65904 datasets are shown in Figure 4A and B (upper and middle figures). The ROC curve results are shown in Figure 4A and B (bottom figure). The KM curve of each dataset is shown in Figure 4C and D, and. These results revealed a link between the actual prognosis and the two risk groups in both datasets.

Figure 4:

Validation of the Gene Risk Model. RS distribution (top), survival time status (middle) and ROC curves (bottom) of TCGA training (A) set and GSE65904 validation set (B). The figures in brackets showed the corresponding ROC curves and sensitivity. TCGA training set (C) and GSE65904 validation set (D) samples on the basis of the RS prediction model and prognostic KM curve. Blue and red curves respectively represented low- and high-risk samples.

Functional assessment

Samples in TCGA were split into high- and low-risk groups in accordance with the RS of the risk model. According to the ‘limma’ package, a total of 2,612 significant DEGs were identified between these two groups, followed by GO and KEGG analyses. A collection of functional categories relevant to immunity was highlighted in the GO terms MF (molecular function), CC (cellular component), and BP (biological process) (Figure 5). In the BP category, “immune response,” “chemokine-mediated signaling pathways,” and “T-cell activation” were the enriched terms. In the CC category, “external side of plasma membrane,” “extracellular region,” “plasma membrane,” and “receptor complex” were also enriched. In the MF category, “receptor binding,” “transmembrane signaling receptor activity,” and “cytokine activity” were the most enriched terms.

Figure 5:

Histogram of GO and KEGG analysis. BP (A), CC (B), MF (C), and KEGG (D) signaling pathways with significant DEGs among risk groups. The vertical axis showed the name of the item, the horizontal axis showed the number of genes, and the color of the column represented significance, the red color represented the most significance and the green color represented lesser significance.

Pathway maps representing reaction networks and molecular interactions were included in the KEGG analysis. The enriched immune pathways included “cytokine–cytokine receptor interaction,” “JAK-STAT signaling pathways,” “natural killer-mediated cytotoxicity,” “T-cell receptor signaling pathway,” and “NF-κB signaling pathways” (Figure 5). Hence, enrichment of immune-related cascades was acquired between the two groups (the top 10 were chosen according to the order of the FDR).

Immune infiltration analysis

CIBERSORT was used to extract the following six forms of immune cells with obvious variances in distribution between the two groups: CD4⁺ activated memory T cells, CD8⁺ T cells, M2 macrophages, activated mast cells, resting NK cells, and eosinophils (Figure 6A).

Figure 6:

Immune infiltration analysis. (A) Immune cells distribution with apparent discrepancy in distribution among the two risk groups. (B) The relationship between the immune cells distribution and the Risk Score characteristic values based on RS model (left: horizontal distribution in different risk groups and right: scatter plot of correlation with Risk Score characteristic values.). (C) The estimated scores of different risk groups. * represents p<0.05, ** represents p<0.01, and *** represents p<0.001.

In comparison with those in the high-risk group, the low-risk group exhibited a higher presence of CD8⁺ T cells and CD4⁺ activated memory T cells, but a lower abundance of resting NK cells, M2 macrophages, eosinophils, and activated mast cells. In R3.6.1, the RS characteristic values according to the RS model were negatively related to the distribution of CD4⁺ activated memory T cells (p=0.0001499) and CD8⁺ T cells (p=0.01181) but positively related to the distributions of resting NK cells (p=0.02363), M2 macrophages (p=0.007852), and activated mast cells (p=0.002235). Furthermore, the group at lower risk exhibited decreased tumor purity but elevated levels of immune, stromal, and ESTIMATE scores (Figure 6B).

Comparison of immune response and drug sensitivity

A total of 148 non-responder and 311 responder samples were predicted, revealing that the proportion of responder samples in the group at lower risk was greater than that in the group at higher risk (Figure 7A). Furthermore, patients in the group at lower risk exhibited increased sensitivity to AZD7762, camptothecin, gefitinib, gemcitabine, nilotinib, rapamycin, and vinblastine (Figure 7B).

Figure 7:

Comparison of immune response and drug sensitivity analysis. (A) Distribution of TIDE prediction of immune response in the two risk groups. (B) Contrast of the IC₅₀ values of various drugs between different risk groups by Kruskal–Wallis test.

Validation of the genes by RT‒qPCR, western blotting, and HPA immunohistochemical analysis

The RT-qPCR and Western blot confirmed that ADH4 and ALDH7A1 expression were downregulated in the M14, A375, A2058, and SK-MEL-28 cell lines, but upregulated in the MC cell line, whereas HADH was increased in the M14, A2058, and SK-MEL-28 cell lines, and HADHA was increased in the SK-MEL-28 and A375 cell lines (Figure 8). These differences were significant (p<0.05). The HPA immunohistochemical results supported the downregulation of ADH4 and ALDH7A1 expression and the upregulation of HADH and HADHA expression in melanoma.

Figure 8:

Gene expression validated by RT-qPCR, Western blot and HPA immunohistochemical analysis. (A) RT-qPCR revealed the relative mRNAs expression levels of ADH4, ALDH7A1, HADH and HADHA. Data are shown as mean±standard deviations, n=3; *p<0.05; **p<0.01; ***p<0.001; ns means p≧0.05. (B) Western blot and their quantification revealed the relative protein levels of ADH4, ALDH7A1, HADH and HADHA. (C) The protein products expression of ADH4, ALDH7A1, HADH and HADHA were obtained by HPA immunohistochemical analysis.