Clinical outcome of melanoma patients is impacted by the immune-related tumor microenvironment
In order to characterize human metastatic melanoma samples according to their associated immune phenotypes, we inferred the tumor microenvironment (TME)-associated immune cell populations [10] of 164 metastatic samples from The Cancer Genome Atlas (TCGA) database (Additional file 1: Table S1). This allowed the identification of three major groups containing at least 20 samples each: G1 (pink) with 37 samples, G2 (green) with 57 samples, and G3 (orange) with 65 samples (Fig. 1a). Bootstrapping showed that G2 and G3 were highly stable clusters, with a clusterwise Jaccard bootstrap mean of 0.56 for G1, 0.75 for G2, and 0.84 for G3. Five samples did not cluster within any of the groups and were excluded from the subsequent analysis (white labels, Fig. 1a).
According to the TME predictions, G1 samples are enriched in naïve, memory and plasma B cells, and depleted in resting natural killer (NK) cells; G2 samples are enriched in CD8+ T cells, monocytes, and macrophages M1; and G3 samples are enriched in M0 macrophages and depleted in plasma cells, CD8+ T cells, memory activated CD4+ T cells, follicular T helper cells, activated NK cells, monocytes, and resting dendritic cells (p ≤ 0.05, Mann–Whitney test—MW) (Fig. 1a and Additional file 1: Figure S1A). Regarding the histopathological characteristics, G1 samples present a higher percentage of tumor cells and, conversely, a small percentage of stromal cells when compared to G2 and G3 samples (Additional file 1: Figure S1B). Compared to G1 samples, samples in G3 have more necrotic cells and samples in G2 have more lymphocyte infiltration (Additional file 1: Figure S1B). Importantly, all samples selected for the analysis had > 50% of tumor cells.
The TCGA consortium has previously classified melanomas in four molecular subtypes based on their genomic profiles (BRAF, RAS, NF1, and triple wild-type), three transcriptional subclasses based on their gene expression signatures (Immune, Keratin, MITF-low) and four methylation clusters based on their methylation patterns (CpG island-methylated, hyper-methylated, hypo-methylated and normal-like) [36]. We did not find any association between the immune-related groups described here and the molecular subtypes (Fisher’s exact test p = 0.4) or the methylation profiles identified (Fisher’s exact test p = 0.12) (Fig. 1a and Additional file 1: Table S1) [36]. However, regarding the transcriptomic classification of melanomas, we found that G3 was enriched with metastatic melanoma samples classified into the MITF-low subclass whereas melanomas clustered into the Immune subclass were more present among G1 and G2 (Fisher’s exact test p = 2.51e−06).
To confirm the identified TME profiles of each group, we validated our findings using an alternative gene signature-based approach [11] (Additional file 1: Figure S1C). We performed non-supervised clustering of the TME populations using the median of the population signature predicted for each group. Five different B-cell types, CD4+ T naive cell type and non-immune cell types such as microvascular endothelial cells showed a strong association with G1. G2 showed a strong signature of all T cell subpopulations, such as CD8+, regulatory T (Treg), Th1 and Th2 cells along with macrophages and basophils. Finally, G3 showed a strong signature of non-immune cell types such as smooth muscle cells, pericytes, neurons, and mesenchymal stem cells. Although this alternative approach does not cover the M0 macrophage signature, the presence of a weaker signature of other immune cell types in G3 and the overlap between the immune cell populations predicted by the two different approaches present in G1 and G2 confirm our initial classification into three distinct immune-related groups.
To investigate whether different TME compositions correlate with clinical outcome, we next assessed the overall survival of melanoma patients grouped according to our TME-based classification (Fig. 1b). Patients clustered in G2 (green line) showed significantly better overall survival than G3 patients (orange line) (p = 0.01, log-rank test, Hazard Ratio (HR) = 0.49, confidence interval (CI) .95 = 0.28–0.85). G1 patients’ overall survival did not differ from the other two groups. To better account for the impact of specific immune cell types on survival, we calculated the Hazard Ratio (HR) of the relative fractions of each cell type population divided into quartiles (Fig. 1c). We identified poorer HR related to the presence of memory B cells (p = 0.002, log-rank test, HR = 2.54, CI.95 = 1.41–4.6) and M0 macrophages (p < 0.001, log-rank test, HR = 1.78, CI.95 = 1.31–2.4), which are both immune cell types highly represented in G1 and G3, respectively. In addition, we analyzed the correlation among the different predicted immune cell types in all samples combined (Additional file 1: Figure S2A). This revealed several connections between cells that could potentially impact tumor growth and help explain the observed prognosis. Examples include a positive correlation between Treg cells and memory B cells or M0 macrophages (r = 0.47 and r = 0.56, respectively, p ≤ 0.05, Pearson correlation test—PCT) and a strong negative correlation between CD8+ T cells and M0 macrophages (r = − 0.41, p ≤ 0.05, PCT). Most cell types were positively correlated, of which the correlation between neutrophils and activated mast cells was the highest observed (r = 0.65, p ≤ 0.05, PCT), followed by CD4+ memory resting T cells and M2 macrophages (r = 0.62, p ≤ 0.05, PCT), and CD8+ T cells and CD4+ memory activated T cells (r = 0.61, p ≤ 0.05, PCT). Since we observed a significant difference between the overall survival of G2 and G3 patients, we next investigated which immune cell types better discriminate these groups by calculating the point-biserial correlation coefficient (Additional file 1: Figure S2B). We observed that monocytes were better correlated to G2 (r = 0.29, p ≤ 0.05, PCT), and M0 macrophages to G3 samples (r = − 0.27, p ≤ 0.05, PCT). Also, a strong and significant negative correlation between these two cell populations was observed (r = − 0.59, p ≤ 0.05, PCT).
Finally, we evaluated the frequency and diversity of T and B-cell receptors (TCR and BCR, respectively) as they correspond, at least in part, to the T and B lymphocyte populations observed. G1 and G2 presented higher numbers of unique alpha/beta chains clonotypes in the TCR repertoire (medians of 28.5 and 17, respectively) in comparison to G3 (median of 4 and p ≤ 0.0001 in both comparisons using MW test), but no significant difference was found between G1 and G2 (Fig. 1d). G2 also presented a greater variety of gamma chain clonotypes (median of 2) when compared to either G1 (median of 1; p = 0.011) or G3 (median of 0; p ≤ 4.9e10−6). The lymphocyte repertoire diversity, assessed by the inverse Simpson diversity index, was higher in G1 and G2 (p ≤ 9e10−8 for G1 and G2) compared to G3 (Fig. 1e). Higher TCR diversity in G1 and G2 supports the notion that more diverse antigen collections are being recognized and contribute to explain the better outcome of G2 patients. We also analyzed the BCR repertoire and found a higher immunoglobulin A (IgA) representation in G1 (median of 15 clones per sample; p ≤ 7.1e10−6) compared to G2 and G3 (median of 7 and 8, respectively) (Fig. 1f). The same was also true for IgG clonotypes which were more represented in G1 (median of 35) in comparison to G2 (p = 0.002, median of 30) and G3 (p = 0.0008, median of 30). G3 presented higher IgM representation (median of 8) in comparison to G2 (median of 5, p = 0.027). The BCR diversity was also higher in G1 compared to G2 and G3 (p ≤ 5.3e10−8) (Fig. 1g). In accordance with the variety of lymphocyte subpopulations, the expression of chemokines that recruit several of the G2-enriched immune cell populations was augmented in this group when compared to G3 (adjusted p-value ≤ 0.001), except for CCL5, which acts on neutrophils and endothelial cells (Fig. 1h). This lack of chemokine expression in G3 can explain, at least in part, the paucity of CD8+ T lymphocytes in these samples, a characteristic that underpins malignant development.
Phenotypic characteristics that favor immune evasion co-occur in melanomas with worse prognosis
Besides avoiding the recruitment of inflammatory cells and lymphocytes to the TME, tumor cells can develop strategies to become “invisible” to the immune system or to inhibit effector responses. These strategies, based on distinct mechanisms of immune evasion, may have a significant impact on clinical outcome of cancer patients. Therefore, we characterized different aspects of known immune evasion mechanisms in the TME-grouped melanoma samples.
First, we identified signatures of effector cells (EC), suppressor cells (SC), checkpoint molecules (CP) and major histocompatibility complex (MHC) expression in the three groups, which were previously shown to predict response to immune checkpoint blockade [19]. Weaker signatures of SC and CP were detected in G3 whereas a stronger signature of EC was observed in G1 and G2 (Fig. 2a). These results reinforce that ongoing effector responses are more likely to occur in G1 and G2 and that suppressor mechanisms are of great importance in G3 and could favor immune evasion.
Tumors with a higher mutation burden tend to generate more neoantigens and, thus, be more immunogenic [37]. We observed significantly higher numbers of somatic mutations per megabase (Mb) in G2 and G3 compared to G1 (p = 0.0069 and 0.00013, respectively, MW test, Fig. 2b). The median of exonic mutations per Mb was 60 for G1, 130 for G2, and 187 for G3 (Additional file 1: Figure S3A) alike the intronic mutations (Additional file 1: Figure S3B). Similarly, G2 and G3 samples presented more neoepitopes compared to G1, of which the median was 51 for G1, 77 for G2 and 120 for G3 (p = 0.026 and p = 0.0037, respectively; MW test, Fig. 2c). Despite observing no difference in neoepitope and mutation burden between G2 and G3, and based on the lower T and B cell diversity observed in G3 (Fig. 1e, g), we searched for features that could impact neoepitope processing and presentation and thus diminish the recognition of a tumor cell. First, we searched for high impact mutations in genes belonging to the antigen processing and presentation pathway (Additional file 1: Table S2) and identified an increase in neoepitope burden according to the number of genes mutated in the pathway in a given sample (Fig. 2d and Additional file 1: Table S3). Two samples from G1 and two from G2 had one mutated gene on this pathway while nine samples from G3 presented at least one gene mutated. Within G3, samples bearing mutations on this pathway had increased neoepitope burden when compared to the non-mutated samples (p = 0.022, MW test). No statistical test could be performed for G1 and G2 samples due to the small number of events.
Another feature that impacts neoepitopes presentation is the predicted affinity between them and the Human Leukocyte Antigen I (HLA-I) alleles so that higher dissociation constants (KD) account for lower stability of the neoepitope:HLA-I complex and therefore a less efficient presentation. The median of neoepitope’s KD was calculated for each sample and compared among groups (Fig. 2e). Neoepitopes’s KD were higher in G3 compared to G2 (median of 149.5 nM and 126.9 nM, respectively; p = 0.042, MW test), supporting the idea that the immune evasion phenotype in G3 samples may also involve selecting less stable neoepitope:HLA-I complexes.
Mutational signatures are consequences of different etiological agents that, in the process of mutational carcinogenesis, favor some specific DNA transversions and transitions [34] and these mutational signatures have been related to patient’s survival [38]. Therefore, we asked if the distinct mutational signatures were differentially represented across the immune-related groups. Two main mutational signatures were identified across the groups: signature 7 (green) and 1A (pink), both corresponding to C>T transitions (Additional file 1: Figure S3C). Signature 7 is associated to ultraviolet radiation and previously described as related to melanoma, while signature 1A is broadly present among tumors [34]. We observed a similar mutational signature profile across the immune-related groups, suggesting that they can not be distinguished by genomic signatures.
Differential miRNA/mRNA expression profiles dictate the clinical outcome of melanoma patients
In order to better understand the complex gene regulatory network related to the differences observed in patients’ survival and the immune evasion profiles of tumors, we compared the gene expression profiles of G2 and G3 samples. We identified 1783 differentially expressed genes (DEGs) (adjusted p ≤ 0.001) (Additional file 1: Figure S4A, B) and 93 differentially expressed miRNAs (DEM) (adjusted p ≤ 0.05) (Additional file 1: Figure S4C, D), of which 641 DEG and 34 DEM were up- and 1142 DEG and 59 DEM were downregulated in G3 (worse prognosis). The top 20 DEGs comprise four upregulated (AL035610.1, ST8SIA5, NRXN1 and FAM131B) and 16 downregulated (IFNG, TMEM155, RP11-109E24.1, CD8A, KLRK1, AC104820.2, FASLG, CCL4, GZMA, RP11-1094M14.8, AKAP5, RP11-1094M14.5, CLEC2D, TRGC2, CTC-303L1.1 and JAKMIP1) genes (Additional file 1: Table S4), while within the top 20 DEMs, 12 were upregulated (mir-206, mir-203a, mir-183, mir-205, mir-6892, mir-675, mir-887, mir-200c, mir-375, mir-1-1, mir-1-2, mir-130b) and 8 were downregulated (mir-142, mir-7702, mir-342, mir-4494, mir-155, mir-4491, mir-150, mir-6842) (Additional file 1: Table S5).
Gene Set Enrichment analysis (GSEA) revealed that the pathways enriched in G3 were involved in epithelial cell processes such as Keratinization and Extracellular Matrix Organization (False Discovery Rate (FDR) < 0.001, test = Permutation test for both) (Additional file 1: Figure S4E). The enriched pathways in G2 participate in immune response processes such as Downstream TCR signaling (FDR = 0.002, Permutation test), Interferon Signaling (FDR < 0.001, Permutation test), and Adaptive Immune System (FDR < 0.001, Permutation test). To analyze the potential impact of miRNA expression on the presence of each predicted cell type, analysis of the correlation of the DEMs and the predicted immune cell types was performed. A significantly negative correlation was found between miRNAs downregulated in G3 and M0 and M2 macrophages, including mir-148a: macrophage M2 (r = − 0.4, p = 8.81e−06, PCT) and mir-29c: macrophage M0 correlations (r = − 0.48, p = 7.68e−08, PCT) (Fig. 3a). Moreover, CD8+ T cell type, enriched in G2, was positively correlated to mir-142 and mir-7702 (r = 0.55 for both, p = 2.14e−10 and p = 1.93e−10, respectively, PCT) as well as to other 11 miRNAs. Dendritic resting cell type, decreased in G3, also correlated to mir-203a (r = 0.43, p = 1.57e−06, PCT) and mir-205 (r = 0.51, p = 7.16e−09, PCT).
In order to illustrate the relationship between miRNAs and the differences in gene expression in G2 versus G3, we determined miRNA-target gene (MTG) pairs. We identified 1111 DEGs predicted as targets of the DEMs (Fig. 3b), with 139 MTG pairs being negatively correlated, suggesting differences in the canonical regulation of mRNAs by miRNAs between the groups. We then sought to identify the interactions between melanoma and TME cells that were most likely to occur. For that, we used a manually curated list of mRNAs and miRNAs putatively expressed by melanoma cells retrieved from the Cancer Cell Line Encyclopedia (CCLE) and a publicly available list of circulating miRNA [29] and integrated this information into a gene expression network that shows the expression pattern of each gene in G3 relative to G2 (Fig. 3c). A hundred and thirteen MTG pairs (81%) involved putative circulating miRNA produced by melanoma cells suggesting that the DEMs are regulating the gene expression both at the intra- and intercellular levels. We also evaluated the impact of the expression of each MTG pair on patients’ outcome. In total, we found 74 MTG pairs with a significant impact on overall survival (p ≤ 0.05, log-rank test, Additional file 1: Figure S5). Of those, 12 pairs (mir-1296/GCH1, mir-1306/ZBP1, mir-1306/ICAM3, mir-1306/SEPT1, mir-142/NR2F6, mir-142/TRIM28, mir-142/ZNF74, mir-1914/HLA-F, mir-323a/PIM2, mir-3619/CD74, mir-4736/ZNF74 and mir-342/PPM1F) did not involve potentially circulating miRNAs.
miRNA-based intercellular communication in the TME impacts patient’s survival
Analysis of bulk tumor samples allowed us to identify potential miRNA-target genes but not to predict which cells within the tumor were expressing the miRNA or the target gene. To refine these analyzes, we took advantage of two publicly available datasets. First, we analyzed the candidate miRNA-target gene expression at a cellular level using Single Cell RNA-Seq (SCRS) data from melanoma tumors [14]. Second, we used a dataset consisting of miRNA expression profiles in extracellular vesicles isolated from the plasma of metastatic melanoma patients (hereafter referred to as bearer patients), healthy individuals, and R0-operated patients whose melanomas were surgically removed (with clear margin). The latter was subdivided into high relapse-risk and low relapse-risk based on tumor staging [21]. When combining these two datasets, we identified circulating miRNAs that are potentially expressed by TME cells and whose targets are mainly expressed in tumor cells, as well as miRNA potentially expressed by tumor cells that may affect gene expression in cells from the microenvironment. We further investigated the impact of these miRNAs in overall survival across different types of malignancy (Additional file 2). We describe below some examples to depict the potential role of circulating miRNAs in the interplay between the TME and malignant cells.
Circulating miRNAs potentially suppressing tumor growth
One of the strongest MTG pairs identified in our analysis with a significant impact on patients’ prognosis was mir-150/HILPDA (Fig. 4a, b). The putative circulating mir-150 was found to be downregulated in G3. Consistently, its potential target, the Hypoxia-inducible lipid droplet-associated (HILPDA) gene, which encodes a protein related to intracellular lipid accumulation, was found to be upregulated in the same group (Additional file 2). analysis of SCRS data revealed that HILPDA was mainly expressed in malignant cells (Fig. 4c). The average expression of the mature form of hsa-miR-150-5p was lower in the extracellular vesicles from bearer patients (Log Fold Change (LFC): − 1.48, p = 0.001, MW test) when compared to healthy individuals (Fig. 4d). Interestingly, comparison of bearer and high-risk patients demonstrated a gain in hsa-miR-150-5p expression after surgical removal of the melanoma (LFC: 0.62, p = 0.03, MW test), although expression was still lower than in healthy individuals (LFC: − 0.62, p = 0.03). Higher levels of mir-150 were also related to better survival in lung adenocarcinoma and ovarian cancer (p = 0.011, log-rank test, HR = 0.67, CI.95 = 0.49–0.91, and p = 0.021, log-rank test, HR = 0.72, CI.95 = 0.53–0.95 respectively. Additional file 2). Altogether, these findings suggest that the release of hsa-miR-150-5p in EVs can be involved in the crosstalk between immune and tumor cells associated with suppression of tumor growth.
With similar characteristics, mir-342 was found to be downregulated in the group with worse prognosis, and its lower expression was associated with a worse overall survival in melanoma (p = 0.002, log-rank test, HR = 0.43, CI.95 = 0.25–0.75, Additional file 1: Figure S6) and breast carcinoma (p = 9.9e−04, log-rank test, HR = 0.51, CI.95 = 0.33–0.76, Additional file 2). Consistently, its target, the protein-coding gene PPM1F, of which higher expression was associated with a worse outcome (p = 0.033, log-rank test, HR = 1.77, CI.95 = 0.25–0.75, Additional file 1: Figure S5), was found to be upregulated in G3 (Additional file 2) and expressed in several cell types according to the SCRS analysis, including malignant cells (Additional file 1: Figure S7). Although not described as circulating miRNA in our initial analysis, hsa-miR-342-3p was downregulated in EVs from bearer patients when compared to healthy individuals (LFC = − 0.80, p = 0.001, MW test). Similar to what was observed for hsa-miR-150-5p, an increase in the expression of hsa-miR-342-3p was detected in the plasma of both high- and low-risk melanoma survivors when compared to metastatic melanoma patients (LFC = 0.61 and p = 0.004, LFC = 0.53 and p = 0.03, respectively, MW test, Additional file 1: Figure S8).
Circulating miRNAs potentially favoring tumor growth
We also identified miRNAs whose high expression were associated with worse outcome across at least three different tumor types, such as the putative circulating mir-130b, which was up-regulated in the group with worse prognosis (LFC = 0.65, padj = 1.75e−03). Seven out of 15 mir-130b targets identified among the DEGs were mainly expressed in non-malignant cells, such as T and B cells. (Additional file 2). Our independent analysis on the circulating miRNA dataset showed that average expression of miR-130b-3p was higher in the plasma EVs of high-risk patients when compared to plasma samples from healthy individuals (LFC = 1.69, p = 0.03, MW test, Additional file 1: Figure S8) also indicating the putative tumor-promoting role of this miRNA.
The putative circulating mir-149, which potentially regulates 39 targets of the 139 MTGs (Fig. 3c and Additional file 2), was upregulated in G3. Although the expression levels of mir-149 was not significantly associated with survival (Fig. 5a), all of its targets were downregulated in G3 and their low expression levels were associated with poor overall survival (Fig. 5b and Additional file 1: Figure S5). According to the melanoma SCRS data, all mir-149 targets identified were expressed in lymphocytes, some of which were exclusively expressed in this cell type and include CD96, CD48, SLAMF7, FASLG, NUGGC (Additional file 1: Figure S7) and NLRC5 (Fig. 5c). Interestingly, NLRC5 encodes a transcription coactivator of genes involved in HLA class I presentation, such as TAP1, B2M and HLA-A/B/C, all of which were downregulated in G3 (Fig. 5d), suggesting an immune evasion phenotype in G3 regulated by mir-149. It is worth to note that, among these genes, only HLA-A/B could not be associated with poor prognosis (Additional file 1: Figure S5). Unfortunately, the restricted expression of hsa-miR-149-3p in only two bearer patients combined with the absence of other mature forms of miR-149 in the microarray chip prevented us from assessing the differential expression of this miRNA in the EV dataset.
Another candidate MTG pair with a potential impact in the melanoma-TME crosstalk identified by our analysis was mir-1914/HLA-F (Fig. 5e, f). mir-1914 was upregulated in G3 (LFC = 1.12, padj = 1.17e−02) and was found to regulate four putative targets (Additional file 2). Low expression of one of them, HLA-F, was associated with a poor outcome (p = 0.027, log-rank test, HR = 0.55, CI.95 = 0.32–0.94) (Fig. 5f). HLA-F codes for a non-classical HLA molecule that forms a heterodimer with B2M [39] and, although its expression was detected in several cell types, it was mainly enriched in TME cells including CAFs, NK and T cells (Fig. 5g). Moreover, high expression of hsa-miR-1914-3p was observed in EVs derived from melanoma bearer and high relapse risk patients (LFC = 2.63 and 1.54, p = 0.004 and 0.045, respectively, Wilcox test, Fig. 5h) when compared to healthy individuals. High levels of mir-1914 were also associated with worse survival in lung adenocarcinoma (p = 0.046, log-rank test, HR = 1.41, CI.95 = 1.004–1.96, Additional file 2).