Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblasts Signature for Predicting Prognosis and Immunotherapeutic Responses in Skin Cutaneous Melanoma

doi:10.21203/rs.3.rs-3330660/v1

Download PDF

Research Article

Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblasts Signature for Predicting Prognosis and Immunotherapeutic Responses in Skin Cutaneous Melanoma

https://doi.org/10.21203/rs.3.rs-3330660/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Objective: Cancer-associated fibroblasts (CAFs) are the most prominent cellular component in the stroma of skin cutaneous melanoma (SKCM). The aim of this study was to explore the CAFs-related genes and to develop a CAFs-related model to predict the prognosis and immunotherapeutic efficacy of SKCM.

Method: We collected mRNA expression and clinical information for 449 SKCM patients from the TCGA (The Cancer Genome Atlas) database and 210 patients in the GSE65904 dataset of the GEO (Gene expression Omnibus) database. CAFs infiltrations were quantified by the estimate the proportion of immune and cancer cells (EPIC) method. Weighted gene coexpression network analysis (WGCNA) was used to identify genes CAFs-related genes. A CAFs risk signature was then developed using the univariate and least absolute shrinkage and selection operator method (LASSO) Cox regression model. Gene set enrichment analysis (GSEA) and single sample gene set enrichment analysis (ssGSEA) were applied to elucidate the molecular mechanisms. The tumor immune dysfunction and exclusion (TIDE) algorithm is further used to assess immunotherapeutic responses.

Results: A six-gene (NOTCH3, HEYL, BGN, COL5A1, SULF1, COL1A1) CAFs prediction model was constructed. SKCM patients were divided into high– and low–CAFs-risk groups according to the median CAFs risk score. patients in the CAFs high-risk group had a significantly worse prognosis. GSEA showed that the ECM receptor-interacting and basal cell oncogene set was highly and significantly enriched in the high–CAFs-risk group. The ssGSEA results showed that CAFs risk scores were significantly positively correlated with basal cell carcinoma, melanogenesis and Notch signaling pathways, and significantly negatively correlated with regulation of autophagy and Toll-like receptor signaling pathways. ssGSEA results also showed that patients in the high–CAFs-risk group had a poorer response to immunotherapy.

Conclusion: The six CAFs signature genes explored in this study not only reliably predicted prognosis, but also assessed clinical immunotherapy response in SKCM patients, which may have important clinical implications for tailoring anti-CAFs therapy and immunotherapy regimen for each SKCM patient.

skin cutaneous melanoma

cancer-associated fibroblasts

weighted gene co-expression network analysis

prognosis

immunotherapy

Skin cutaneous melanoma (SKCM) is a highly heterogeneous and aggressive malignant tumor(1). Although SKCM is not as common as other skin cancers (basal cell carcinoma, squamous cell carcinoma), accounting for less than 5% of all skin cancers diagnosed in the United States each year. However, its lethality rate is very high. SKCM is responsible for approximately 73% of all skin cancer-related deaths(2, 3). According to statistics from 2014 to 2018, although the incidence of SKCM decreased by about 1% per year, the number of new cases and the death cases per year remained large(4). The number of new cases of SKCM worldwide in 2018 was 288,000 with 61,000 death cases (5). In recent years markers such as BRAF and NRAS mutations and immune checkpoints such as PD-L1 have contributed significantly to the accurate diagnosis and prognosis of SKCM patients(6), and anti-CTLA4 combined with anti-PD-1 checkpoint inhibition has improved the 5-year relative survival rate of advanced SKCM from 15–30% in 2004(4), but we still need to actively develop novel SKCM biomarkers that can predict SKCM prognosis and immune therapy response with good predictive performance to improve the prognosis and treatment of advanced SKCM.

The tumor immune microenvironment (TME) is the cellular environment in which tumor cells exist. TME includes malignant tumor cells, immune cells, blood vessels, extracellular matrix, fibroblasts, lymphocytes, myeloid-derived inflammatory cells, and signaling molecules that are involved in the tumorigenesis and cancer progression(7). Among the cells in TME, cancer-associated fibroblasts (CAFs) are a mesenchymal-derived cell type that maintain the structural framework of tissues(8). Studies have shown that most CAFs are produced by the activation of normal fibroblasts (NFs). In contrast to NFs, CAFs have increased proliferation, enhanced extracellular matrix (ECM) production, and unique cytokine secretion (e.g. stromal cell-derived factor 1, vascular endothelial growth factor, platelet-derived growth factor, hepatocyte growth factor) and are the most prevalent and essential stromal cells associated with poor prognosis in many cancers(8, 9). During SKCM pathogenesis, CAFs interact with SKCM cells to promote tumor growth, invasion, dense structure formation, angiogenesis, metastasis and even drug resistance(7, 10).

Weighted gene co-expression network analysis (WGCNA) enables combining highly coordinately expressed genes into several gene modules and investigating the relationship between the modules and the phenotype of interest(11). WGCNA has been successfully applied for identifying CAFs markers(12–14). As of now, however, CAFs and stromal infiltrates have not been analyzed by WGCNA in SKCM patients. We used WGCNA on both transcriptomic datasets collected from the GEO and TCGA databases to detect the hub modules most associated with stroma and infiltration of CAFs. Then by univariate and LASSO Cox regression analysis, we identified NOTCH3, HEYL, BGN, COL5A1, SULF1, and COL1A1 as CAFs signature genes predicting prognosis in SKCM patients and constructed models capable of predicting prognosis and treatment response in SKCM. Our results suggest that CAFs signature genes can be used as targets for new anti-CAFs therapy in SKCM.

2.1 Data Acquisition and Methodology

472 SKCM cases were obtained in the TCGA database (https://portal.gdc.cancer.gov/), of which 449 had complete survival data. 241 SKCM cases were obtained in the GSE65904 dataset(15, 16) of the GEO database (https://www.ncbi.nlm.nih.gov/geo/), of which 210 patients with complete survival information were used as an independently validated external test set.

2.2 Infiltration estimation of CAFs and analysis of tumor immune microenvironment

CAFs infiltration score were separately estimated via four methods: cell-type deconvolution (constrained least square optimization)–based Estimate the Proportion of Immune and Cancer cells (EPIC) algorithm(17), gene signature enrichment–based xCell algorithm(18), marker genes expressions–based microenvironment cell populations counter (MCP-counter)(19), and Tumor Immune Dysfunction and Exclusion (TIDE) algorithms(20).The first three methods are respectively achieved by the R packages "EPCI(1.1.5)", "xCell(1.1.0༉", "MCPcounter(1.2.0༉". The TIDE method was implemented through the website http://tide.dfci.harvard.edu/.

the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was applied to calculate the stromal score via R package "ESTIMATE (1.0.13)", which indicates the stromal infiltrating levels of each tumor sample(21).

2.3 Weighted gene co-expression network analysis

WGCNA is a systematic bioinformatics algorithm capable of integrating highly coordinated expressed genes into several gene modules and investigating the relationship of modules to phenotypes of interest. An appropriate soft power threshold (β) was chosen to find the best balance to generate the largest number of modules without loss of gene module membership (MM).WGCNA was performed using the R package "WGCNA(1.72-1)" (11)to screen out the CAFs-related hub genes.

2.4 Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Analyses

To explore the biological functions and pathways of the CAFs-related hub genes, we performed gene ontology (GO) analysis and Kyoto Gene and Genome Encyclopedia (KEGG) pathway enrichment analysis using the R package "clusterProfiler (4.6.0)"(22, 23). p < 0.05 was considered statistically enriched.

2.5 Prediction model construction and validation based on CAFs-related hub genes

The TCGA-SKCM cohort was selected for CAFs risk model construction, while GEO cohort were assigned to the validation cohort. The univariate Cox regression model was performed to identify prognostic CAFs hub genes on overall survival (OS); genes with p < 0.05 were subsequently put into LASSO Cox regression analysis with 1,000 iterations for gene reduction via the R package "glmnet(4.1-6)"(24). Then the risk model of CAFs was constructed and the CAFs-related risk score (RS) was calculated, and RS was used as a new prognostic feature. RS calculation formula: \(\text{R}\text{S}={\sum }_{i}Cofficient（{miRNA}_{i}）\times Expression\left({miRNA}_{i}\right)\). Patients were divided into low- and high-CAFs-risk groups according to RS, and survival analysis was performed by the R package "survival (3.5-0)", and multiple Cox regression analysis was performed to verify that RS is an independent risk factor for the prognosis of SKCM. To verify the prognostic value and accuracy of the model, we used the R packages "survminer (0.4.9)", "timeROC (0.4)", and "pheatmap (1.0. 12)" package to plot K-M curves, risk curves and receiver operating characteristic (ROC) curves. Similarly, the model was validated in the GSE65904 validation cohort.

2.6 Enrichment analysis of CAFs signature genes

To investigate the differences in biological functions between low- and high-CAFs-risk groups, we performed GSEA and ssGSEA analysis using the R packages "clusterProfiler (4.6.0)" and "GSVA (1.46.0)"(25). The internal reference gene set (c2.cp.kegg.v7.0.symbols.gmt) was downloaded from the Molecular Characterization Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb/)(26). FDR (false discovery rate) < 0.25, P < 0.05, considered pathways statistically enriched. Spearman correlation analysis was used to assess the correlation between CAFs risk scores and gene set enrichment scores.

2.7 Immune infiltration and immune checkpoint analysis

Immune cell infiltration was assessed by ssGSEA and CIBERSORT algorithms(27) based on the R packages "estimate (1.013)", "limma (3.54.1)", "ggpubr (0.50)", "GSVA (1.46.0)". "GSVA (1.46.0)" Immune cell infiltration was assessed by ssGSEA and CIBERSORT algorithms. The immune infiltration status of the samples was determined by calculating the normalized enrichment score using single sample gene set enrichment analysis (ssGSEA). The relative proportions of the 22 immune cell types in the TME were determined by the CIBERSORT method. The genetic marker for these 22 immune cell types is called LM22, which contains 547 genes capable of distinguishing between different immune cell types. The correlation between RS and immune checkpoints was further analyzed based on the R package "ggExtra (0.10.0)" and "corrplot (0.92)". Immune cell genes were obtained from TISIDB, an integrated tumor and immune database.

2.8 Immunotherapy response prediction

To further validate the predictive value of CAFs signature genes for immunotherapy response, the TIDE algorithm was applied to assess the immunotherapy sensitivity of SKCM patients. The TIDE (http://tide.dfci.harvard.edu/) online algorithm was used for immune checkpoint blockade treatment response prediction. The difference in efficiency between the low- and high-CAFs-risk groups was then tested using a chi-square test, and the predicted efficacy of the risk profile of CAFs was assessed by ROC curve and area under the curve (AUC) values.

2.9 Validated by Cancer Cell Line Encyclopedia (CCLE) and Human Protein Atlas (HPA) databases

For cellular level validation, we downloaded mRNA expression of target genes in fibroblast and skin cell lines from the CCLE database (https://portals.broadinstitute.org/ccle) (28)and analyzed the differences between them by heat map and Wilcoxon assay. For protein level validation, we downloaded immunohistochemical images of SKCM tissues and normal skin tissues from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/)(29), which allowed direct observation of target protein localization and differences in target gene expression between them.

2.10 Statistical Analysis

All statistical analyses and visualizations were performed using R software v4.2.2 (https://www.r-project.org/).The median CAFs risk score was the cut-off value for classifying SKCM patients into low- and high-CAFs-risk groups, and the Wilcoxon test was applied for pairwise comparisons. The Kaplan–Meier curve with the log-rank test was adopted for overall survival comparisons via R package "survival(3.5-0)" and "survminer(0.4.9)". p < 0.05 was regarded as statistically significant.

3.1 Patients with SKCM with higher infiltration of CAFs and lower stromal scores had worse OS

CAFs infiltrations were multiply predicted by EPIC, xCell, MCP-counter, and TIDE methods, and the stromal score was calculated by the estimate algorithm. Their prognostic values on OS were evaluated via log-rank tests; Kaplan–Meier curves indicated that higher CAFs infiltrations and lower stromal scores were notably correlated with poorer OS of SKCM patients in both TCGA-SKCM (Fig. 1A) and GSE65904 (Fig. 1B) cohorts, which highlighted the importance of further exploration of CAFs- and stromal-related genes for SKCM. In this paper, the infiltration of CAFs estimated and stromal scores by the MCP-counter method were randomly used as phenotypic data for subsequent analysis. The infiltration data of CAFs evaluated by the other 3 methods were used for external validation.

3.2 Co-expression network of CAFs and stromal scores

To screen for modules closely related to CAFs, co-expression networks and central genes for infiltration and stroma scoring of CAFs were constructed and detected by the R package "WGCNA(1.72-1)". In this study, we used WGCNA to analyze RNA sequencing data, CAFs infiltration and stromal score of each SKCM patient in TCGA-SKCM cohort and GSE65904 cohort. Genes with the top 5,000 of median absolute deviation (MAD) were first chosen as the input genes for network constructions in both TCGA-STAD and GSE84437 cohorts. Then, the Pearson’s correlation similarity matrix between any gene pairs was calculated and increased to soft-thresholding power β based on the scale-free topology network criterion. Subsequently, the adjacency matrix was clustered using topological overlap measure (TOM) and dissimilarity (1-TOM) between genes, and the dynamic tree cut algorithm was applied to the dendrogram for module identification with the mini-size of module gene numbers set as 30 and a cut height of 0.9. The soft threshold༈β༉in the TCGA-SKCM cohort was determined to be 5 (Fig. 2A), and the hierarchical clustering tree showed that seven co-expression models were clustered (Fig. 2C), with Module brown being strongly associated with CAFs infiltration score (P = 6e-163, R = 0.89) and stroma score (P = 6E-69,R = 0.7) (Fig. 2E). The soft threshold༈β༉in the GSE65904 cohort was determined to be 4 (Fig. 2B), and the hierarchical clustering tree showed that seven co-expression models were clustered (Fig. 2D), where Module green was strongly associated with CAFs infiltration score (P = 3e-60, R = 0.85) and stroma score (P = 2E-26,R = 0.64) (Fig. 2F).There are 550 genes in Module brown and 418 genes in Module green, and we measured gene significance (GS) for traits and module membership (MM) (MM indicates the correlation between module and gene expression) of individual genes in these two hub modules. In Module brown, the scatter plot illustrates a strong correlation between MM and GS of CAFs genes (Cor = 0.92, P < 1e-200) (Fig. 2G). In Module green, there was also a strong correlation between MM and GS of CAFs genes (Cor = 0.9, P = 4.2e-152) (Fig. 2H). Then, 207 core genes were screened from the TCGA-SKCM cohort and 123 core genes from the GSE65904 cohort using MM > 0.6 and GS > 0.4 as threshold criteria, respectively.

3.3 Enrichment analysis of hub genes

As shown in Fig. 3, 60 genes were detected by the crossover of two hub gene sets and visualized by Venn diagrams. Subsequently, we performed GO and KEGG analysis on these 60 genes. Subsequently, we performed GO and KEGG analysis on these 60 genes. extracellular matrix organization, extracellular structural organization and external encapsulating structure organization were the major biological process enrichment terms. Collagen-containing extracellular matrix and extracellular matrix structural constituent were the major cellular component and molecular function enrichment terms, respectively. Protein digestion and absorption was the main pathway for KEGG enrichment.

3.4 6 CAFs signature genes used to construct predictive models

The TCGA-SKCM cohort was used as the test group and the GSE65904 cohort as the external validation group. Eight OS-related genes with p < 0.05 were screened by univariate Cox regression analysis of 60 core genes (Fig. 4A). LASSO Cox regression analysis was subsequently performed (Fig. 4B,4C). A prediction model consisting of six CAFs signature genes, namely NOTCH3, HEYL, BGN, COL5A1, SULF1, and COL1A1, was finally constructed (Table 1). Of the six genes, the SULF1 gene is a protective factor and the remaining five genes are risk factors. RS calculation formula: \(\text{R}\text{S}={\sum }_{i}Cofficient（{miRNA}_{i}）\times Expression\left({miRNA}_{i}\right)\). Multifactorial Cox regression analysis verified that RS could be used as an independent risk factor to predict the prognosis of SKCM patients (P < 0.01) (Fig. 4E). Risk curves showed that mortality in SKCM patients increased with increasing RS and the SULF1 gene was highly expressed in the low-CAFs-risk group and the other five target genes were highly expressed in the high-CAFs-risk group (Fig. 5A). Patients were divided into high- and low-CAFs-risk groups according to RS values, and Kaplan-Meier curves showed significantly shorter OS in the high-CAFs-risk group (HR = 2.162,95%CI:1.646–2.841, log-rank P < 0.01) (Fig. 5C). And the area under curve (AUC) values (0.601,0.627,0.659)of the ROC curve for 1-,2- and 3-year survival of this model indicated its high sensitivity and specificity for survival prediction (Fig. 5E). Also, the model was validated in the GSE65904 cohort (Fig. 5B,5D,5F). These results suggest that marker genes associated with CAFs are important prognostic markers for SKCM.

Table 1

CAFs signature genes that constitute the prediction model
Gene	Coef
NOTCH3	0.131826493
HEYL	0.076045351
BGN	0.016292909
COL5A1	0.087442816
SULF1	-0.256555994
COL1A1	0.051149458

3.5 CAFs signature genes are positively associated with CAFs infiltration and CAFs markers

To further validate the robustness of the CAFs model as a predictor of CAFs infiltration, we performed Spearman correlation analysis of the RS with the stromal score and the CAFs infiltration estimated by the four methods. We could observe that the RS in the TCGA-SKCM cohort (Fig. 6A) and the GSE65904 cohort (Fig. 6B) are positively correlated with the infiltration of CAFs estimated by multiple methods. RS is negatively correlated with stroma scores in the TCGA-SKCM cohort (R=-0.017) (Fig. 6A) and positively correlated with stroma scores in the GSE65904 cohort (R = 0.084) (Fig. 6B). With addition, we observed that RS and the expression levels of six CAFs signature genes in the TCGA-SKCM cohort (Fig. 6C,6E) and GSE65904 cohort (Fig. 6D,6F) were highly positively correlated with the large number of CAFs markers collected from published papers.

3.6 GSEA and ssGSEA analysis of 6 CAFs signature genes

To further elucidate the functional enrichment of CAFs signature genes, GSEA analysis and ssGSEA analysis were performed between high- and low-CAFs-risk groups in the TCGA-SKCM cohort. KEGG enrichment analysis showed that the high-CAFs-risk group was mainly enriched in extracellular matrix receptor interactions, basal cell carcinoma, and 3 cardiac diseases (Fig. 7A). The low-CAFs-risk group was mainly enriched in antigen processing and presentation, chemokine signaling pathways, cytokine cytokine receptor interactions, graft-versus-host disease, and allograft rejection (Fig. 7B). In addition, ssGSEA results showed that RS were significantly positively correlated with basal cell carcinoma (R = 0.46, p < 2.2e-16), melanogenesis (R = 0.43, p < 2. 2e-16), Notch signaling pathway (R = 0.49, p < 2.2e -16); significantly negatively correlated with regulation of autophagy (R = -0.46, p < 2.2e-16), Toll-like receptor signaling pathway (R = -0.43, p < 2.2e-16) (Fig. 7C).

3.7 Correlation analysis of RS with tumor immune microenvironment and immune checkpoints

The percentage of immune and stromal cells in the tumor tissue of SKCM patients was first calculated using the R package "estimate (1.0.13)" and analyzed for differences between high- and low-CAFs-risk groups. The results of the analysis showed that the low-CAFs-risk group was accompanied by higher immune score and composite score (P < 0.05), the high-CAFs-risk group was accompanied by higher tumor purity score (P < 0.05), and there was no statistical significance in stromal score between the high- and low-CAFs-risk groups (P = 0.7) (Fig. 8A). This suggests that the low-CAFs-risk group had more active immune function than the high-CAFs-risk group. For further immunological analysis, we investigated the infiltration of different immune cells in 22 using the CIBERSORT algorithm (Fig. 8B) and showed that RS was significantly associated with the following immune cells: Macrophages M0 (R = 0.34, p = 1.2e-08); NK cells resting (R = 0.2, p = 0.00085); Macrophages M1(R = -0.27, p = 4.4e-06); T cells CD8 (R = -0.34, p = 7.9e -09); T cells CD4 memory activated (R = -0.33, p = 2e-08) (Fig. 8C). In addition, ssGSEA analysis scored immune cells and immune-related functions and performed a difference analysis of the scores, which showed that the vast majority of immune cells and immune functions differed statistically significantly (p < 0.05) in the high- and low-CAFs-risk groups, and scored higher in the low-CAFs-risk group, while mast cells had higher scores in the high-CAFs-risk group (Fig. 9A), which is more visually illustrated in Fig. 9B. We analyzed the correlation between RS and immune checkpoints using the R package "ggExtra" and "corrplot". The results showed that all 34 common immune checkpoints differed statistically significantly (p < 0.05) in the high- and low-CAFs-risk groups, and the genes were enriched in the low-CAFs-risk group. Negative correlations were found between RS and various immune checkpoints such as CD274, IDO1, PDCD1LG2, CTLA4, HAVCR2, and PDCD1 (p < 0.05) (Fig. 9C,9D).

3.8 Analysis of the correlation between RS and immunotherapy

Immunotherapy with immune checkpoint inhibitors offers hope for SKCM patients. We applied the TIDE method to assess whether RS could be used as a predictor of immunotherapy in SKCM patients. The results showed that the difference in TIDE scores was statistically significant between the high- and low-CAFs-risk groups (p < 0.001), with higher TIDE scores in the high-CAFs-risk groups than in the low-CAFs-risk groups. Patients in the low-CAFs-risk groups were much more sensitive to immunotherapy than in the high-CAFs-risk groups (p = 0.005). The area under the ROC curve AUC value was 0.781(95% CI:0.729–0.826), indicating that our model of CAFs performed well in predicting response to immunotherapy.

3.9 Multidimensional validation of CAFs signature genes in CCLE and HPA databases

Based on the CCLE database, we verified that the mRNA expression of six CAFs signature genes (NOTCH3, HEYL, BGN, COL5A1, SULF1, COL1A1) was higher in fibroblast cell lines than in SKIN cell lines (Wilcoxon test, p < 0.001) (Fig. 10A,10B). In addition, to determine the protein expression profile of these CAFs signature genes, we analyzed immunohistochemical images from the HPA database. The images showed that the SULF1 protein was deeply stained in normal skin tissue and the other five proteins were deeply stained in the SKCM matrix (Fig. 10C). These validations suggest that these 6 genes may be markers specific for CAFs.

Studies have shown an interactive relationship between SKCM and CAFs. SKCM mainly contributes to the formation of CAFs in a paracrine manner(30). As SKCM develops, SKCM cells secrete a variety of chemokines, such as platelet-derived growth factor CC(31, 32), transforming growth factor-β(33), and Nodal protein, which induce NFs cells to migrate toward, surround, and infiltrate the growing SKCM tumor. These signals also induce epigenetic and biochemical changes of NFs in the stroma, resulting in the formation of CAFs. Extracellular vesicles secreted by SKCM cells can also induce reprogramming of NFs into CAFs. Extracellular vesicles are extracellular organelles, such as exosomes and cellular microvesicles, released or shed by eukaryotic cells into the tumor microenvironment. They contain various bioactive proteins, messenger RNAs and microRNAs (microRNA, miRNA/miR), which can be taken up by each other between cancerous cells, CAFs, mesenchymal stem cells and other cells, and become the key to intercellular communication(34). SKCM cells secrete extracellular vesicles that mediate the conversion of NFs to CAFs and regulate the expression of vascular cell adhesion molecules, thereby enhancing the adhesion of SKCM cells to CAFs. At the same time, CAFs in turn directly or indirectly affect SKCM progression. Various in vivo and ex vivo models have demonstrated the promotional effect of CAFs on SKCM growth. Studies have shown that interleukin-6, interleukin-8, CXC chemokine ligand 1, insulin-like growth factor 2, fibroblast growth factor 7 and transforming growth factor-β₃ are involved in the interaction between CAFs and SKCM cells, and that the above cytokines are involved in immune responses and induce inflammation and stimulate SKCM growth(10, 32, 35). CAFs have a critical role in tumorigenesis, growth, immunity, energy metabolism, angiogenesis, recurrence and metastasis(36).

In this study, a six-gene (NOTCH3, HEYL, BGN, COL5A1, SULF1, COL1A1) prediction model for CAFs was constructed, which had prognostic accuracy. Each of the six CAFs signature genes involved in the construction of the prediction model had a different role. NOTCH3 is one of the mammalian NOTCH proteins(37), and HEYL as a downstream target molecule of NOTCH3 (38). NOTCH3 dysregulation is associated with a variety of cancers, and it has been shown to affect tumor aggressiveness and resistance to chemotherapy(37). Studies have shown that squamous cell carcinoma has higher NOTCH3 expression compared to surrounding normal epithelial tissue, and higher NOTCH3 expression is associated with poorer prognosis in squamous cell carcinoma(39). COL1A1 and COL5A1 are members of the collagen fibronectin family, type I collagen α1 and type V collagen α1, respectively. as the most important components of the ECM, the collagen family is the basis for maintaining the structural integrity and tensile strength of human tissues and organs(40, 41). In malignant tumors, members of the collagen family are able to regulate the polarity, migration and signaling of tumor cells by modulating the biochemical and physical properties of TME(42). BGN(the proteoglycan-I gene), which encodes a protein modified to form the core of a glycoprotein, is a key component of the ECM, involved in scaffolding collagen protofibrils and mediating cellular signaling(43). Studies have already demonstrated the role of BGN in tumor proliferation, adhesion and invasion(44). Enforced expression of SULF1(human sulfate esterase 1), decreases cell proliferation, migration and invasion and is therefore downregulated in most cancer cells. SULF1 also promotes drug-induced apoptosis of cancer cells in vitro and inhibits tumorigenesis and angiogenesis in vivo.(45)

GSEA showed that ECM receptor interactions, basal cell carcinoma gene set were highly and significantly enriched in the high CAFs risk group. The ssGSEA results also showed that CAFs risk score was significantly positively correlated with basal cell carcinoma, melanogenesis, and Notch signaling pathway, and significantly negatively correlated with regulation of autophagy, Toll-like receptor signaling pathway. CAFs are the most efficient cells for deposition and remodeling of ECM in the tumor microenvironment. CAFs can release many growth factors and proteases (e.g. matrix metalloproteinases) and increase ECM protein expression, resulting in a strong ECM biosynthesis and deposition capacity. Thus, the dense, fibrotic nature of solid tumors is thought to be a direct result of fibroblast infiltration and ECM protein deposition. This result could prevent therapeutic agents from reaching cells at the tissue core of SKCM, thus contributing to tumor drug resistance(9). Studies have shown that NOTCH signaling is involved in the differentiation of CAFs. In keratin-forming cell tumors, NOTCH signaling deficiency promotes differentiation of CAFs and further tumorigenesis. However, in colon and prostate cancers, CAFs differentiation is caused by elevated NOTCH signaling. Furthermore, CAFs activate NOTCH signaling in cancer cells to promote various malignant behaviors, including cancer stem cell phenotype, chemoresistance, metastasis, and disease recurrence(46). Cellular autophagy is an important intracellular biometabolic pathway with a highly conserved evolutionary process, and it degrades its own cellular contents through lysosomes, which are part of the lysosomal degradation pathway. The organism removes, degrades and recycles waste or damaged biomolecules and organelles in the cell through cellular autophagy. The presence of autophagy facilitates the recycling of intracellular components and plays an important role in maintaining intracellular homeostasis(47, 48). In the early stages of tumorigenesis, autophagy may help control or kill cancer cells. In contrast, for advanced tumor cells, autophagy can promote tumor progression(49). The Toll-like receptor (TLR) family is a highly conserved and critically important component of the innate immune system. TLRs expressed extracellularly (e.g. TLR1/2/4/5/6) can be activated by recognition of extracellular irritations. Activation of TLRs not only promotes the expression and release of pro-inflammatory cytokines and mediators, but also directly affects tumor progression(50). Among them, TLR3 is the main molecule involved in activating Th1 immunity, secreting cytokines and chemokines, and inducing apoptosis, which can promote cell death in a variety of tumor cells, including SKCM, breast cancer, colon cancer, bladder cancer, head and neck cancer, pharyngeal cancer, kidney hepatocellular carcinoma, and lung cancer(51).

To reveal the relationship between RS and immune infiltration, we performed a series of immune correlation analyses. The results showed that immune cell infiltration was higher in the low-CAFs-risk group than in the high-CAFs-risk group, and RS was significantly positively correlated with macrophage M0 and negatively correlated with CD8 + T cells. It has been suggested that increased infiltration of macrophage M0 is associated with decreased OS and increased tumor stage(52). Macrophage M0 can express some cytokines and enzymes to inhibit the recruitment and activation of T cells, thus enhancing the resistance of tumor tissue to immunotherapy(53). Therefore, macrophage M0 may play an important role in the development of SKCM and may have therapeutic targeting. In addition, ssGSEA analysis showed that mast cells were enriched in the high-CAFs-risk group. Mast cells produce various mediators that stimulate angiogenesis, induce extracellular matrix breakdown, and stimulate tumor growth(54). This study also showed a negative correlation between RS and various immune checkpoints such as CD274, IDO1, PDCD1LG2, CTLA4, HAVCR2, and PDCD1 (p < 0.05), and the TIDE scores were higher in the high-CAFs-risk group than in the low-CAFs-risk group, implying that immune checkpoint inhibitors were more effective in treating the low-CAFs-risk group.

The present study has some limitations. This study is a retrospective bioinformatics analysis based on two publicly available gene expression data from TCGA and GEO, and the true prognostic and therapeutic value of the CAFs prediction model should be cross-validated in multicenter and perspective data. In addition, the specific biological roles of the six genes involved in the construction of the CAFs prediction model in SKCM need to be validated by molecular and animal experiments.

In summary, we confirmed higher infiltration of CAFs in TME was correlated with worse prognosis in SKCM and identified NOTCH3, HEYL, BGN, COL5A1, SULF1 and COL1A1 as novel prognostic CAFs biomarkers by constructing an integrated co-expression network of infiltrated CAFs and stromal scores. A CAFs model comprising of these six markers was constructed and validated, which could precisely predict prognosis, CAFs infiltrations and immunotherapy responses in SKCM patients, which might provide novel insights for targeting CAFs therapeutic strategies in SKCM.

Acknowledgments

The authors are grateful to all researchers who uploaded experimental data in the TCGA database and the GEO database, and to the staff who maintain the database operations.

Authors^，contributions

QW and JL designed the study. JL participated in data selection, assembly and analysis. JL, JH and YH were involved in drafting the manuscript. QW and ZY revised the manuscript critically for important intellectual content. All authors read and approved the final manuscript.

Funding

None.

Availability of data and materials

Publicly available datasets were analyzed in this study. These data can be found here: GDC TCGA database (https://portal.gdc.cancer.gov/); GEO datasets (https://www.ncbi.nlm.nih.gov/geo/): GSE65904.

Ethics approval and consent to participate

This study was based on two publicly available databases so we did not require Institutional Review Board approval. This study was conducted in concordance with the Helsinki Declaration.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

All authors are from the Plastic Surgery Department of the First Affiliated Hospital.

Yan Y, Leontovich AA, Gerdes MJ, Desai K, Dong J, Sood A, et al. Understanding heterogeneous tumor microenvironment in metastatic melanoma. PLOS ONE. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2019 2019-1-20;14(6):e216485.
Gershenwald JE, Guy GP. Stemming the Rising Incidence of Melanoma: Calling Prevention to Action. JNCI: Journal of the National Cancer Institute. 2016 2016-1-1;108(1).
Carr S, Smith C, Wernberg J. Epidemiology and Risk Factors of Melanoma. SURG CLIN N AM. 2020 2020-1-1;100(1):1-12.
Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA-CANCER J CLIN. [Journal Article]. 2022 2022-1-1;72(1):7-33.
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer Journal for Clinicians. 2018;68(6):394-424.
Thomas NE, Edmiston SN, Alexander A, Groben PA, Parrish E, Kricker A, et al. Association Between NRAS and BRAF Mutational Status and Melanoma-Specific Survival Among Patients With Higher-Risk Primary Melanoma. JAMA ONCOL. [Journal Article; Multicenter Study; Research Support, N.I.H., Extramural]. 2015 2015-6-1;1(3):359-68.
Arneth B. Tumor Microenvironment. MEDICINA-LITHUANIA. [Journal Article; Review]. 2019 2019-12-30;56(1).
Junttila MR, de Sauvage FJ. Influence of tumour micro-environment heterogeneity on therapeutic response. NATURE. [Journal Article; Review]. 2013 2013-9-19;501(7467):346-54.
Sahai E, Astsaturov I, Cukierman E, DeNardo DG, Egeblad M, Evans RM, et al. A framework for advancing our understanding of cancer-associated fibroblasts. NAT REV CANCER. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Research Support, U.S. Gov't, Non-P.H.S.; Review]. 2020 2020-3-1;20(3):174-86.
Zhou L, Yang K, Andl T, Wickett RR, Zhang Y. Perspective of Targeting Cancer-Associated Fibroblasts in Melanoma. J CANCER. [Journal Article; Review]. 2015 2015-1-20;6(8):717-26.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC BIOINFORMATICS. [Journal Article; Research Support, N.I.H., Extramural]. 2008 2008-12-29;9:559.
Liu B, Chen X, Zhan Y, Wu B, Pan S. Identification of a Gene Signature for Renal Cell Carcinoma-Associated Fibroblasts Mediating Cancer Progression and Affecting Prognosis. FRONT CELL DEV BIOL. [Journal Article]. 2020 2020-1-20;8:604627.
Liu B, Zhan Y, Chen X, Hu X, Wu B, Pan S. Weighted gene co-expression network analysis can sort cancer-associated fibroblast-specific markers promoting bladder cancer progression. J CELL PHYSIOL. [Journal Article; Research Support, Non-U.S. Gov't]. 2021 2021-2-1;236(2):1321-31.
Zheng H, Liu H, Li H, Dou W, Wang X. Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblast Signature for Predicting Prognosis and Therapeutic Responses in Gastric Cancer. FRONT MOL BIOSCI. [Journal Article]. 2021 2021-1-20;8:744677.
Cabrita R, Lauss M, Sanna A, Donia M, Skaarup LM, Mitra S, et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. NATURE. [Journal Article; Research Support, Non-U.S. Gov't]. 2020 2020-1-1;577(7791):561-5.
Cirenajwis H, Ekedahl H, Lauss M, Harbst K, Carneiro A, Enoksson J, et al. Molecular stratification of metastatic melanoma using gene expression profiling: Prediction of survival outcome and benefit from molecular targeted therapy. Oncotarget. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2015 2015-5-20;6(14):12297-309.
Racle J, de Jonge K, Baumgaertner P, Speiser DE, Gfeller D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. ELIFE. [Journal Article; Validation Study]. 2017 2017-11-13;6.
Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. GENOME BIOL. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2017 2017-11-15;18(1):220.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. GENOME BIOL. [Journal Article]. 2016 2016-10-20;17(1):218.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. NAT MED. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2018 2018-10-1;24(10):1550-8.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. NAT COMMUN. [Journal Article; Research Support, N.I.H., Extramural]. 2013 2013-1-20;4:2612.
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). [Journal Article]. 2021 2021-8-28;2(3):100141.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. [Journal Article; Research Support, Non-U.S. Gov't]. 2012 2012-5-1;16(5):284-7.
Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J STAT SOFTW. [Journal Article]. 2010 2010-1-20;33(1):1-22.
Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC BIOINFORMATICS. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2013 2013-1-16;14:7.
Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. CELL SYST. [Journal Article]. 2015 2015-12-23;1(6):417-25.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. NAT METHODS. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Research Support, U.S. Gov't, Non-P.H.S.]. 2015 2015-5-1;12(5):453-7.
Ghandi M, Huang FW, Jane-Valbuena J, Kryukov GV, Lo CC, McDonald ER, et al. Next-generation characterization of the Cancer Cell Line Encyclopedia. NATURE. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2019 2019-5-1;569(7757):503-8.
Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue-based map of the human proteome. SCIENCE. [Journal Article; Research Support, Non-U.S. Gov't]. 2015 2015-1-23;347(6220):1260419.
Flach EH, Rebecca VW, Herlyn M, Smalley KS, Anderson AR. Fibroblasts contribute to melanoma tumor growth and drug resistance. MOL PHARMACEUT. [Journal Article; Research Support, N.I.H., Extramural]. 2011 2011-12-5;8(6):2039-49.
Jobe NP, Rosel D, Dvorankova B, Kodet O, Lacina L, Mateu R, et al. Simultaneous blocking of IL-6 and IL-8 is sufficient to fully inhibit CAF-induced human melanoma cell invasiveness. HISTOCHEM CELL BIOL. [Journal Article]. 2016 2016-8-1;146(2):205-17.
Dvorankova B, Szabo P, Kodet O, Strnad H, Kolar M, Lacina L, et al. Intercellular crosstalk in human malignant melanoma. PROTOPLASMA. [Journal Article; Review]. 2017 2017-5-1;254(3):1143-50.
Li Z, Zhang J, Zhou J, Lu L, Wang H, Zhang G, et al. Nodal Facilitates Differentiation of Fibroblasts to Cancer-Associated Fibroblasts that Support Tumor Growth in Melanoma and Colorectal Cancer. CELLS-BASEL. [Journal Article; Research Support, Non-U.S. Gov't]. 2019 2019-6-4;8(6).
Colombo M, Raposo G, Thery C. Biogenesis, secretion, and intercellular interactions of exosomes and other extracellular vesicles. ANNU REV CELL DEV BI. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2014 2014-1-20;30:255-89.
Chen Y, McAndrews KM, Kalluri R. Clinical and therapeutic relevance of cancer-associated fibroblasts. NAT REV CLIN ONCOL. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Review]. 2021 2021-12-1;18(12):792-804.
Wu F, Yang J, Liu J, Wang Y, Mu J, Zeng Q, et al. Signaling pathways in cancer-associated fibroblasts and targeted therapy for cancer. SIGNAL TRANSDUCT TAR. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2021 2021-6-10;6(1):218.
Hosseini-Alghaderi S, Baron M. Notch3 in Development, Health and Disease. BIOMOLECULES. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2020 2020-3-23;10(3).
Bodas M, Subramaniyan B, Moore AR, Metcalf JP, Ocanas SR, Freeman WM, et al. The NOTCH3 Downstream Target HEYL Is Required for Efficient Human Airway Basal Cell Differentiation. CELLS-BASEL. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't]. 2021 2021-11-18;10(11).
Aburjania Z, Jang S, Whitt J, Jaskula-Stzul R, Chen H, Rose JB. The Role of Notch3 in Cancer. ONCOLOGIST. [Journal Article; Review]. 2018 2018-8-1;23(8):900-11.
van der Rest M, Garrone R. Collagen family of proteins. FASEB J. [Journal Article; Review]. 1991 1991-10-1;5(13):2814-23.
Ricard-Blum S. The collagen family. CSH PERSPECT BIOL. [Journal Article; Review]. 2011 2011-1-1;3(1):a4978.
Gilkes DM, Semenza GL, Wirtz D. Hypoxia and the extracellular matrix: drivers of tumour metastasis. NAT REV CANCER. [Review]. 2014 2014-6-1;14(6):430-9.
Chen D, Qin Y, Dai M, Li L, Liu H, Zhou Y, et al. BGN and COL11A1 Regulatory Network Analysis in Colorectal Cancer (CRC) Reveals That BGN Influences CRC Cell Biological Functions and Interacts with miR-6828-5p. CANCER MANAG RES. [Journal Article]. 2020 2020-1-20;12:13051-69.
Yousefi H, Vatanmakanian M, Mahdiannasser M, Mashouri L, Alahari NV, Monjezi MR, et al. Understanding the role of integrins in breast cancer invasion, metastasis, angiogenesis, and drug resistance. ONCOGENE. [Journal Article; Review]. 2021 2021-2-1;40(6):1043-63.
Lai JP, Sandhu DS, Shire AM, Roberts LR. The tumor suppressor function of human sulfatase 1 (SULF1) in carcinogenesis. J GASTROINTEST CANC. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Review]. 2008 2008-1-20;39(1-4):149-58.
Zhou B, Lin W, Long Y, Yang Y, Zhang H, Wu K, et al. Notch signaling pathway: architecture, disease, and therapeutics. SIGNAL TRANSDUCT TAR. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2022 2022-3-24;7(1):95.
Li J, Liu G, Li L, Yao Z, Huang J. Research progress on the effect of autophagy-lysosomal pathway on tumor drug resistance. EXP CELL RES. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2020 2020-4-15;389(2):111925.
Patel NH, Sohal SS, Manjili MH, Harrell JC, Gewirtz DA. The Roles of Autophagy and Senescence in the Tumor Cell Response to Radiation. RADIAT RES. [Journal Article; Research Support, N.I.H., Extramural; Research Support, Non-U.S. Gov't; Research Support, U.S. Gov't, Non-P.H.S.; Review]. 2020 2020-8-1;194(2):103-15.
Zamame RJ, Romagnoli GG, Falasco BF, Gorgulho CM, Sanzochi FC, Dos SD, et al. Blocking drug-induced autophagy with chloroquine in HCT-116 colon cancer cells enhances DC maturation and T cell responses induced by tumor cell lysate. INT IMMUNOPHARMACOL. [Journal Article]. 2020 2020-7-1;84:106495.
Vijay K. Toll-like receptors in immunity and inflammatory diseases: Past, present, and future. INT IMMUNOPHARMACOL. [Journal Article; Review]. 2018 2018-6-1;59:391-412.
Akesolo O, Buey B, Beltran-Visiedo M, Giraldos D, Marzo I, Latorre E. Toll-like receptors: New targets for multiple myeloma treatment? BIOCHEM PHARMACOL. [Journal Article; Research Support, Non-U.S. Gov't; Review]. 2022 2022-5-1;199:114992.
Wang Y, Chen L, Yu M, Fang Y, Qian K, Wang G, et al. Immune-related signature predicts the prognosis and immunotherapy benefit in bladder cancer. CANCER MED-US. [Journal Article; Research Support, Non-U.S. Gov't]. 2020 2020-10-1;9(20):7729-41.
Pan Y, Yu Y, Wang X, Zhang T. Tumor-Associated Macrophages in Tumor Immunity. FRONT IMMUNOL. [Journal Article; Review]. 2020 2020-1-20;11:583084.
Komi D, Redegeld FA. Role of Mast Cells in Shaping the Tumor Microenvironment. CLIN REV ALLERG IMMU. [Journal Article; Review]. 2020 2020-6-1;58(3):313-25.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Weighted Gene Co-expression Network Analysis Identifies a Cancer-Associated Fibroblasts Signature for Predicting Prognosis and Immunotherapeutic Responses in Skin Cutaneous Melanoma

Status:

Version 1

Abstract

Figures

1. Background

2. Materials and Methods

2.1 Data Acquisition and Methodology

2.2 Infiltration estimation of CAFs and analysis of tumor immune microenvironment

2.3 Weighted gene co-expression network analysis

2.4 Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Analyses

2.5 Prediction model construction and validation based on CAFs-related hub genes

2.6 Enrichment analysis of CAFs signature genes

2.7 Immune infiltration and immune checkpoint analysis

2.8 Immunotherapy response prediction

2.9 Validated by Cancer Cell Line Encyclopedia (CCLE) and Human Protein Atlas (HPA) databases

2.10 Statistical Analysis

3. Results

3.1 Patients with SKCM with higher infiltration of CAFs and lower stromal scores had worse OS

3.2 Co-expression network of CAFs and stromal scores

3.3 Enrichment analysis of hub genes

3.4 6 CAFs signature genes used to construct predictive models

3.5 CAFs signature genes are positively associated with CAFs infiltration and CAFs markers

3.6 GSEA and ssGSEA analysis of 6 CAFs signature genes

3.7 Correlation analysis of RS with tumor immune microenvironment and immune checkpoints

3.8 Analysis of the correlation between RS and immunotherapy

3.9 Multidimensional validation of CAFs signature genes in CCLE and HPA databases

4. Discussion

5. Conclusion

Abbreviations

Declarations

References

Additional Declarations

Status:

Version 1