Next Article in Journal
Tumor-Associated Microglia and Macrophages in the Glioblastoma Microenvironment and Their Implications for Therapy
Next Article in Special Issue
Whole-Genome Analysis of De Novo Somatic Point Mutations Reveals Novel Mutational Biomarkers in Pancreatic Cancer
Previous Article in Journal
The Synergistic Anti-Tumor Activity of EZH2 Inhibitor SHR2554 and HDAC Inhibitor Chidamide through ORC1 Reduction of DNA Replication Process in Diffuse Large B Cell Lymphoma
Previous Article in Special Issue
Identification of Cancer Hub Gene Signatures Associated with Immune-Suppressive Tumor Microenvironment and Ovatodiolide as a Potential Cancer Immunotherapeutic Agent
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cancers
Volume 13
Issue 17
10.3390/cancers13174250
cancers-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptional Reprogramming and Constitutive PD-L1 Expression in Melanoma Are Associated with Dedifferentiation and Activation of Interferon and Tumour Necrosis Factor Signalling Pathways

by
Antonio Ahn
1,
Euan J. Rodger
1,2,
Jyoti Motwani
1,
Gregory Gimenez
1,
Peter A. Stockwell
1,
Matthew Parry
3,
Peter Hersey
4,
Aniruddha Chatterjee
1,2,*,† and
Michael R. Eccles
1,2,*,†
1
Department of Pathology, Dunedin School of Medicine, University of Otago, 270 Great King Street, Dunedin 9054, New Zealand
2
Maurice Wilkins Centre for Molecular Biodiscovery, Level 2, 3A Symonds Street, Auckland 1010, New Zealand
3
Department of Mathematics & Statistics, University of Otago, 710 Cumberland Street, Dunedin 9054, New Zealand
4
Melanoma Immunology and Oncology Group, The Centenary Institute, University of Sydney, Camperdown, NSW 2050, Australia
*
Authors to whom correspondence should be addressed.
†
Joint senior authors.
Cancers 2021, 13(17), 4250; https://0-doi-org.brum.beds.ac.uk/10.3390/cancers13174250
Submission received: 15 June 2021 / Revised: 7 August 2021 / Accepted: 13 August 2021 / Published: 24 August 2021
(This article belongs to the Special Issue Bioinformatics, Big Data and Cancer)
Browse Figures
Versions Notes

Abstract

:

Simple Summary

Melanoma, an aggressive form of skin cancer, is frequently associated with drug resistance in the advanced stages. For instance, frequently resistance is observed in sequential treatment of melanoma with targeted therapy and immunotherapy. In this research, the authors investigated whether potential transcriptional mechanisms and pathways associated with PD-L1 protein expression could underlie targeted therapy drug resistance in melanoma. The authors found a PD-L1 expression transcriptional pattern underlies resistance to targeted therapy in a subgroup of melanomas. These melanomas were markedly dedifferentiated, as compared to melanomas that were not drug resistant. Understanding changes in transcription and molecular pathways that lead to drug resistance could allow researchers to develop interventions to prevent drug resistance from occurring in melanoma, which could also be relevant to other cancer types.

Abstract

Melanoma is the most aggressive type of skin cancer, with increasing incidence worldwide. Advances in targeted therapy and immunotherapy have improved the survival of melanoma patients experiencing recurrent disease, but unfortunately treatment resistance frequently reduces patient survival. Resistance to targeted therapy is associated with transcriptomic changes and has also been shown to be accompanied by increased expression of programmed death ligand 1 (PD-L1), a potent inhibitor of immune response. Intrinsic upregulation of PD-L1 is associated with genome-wide DNA hypomethylation and widespread alterations in gene expression in melanoma cell lines. However, an in-depth analysis of the transcriptomic landscape of melanoma cells with intrinsically upregulated PD-L1 expression is lacking. To determine the transcriptomic landscape of intrinsically upregulated PD-L1 expression in melanoma, we investigated transcriptomes in melanomas with constitutive versus inducible PD-L1 expression (referred to as PD-L1CON and PD-L1IND). RNA-Seq analysis was performed on seven PD-L1CON melanoma cell lines and ten melanoma cell lines with low inducible PD-L1IND expression. We observed that PD-L1CON melanoma cells had a reprogrammed transcriptome with a characteristic pattern of dedifferentiated gene expression, together with active interferon (IFN) and tumour necrosis factor (TNF) signalling pathways. Furthermore, we identified key transcription factors that were also differentially expressed in PD-L1CON versus PD-L1IND melanoma cell lines. Overall, our studies describe transcriptomic reprogramming of melanomas with PD-L1CON expression.

Graphical Abstract

1. Introduction

Melanoma is the most deadly form skin cancer, as it frequently presents with highly aggressive features, including high propensity to metastasise and innate drug resistance. Moreover, these features frequently occur at a relatively early stage in the growth of the tumour [1]. Treatment of metastatic melanoma has been revolutionized over the last decade, as greater understanding has emerged of two critical hallmarks of melanoma. Firstly, a large proportion of melanomas (50–65%) are addicted to MAPK signalling through BRAF or NRAS mutations. In keeping with this, inhibition of the oncogenic BRAF protein has resulted in significant response rates in BRAF mutant melanomas. Secondly, irrespective of mutation status, melanoma is frequently dependent on immune suppression through programmed death 1 (PD1) signalling upon the binding of ligand, either PD-L1 or PD-L2 [2]. Anti-PD1 therapy, which inhibits binding of PD-L1 or PD-L2 to the PD1 receptor, reactivates immune responses and has greatly improved melanoma patient survival [3]. Unfortunately, for both of these therapies, resistance inevitably develops, and currently no robust biomarker has been identified that is able to predict patient response. Neither PD1 nor PD-L1/PD-L2 expression accurately predict response, and the basis for intrinsic resistance to anti-PD1 treatment of melanoma is incompletely understood. Nevertheless, PD-L1 is constitutively expressed in some melanomas despite the absence of immune cell infiltration in the tumour [4,5,6]. Furthermore, PD-L1 has also been shown to be upregulated upon development of resistance to MAPK pathway inhibitors, and it is accompanied by tumour cell-intrinsic transcriptomic reprogramming [4,5]. However, it is currently unclear how intrinsic transcriptomic reprogramming occurs.
As we have previously described [6,7] PD-L1 expression in melanoma can broadly be categorised into mechanisms that are mediated by the presence or absence of tumour infiltrating lymphocytes (TILs). In the presence of TILs, tumour-associated PD-L1 expression is predominantly induced via interferons (IFN) [8,9,10] and/or cytokines, such as tumour necrosis factor (TNF) [11] secreted from TILs. We refer to melanoma cells with inducible PD-L1 expression as PD-L1IND. In contrast, PD-L1 expression without TILs is largely mediated cell-intrinsically (or constitutively) via genetic or epigenetic mechanisms, and we refer to these types of melanoma cells as PD-L1CON [2,12].
In this study, we describe transcriptomic features of melanoma cell lines with constitutive high expression of PD-L1 (PD-L1CON) compared to melanoma cell lines with low inducible levels of PD-L1 expression (PD-L1IND). We found that PD-L1CON expression was associated with a reprogrammed transcriptomic state, inclusive of dedifferentiation, and active innate inflammatory pathways associated with IFN and TNF signalling and reduced oxidative phosphorylation. Overall, changes in the expression of key transcription factors were observed that drive dedifferentiation (such as the loss of MITF and SOX10), as well as key transcription factors that enhance IFN and TNF signalling and PD-L1 expression (such as IRF1, JUN, and FOSL2, in which the latter two encode an AP-1 protein complex). We additionally found that the altered expression of transcription factors and transcriptional reprogramming correlated with the PD-L1CON mRNA expression, and with factors associated with resistance to MAPK pathway inhibitors.

2. Materials and Methods

2.1. Selection and Culture of Melanoma Cell Lines

We analysed transcriptome feature of seven PD-L1CON lines, which includes four (CM143.pre, CM143.post, NZM9, NZM40) lines from our previous study [6] and three additional new PD-L1CON cell lines (MM127, MM595, and COLO239F) that were included in this study. These three cell lines were kindly provided by professor Glen Boyle from the QIMR Berghofer Medical Research Institute. For PD-L1IND group, we have analysed 10 cell lines, which includes six (CM138, CM150.post, CM145.pre, CM145.post, NZM22, NZM42) lines from our previous study [6] and four new cell lines (NZM12, NZM15, WM115, and WM2664). The CM138, CM145.pre, CM145.post, CM150.post, CM143.pre, and CM143.post lines were cultured in DMEM medium (Invitrogen) supplemented with 10% foetal bovine serum (FBS) and 1% penicillin–streptomycin. WM115 and WM2664 were cultured in Minimum Essential Medium (MEM-α) (Invitrogen) supplemented with 1% penicillin–streptomycin (Gibco, NY, USA) and 10% FBS. NZM9, NZM40, NZM12, NZM15, and NZM42 were cultured in MEM-α media supplemented with 1% penicillin-streptomycin, 5% FBS, and 0.1% Insulin-transferrin-selenium (Roche) [13]. MM127, MM595, and COLO239F were grown in RPMI 1640 (Thermo Fisher Scientific, Waltham, MA, USA) Medium supplemented with 10% FBS and 1% penicillin–streptomycin. All cells were grown under standard cell culture conditions (5% CO2, 21% O2, 37 °C, humidified atmosphere), except WM115, which were cultured at 35 °C.

2.2. Flow Cytometry Analysis

Fixable viability stain 450 (FVS450, BD horizon, catalog#: 562247, clone: 29E.2A3) was used to stain dead cells in order to selectively analyse live cells. The FVS450 stain has a fluorescence emission maximum at 450 nm. The PE anti-human CD274/PD-L1 (Biolegend, catalog# 329706) antibody and the isotype control antibody (PE Mouse IgG2b, Biolegend, catalog#:400314) have a maximum excitation at 575 nm. No overlap in fluorescence emission was detected between the FVS450 and the anti-PDL1 fluorophore or isotype control antibodies. The PD-L1 expression levels of melanoma cell lines were determined using BD FACS CantoII. All analyses were performed using the Kaluza (Beckman Coulter, version 2.0, Carlsbad, CA, USA) software. Approximately 10,000 events/cells were measured for each sample. Gating strategy was used to exclude dead cells and doublet cells. The median fluorescence intensity (MFI) for the isotype control and anti-PDL1 was obtained. The MFI for PD-L1 staining was normalised for background absorbance by subtracting out the isotype fluorescence value. We used an arbitrary MFI cut-off value of <500 for the PD-L1IND group and >10,000 for PD-L1CON group. Gene expression values confirmed that there were significantly higher levels of CD274 expression in the PD-L1CON group, with an increased log2 fold change of 7.2. For IFN-γ induction of PD-L1, between 50,000 to 100,000 cells were seeded in a single well of a 24-well plate overnight with 1 mL of media. For each sample, around 10 wells were seeded in order to obtain a total amount of between 500,000 to 1 million cells. The following day, the media was removed and fresh media with IFN-γ (final concentration of 100 ng/mL, prospec, catalog#: CYT-206) was added to the cells. After 1 day of IFN-γ induction, flow cytometry was used to assess PD-L1 expression.

2.3. RNA Extraction and Reverse Transcription

Total RNA was isolated from melanoma cell lines using the RNeasy Mini Kit (Qiagen, catalog#:74106) following the protocol manual. This involved cell lysis, homogenisation of the lysate using the QIAshredder (Qiagen, catalog#:79656), and using a spin column to selectively purify RNA. DNase (RNase-Free DNase Set, Qiagen, catalog#:79254) was used to degrade DNA during the extraction as outlined in the RNeasy Mini Handbook. Quality control was first performed on the Nanodrop 2000 Spectrophotometer (Thermo Fisher Scientific) to assess the RNA purity using the ratio of absorbance at 260 to 280 nm higher than 1.8. The RNA integrity was assessed using the Agilent Bioanalyzer 2100 with the RNA integrity number (RIN) higher than nine. Reverse transcription from RNA to complementary DNA (cDNA) was performed using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, catalog#:4368814).

2.4. Additional Cell Line Cohorts and TCGA Data for Melanoma

Melanoma cell line gene expression data was obtained from external sources to validate our results. The getGEO function from the GEOquery package was used to download four gene expression datasets with a total of 175 melanoma cell lines. The datasets included: GSE7127 (n = 63) [14], GSE4843 (n = 45) [15], GSE61544 (n = 13) [16], GSE80829 (n = 54) [17]. Normalisation was performed where the unprocessed microarray data or raw count RNA-seq matrix was available. For GSE7127 (Affymetrix U133 Plus 2 microarray platform), the data was normalised using Robust Multichip Average (RMA, from the affy package). For GSE61544, the raw count data was downloaded and normalised using TMM (edgeR package) [18]. For GSE4843 the normalized data (MAS5.0 normalization) and GSE80829 (FPKM using conditional quantile normalization) data was downloaded used for analysis. RNA-seq data (FPKM) of melanoma cell lines prior to and following acquired resistance upon treatment with MAPK inhibitors as well as patient tumours that were on MAPKi treatment (4 weeks) was obtained from GSE75313 [5]. Patient tumour RNA-seq data containing baseline, and after tumour progression, was obtained from GSE65186 (FPKM normalised gene expression matrix) (n = 70) [19]. RNA-seq data (FPKM normalised gene expression matrix) of melanoma cell lines (n = 8), treated with IFN-γ for 2 to 5 weeks or TNF for 3 days, was acquired from GSE152755 [20].

2.5. TCGA SKCM Data Download

The RNA-seq count matrix was downloaded from the harmonized TCGA SKCM dataset (GRCh38). This was done using the RTCGAbiolinks package in R which downloads the data from the Genomic Data Commons (GDC) database [21]. The GDCquery function was used with the parameters: project = “TCGA-SKCM”, data.category = “Transcriptome Profiling”, data.type = “Gene Expression Quantification”, workflow.type = “HTSeq-Counts”). The TMM method was used to normalise the count matrix.

2.6. Generation and Processing of Transcriptome Data

RNA-seq processing for PD-L1IND and PD-L1CON melanoma cell lines was performed with poly-A-tail selection, paired-end reads, read length of 2 × 100 bps, and a total of 40 million reads. Adaptor trimming was done using cleanadaptors from the DMAP package [22,23]. The RNA-Seq design and analysis that was used in this work has been described in detail previously [24]. Pseudoalignment methods such as Kallisto and Salmon were found to outperform other alignment methods when measuring lncRNA expression abundance [25]. Therefore, we mapped our reads to the hg38 reference genome using Kallisto (version 0.44.0) [26]. Each sample was run with 100 bootstraps and with the bias argument to correct for potential sequence-based bias. Annotations were acquired from GENECODE (Release 28 GRCh38.p12), which entailed nucleotide sequences of all transcripts (protein-coding and lncRNA transcripts) on the reference chromosomes. Tximport was used to import the kallisto gene-level counts data into R [27].

2.7. Statistical Analysis

Following filtering our genes with low counts (lower than 1 in at least seven samples), the TMM method was used to normalise for sequencing depth and RNA composition bias using the edgeR package [18]. Differential expression analysis was performed using edgeR quasi-likelihood method [28]. A False Discovery Rate (FDR) adjusted p value threshold of 0.05 was used to call significant genes. Gene Set Enrichment Analysis (GSEA) was performed using gene sets available in the Broad Institute Molecular Signatures Database (MsigDB) which included H1 (hallmarK), C2 (curated), and C5 (gene ontology) [29]. CAMERA test was performed as available from the edgeR package and a FDR adjusted 5 × 10−5 p value was used as the statistically significant threshold [30]. For generating a gene-set score for each sample, single sample GSEA (ssGSEA) [31] was used from the GSVA Bioconductor R package [32]. The IFN score, TNF score, differentiation score, and oxidative phosphorylation score were obtained from MSigDB from the “MOSERLE_IFNA_RESPONSE”, “PHONG_TNF_TARGETS_UP”, “GO_MELANOCYTE_DIFFERENTIATION”, and “KEGG_OXIDATIVE_PHOSPHORYLATION” gene sets, respectively. The viral mimicry score was self-curated from published articles [33,34] and included DDX58, DDX41, IFIH1, OASL, IRF7, IRF1, ISG15, MAVS, IFI27, IFI44, IFI44L, and IFI16. Differentiation signature genes from Tsoi and colleagues [17] (from Table S3) were acquired to further evaluate the dedifferentiation signature in our melanoma cell lines. These genes included 6 groups, in order, from least differentiated to the most differentiated: (1) Undifferentiated, (2) Undifferentiated-Neural crest-like, (3) Neural crest-like, (4) Transitory, (5) Transitory-Melanocytic, (6) Melanocytic. Z-scores were generated for each gene and unsupervised hierarchical clustering was performed. To infer cytotoxic immune activity from gene expression data, the CYT-score was calculated by finding the geometric mean from GZMA (granzyme A) and PRF1 (perforin) expression values [35]. The absolute abundance of eight immune and two stromal cell populations were estimated from normalised RNA-seq data using MCPcounter [36]. R codes to calculate the moving average and to generate figures were acquired from Riesenberg and colleagues [37].

3. Results

3.1. PD-L1CON Melanoma Cell Lines (High PD-L1 Group) Have a Distinct Gene Expression Profile Compared to PD-L1IND Melanoma Cell Lines (Low PD-L1 Group)

As we have previously suggested [6], melanomas can be categorized into four subgroups based on high or low levels of tumour-associated PD-L1 protein expression, and the presence or absence of TILs in melanoma tissue. In this study, we defined these subgroups as TIL+/PDL1+ (group 1), TIL-/PDL1- (group 2), TIL+/PDL1- (group 3), and TIL-/PDL1+ (group 4) (Figure 1A). Factors leading to PD-L1 expression in melanoma tumour tissues are complicated by the fact that cytokines secreted from TILs frequently induce PD-L1 expression in tumour cells, which masks whether cytokine-independent, tumour cell-intrinsic PD-L1 expression (i.e., PD-L1 constitutive expression, or PD-L1CON) has occurred. Melanoma cell lines grown in vitro do not have TILs present, and yet some melanoma cell lines constitutively express high PD-L1 levels, while other melanoma cell lines express low baseline levels of PD-L1, but the latter subgroup can be induced to express higher PD-L1 levels with IFNγ stimulation. Therefore, we used melanoma cell lines to begin our investigation, as they lack TILs, and high levels of constitutive PD-L1 expression in these cell lines consequently represent tumour cell-intrinsic mechanisms. Assignment of a panel of melanoma cell lines to the high and low PD-L1 expression groups was accomplished by characterizing cell surface PD-L1 protein levels and mRNA expression levels, which were assessed using flow cytometry and RNA-Seq [38] (Figure S1A,B, Table S1). We found that, as opposed to a continuum from low to high of CD274 (PD-L1 mRNA) expression, the CD274 expression in the melanoma cell lines was either at a low level (yet inducible, hence PDL1IND), or at a constitutively high level (PDL1CON). The CD274 mRNA levels (as quantified by RNA-Seq) strongly correlated with the PD-L1 surface protein expression levels (Pearson R = 0.92, p value = 1.4 × 10−7 Figure S1C) in the analysed 17 cell lines (seven PD-L1CON cell lines and ten PD-L1IND cell lines). Given that melanoma cell lines and patient tumours occasionally exhibit dysfunctional IFN signalling due to genetic defects [39,40,41,42,43], we additionally verified, by PD-L1 upregulation upon IFNγ treatment, that ten low PD-L1IND cell lines in our panel were inducible, which confirmed them as being able to respond to IFNγ induction (Figure S1D), while PD-L1CON weren’t inducible [38].
Principal component analysis of the top 500 protein coding (Figure 1B) and long non-coding RNA (Figure 1C) genes with the highest variance segregated the cell lines into two distinct clusters. Unsupervised hierarchical clustering showed the same result with a clear segregation of the two groups (Figure 1D,E), which suggests there is a distinct transcriptional state between PD-L1CON and PD-L1IND melanoma subgroups. Differential expression analysis identified 462 upregulated and 298 downregulated protein coding genes in the PD-L1CON group, compared to the PD-L1IND group (FDR adjusted p value < 0.05) (Table S2). As expected, CD274 expression was significantly increased in the PD-L1CON group (FDR adjusted p value = 0.001, log2 fold-change = 7.1, Figure S2A,C). PDCD1LG2 (which encodes the PD-L2 protein, and is also located in the same chromosome location as CD274, i.e., 9p24.1) was also significantly upregulated in the PD-L1CON group (FDR adjusted p value = 0.03, log2FC = 5.0, Figure S2A,D). For the long non-coding RNA (lncRNA) analysis, we identified 106 upregulated and 71 downregulated lncRNA in PD-L1CON lines compared to the PD-L1IND lines (Figure S2B, Table S2). The upregulated lncRNAs included 75 lincRNAs, 26 antisense, 3 bidirectional promoter and 2 sense overlapping lncRNAs, while the downregulated lncRNAs included 37 lincRNAs, 31 antisense, and 3 sense intronic lncRNAs. Next, given that lncRNAs can regulate gene expression via cis mechanisms, we assessed all differentially expressed lncRNAs with loci mapped adjacent to (within 100 kb) a differentially expressed protein-coding gene, which identified 54 mRNA–lncRNA pairs. 51 out of the 54 (94.4%) mRNA–lncRNA pairs showed expression changes in the same direction (either up or down), with 33 pairs being upregulated (top-right quadrant, Figure 1F), while 18 pairs showing downregulation (bottom-left quadrant, Figure 1F in the PD-L1CON lines compared to the PD-L1IND group. These results suggest coordinated expression changes for the mRNA and lncRNA expression. Overall, PD-L1CON cell lines had a distinct expression profile for protein-coding genes and lncRNAs, and moreover the vast majority of the differentially expressed mRNAs and lncRNAs located in genomic proximity to each other were altered in the same direction.

3.2. PD-L1CON Cell Lines Exhibit a Distinct Transcriptome That Represents a State of Dedifferentiation, Enhanced IFN, and TNF Signalling Pathways and Reduced Oxidative Phosphorylation

To determine which biological processes are altered in the PD-L1CON cell lines compared to PD-L1IND cell lines, we performed gene-set enrichment analysis using CAMERA [30] and C2, C5, and H gene set collections from the Molecular Signature Database (MSigDB) [29]. Thirty-seven gene sets were found to be significantly altered (FDR adjusted p value threshold of 5 × 10−5). Genes involved in interferon-alpha (IFN-α) and tumour necrosis factor (TNF) signalling were among the most significantly upregulated processes in the PD-L1CON group ((FDR adjusted p value = 3.0 × 10−10 (MOSERLE_IFNA_RESPONSE) and 1.1 × 10−9 (PHONG_TNF_TARGETS_UP), Figure 2A and Figure S3A). Given that the PD-L1CON cell lines do not harbour immune cells, this suggests that the IFN and TNF signalling could be activated via the expression of dsRNA derived from endogenous retroviral elements (ERVs). Expression of ERVs can trigger dsRNA sensors and a downstream signalling cascade of interferon response genes (also referred to as a viral mimicry response) [33,34]. Genes involved as dsRNA sensors, as well as the viral mimicry response genes, were also highly upregulated in PD-L1CON group (FDR adjusted p value = 7.5 × 10−4 and 2.1 × 10−4, respectively) (Figure S3A). These findings suggest that PD-L1CON samples may have enhanced activation of viral mimicry pathways. This is also supported in our previous findings, in that ERV genomic regions were hypomethylated in PD-L1CON melanoma cells [6], and were a characteristic feature associated with PD-L1 expression. The downregulated genes in the PD-L1CON samples were significantly enriched for genes involved in melanocyte differentiation (FDR adjusted p value = 2.7 × 10−6, Figure 2A and Figure S3A), suggestive of a dedifferentiated state. Recently, Tsoi and colleagues found that melanomas can exhibit transcriptomic states that are coupled to a differentiation trajectory, which consists of four progressive steps [17], corresponding to (1) undifferentiated, (2) neural crest-like, (3) transitory, and (4) melanocytic. Unsupervised hierarchical clustering of the cell lines based on these melanoma differentiation genes confirmed that the PD-L1CON cell lines were enriched for genes characteristic of the undifferentiated state (Figure 2B). The PD-L1IND cell lines were further clustered into two subgroups from their transcriptomic profiles, with one group corresponding to a neural crest-like gene signature, and the other group to a melanocytic state. Moreover, there was a downregulation of genes involved in oxidative phosphorylation in the PD-L1CON group, suggesting metabolic reprogramming had occurred (FDR adjusted p value = 7.2 × 10−8 (KEGG_OXIDATIVE_PHOSPHORYLATION)) (Figure 2B and Figure S3A). Given that dedifferentiation is closely associated with enhanced invasiveness we investigated whether the PD-L1CON samples were enriched for an invasive gene expression profile. Using two independent invasive gene expression gene sets [44,45], we found PD-L1CON cell lines had higher levels of expression of genes involved in melanoma invasion (Figure S4A,B) although three PD-L1IND cell lines (NZM22, WM115, CM145.pre) had weak expression levels of invasive genes, whereas CM138 clustered more closely to the PD-L1CON cell lines.

3.3. Validation of Upregulated IFN and TNF Pathways, and Downregulated Differentiation and Oxidative Phosphorylation Pathways in Association with Constitutive CD274 Expression in Melanoma Cell Lines

To validate the gene expression signature of an upregulated IFN and TNF pathway, and downregulation of differentiation and oxidative phosphorylation genes in PD- L1CON melanoma cells, we utilised four external gene expression datasets containing a total of 175 melanoma cell lines. A score was generated for the 462 upregulated and 298 downregulated genes in the PD-L1CON samples (referred to as the up and down PD-L1CON score respectively) as well as the genes involved in IFN signalling, TNF signalling, differentiation, and oxidative phosphorylation using single sample gene set enrichment analysis (ssGSEA) [31,32]. In all four external gene expression datasets, CD274 expression was significantly positively correlated with IFN and TNF scores and negatively correlated with the differentiation score (Figure 2C). However, the oxidative phosphorylation score did not show consistent correlation with CD274 expression in the external gene expression datasets (Figure S3B).
We have previously found that blocking IFN signalling in the PD-L1CON samples did not alter the PD-L1CON expression in two of the PD-L1CON samples that were also used in this study [38], which demonstrates that PD-L1CON expression is associated with a cell intrinsic regulated IFN response [6]. However, in tumours, immune cell infiltrates and their secreted cytokines are thought to be mainly responsible for extrinsically inducing IFN signalling, which in turn upregulates PD-L1 expression [9]. To understand whether immune cell independent PD-L1CON expression is similar to that in tumours with an immune cell dependent PD-L1 expression, we asked whether the PD-L1CON samples share an expression signature with those tumours with an elevated immune infiltration. To this end, the TCGA SKCM dataset which contains RNA-seq data for 472 tumours was analysed. A large proportion of significantly upregulated transcripts (135 out of 462 or 29%) in the PD-L1CON samples were positively correlated (Pearson correlation value of higher than 0.25) with the CYTscore (Figure 2D), which is a well-established index that reflects cytolytic immune activity [35]. Moreover, the upregulated and downregulated PD-L1CON genes were positively and negatively correlated, respectively, with the absolute abundance of a large range of immune subtypes as estimated by MCPcounter [36] (Figure 2D), demonstrating that PD-L1CON melanoma cell lines have an immune cell independent active cytokine signalling pathway similar to melanoma tumours where immune infiltrates stimulate this pathway.
IFN gamma (IFN-γ) and TNF stimulation have been shown to induce both PD-L1 expression and a dedifferentiation phenotype [20,46,47]. To investigate whether IFN-γ and TNF treatment can also activate the PD-L1CON gene expression signature in cell lines, we analysed an independent RNA-seq dataset (GSE152755) [20] where eight melanoma cell lines were treated with either IFN-γ for 2 to 5 weeks, or TNF for 3 days. Both IFN-γ and TNF treatment demonstrated gene expression changes towards the PD-L1CON signature (Figure S5A,B). It is important to note that the cytokine induced signatures in these cells are reversible upon removal of stimulation [46,47], whereas in PD-L1CON samples, these signatures are not dependent on any external stimulation.

3.4. Validation of Upregulated IFN and TNF Pathways, and Downregulated Differentiation and Oxidative Phosphorylation Pathways in Association with Constitutive CD274 Expression in Melanoma Tumour Tissues

To validate the association of PD-L1CON expression with active IFN and TNF signalling, and reduced differentiation and oxidative phosphorylation expression signatures in melanoma tumours, we used the SKCM TCGA dataset consisting of 458 melanoma tumour samples. As mentioned above, CD274 expression in melanoma tissues is typically stimulated by immune infiltrates rather than an intrinsic mechanism, therefore we used the significantly differentially expressed genes that were found in the PD-L1CON cell lines as a surrogate for PD-L1CON expression. A large proportion of the upregulated protein-coding genes (n = 462) found in the PD-L1CON samples were positively correlated (Pearson r > 0.25) with CD274 expression (136 out of 462 or 29%), IFN score (136 out of 462 or 29%), and TNF score (213 out of 462 or 46%). Furthermore, the upregulated protein-coding genes were negatively correlated (lower than −0.25 Pearson correlation value) with the differentiation (137 out of 462 or 30%) and oxidative phosphorylation scores (151 out of 462 or 33%) (Figure 3A and Figure S3C). In contrast, a large proportion of the downregulated protein-coding genes (n = 298) were positively correlated (Pearson r > −0.25) with the differentiation-score (127 out of 298 or 43%) and the oxidative phosphorylation score (108 out of 298 or 37%) (Figure 3A and Figure S3C). The same pattern was also seen for the lncRNAs (Figure 3B). Collectively, these results support the finding that PD-L1CON expression is associated with an enhanced IFN and TNF signalling, and reduced differentiation and oxidative phosphorylation.

3.5. Lineage Specific, TNF, and IFN Associated Transcription Factors Are Differentially Expressed in the PD-L1CON Samples

Key transcription factors (TFs) can play a role in modulating the transcriptomic program. Therefore, we investigated what type of TFs are altered in expression in the PD-L1CON cell lines. To this end, we analysed 1107 TF mRNAs with known DNA binding motifs [48] and identified 25 significantly upregulated TF mRNAs (out of 462 significantly upregulated genes) and 19 significantly downregulated TF mRNAs (out of 298 significantly downregulated genes) in the PD-L1CON cell lines (Figure 4A). To determine whether the expression patterns of these TF mRNAs could be generalised to melanoma, we analysed the four external gene expression datasets and the TCGA SKCM RNA-seq dataset in the TF context. The 25 upregulated TF and 19 downregulated TFs were predominantly positively and negatively correlated with CD274 expression in all five datasets (Figure 4B). Furthermore, we identified POU2F2, IRF1, and FOSL2 as the most highly correlated TFs with CD274 expression in melanoma across all five datasets (Figure 4C). IRF1 and FOSL2 have known roles in upregulating PD-L1 expression [8,49,50]. In contrast, SOX10, RXRG, and SOX5 had the lowest correlation with CD274 expression amongst all TFs investigated (Figure 4C). SOX10 and MITF are key transcription factors that regulate the development of melanocytes. SOX10 expression showed a near complete loss (log2FC = −10.7, FDR adjusted p value = 0.0003) and MITF expression was significantly reduced in expression (log2FC= −3.8, FDR adjusted p value = 0.04) in the PD-L1CON cell lines (Figure S6E,F). The upregulation of these TF expression further supports our observation that the PD-L1CON melanoma subgroup is associated with substantial transcriptomic reprogramming.
Additionally, we identified all differentially expressed lncRNAs whose loci were located in close proximity/adjacent to differentially expressed TF mRNAs (within 100 kb of each other). Correspondingly, eight TF mRNA and lncRNA differentially expressed pairs were found (Figure 4D). Out of these eight TFs, four were identified to be closely associated with CD274 expression in the five external datasets. This included AL031587.1 antisense lncRNA near SOX10 (distance of 44,743 bp), the AC116366.2 intergenic lncRNA near IRF1 (distance of 54,895 bp), and FLJ31356 antisense lncRNA near FOSL2 (distance of 0 bp, as they overlap) (Figure 4D).

3.6. The PD-L1CON Expression Signature Is Associated with Transriptomic Reprogramming, and Correlates with MAPK Inhibitor Resistance in Melanoma Cell Lines

Key transcriptional changes that were shown to drive acquired resistance include increased expression of cMET, reduced expression of LEF1, and an enrichment of YAP1 signature genes [19]. PD-L1CON cell lines had increased expression of cMET (log2FC = 2.6, FDR adjusted p value = 0.07), reduced expression of LEF1 (log2FC = −3.2, FDR adjusted p value = 0.017), and enrichment of YAP1 signature genes (t.test p value = 0.012) (Figure S6B–D). Other gene expression changes driving drug resistance include reduced levels of MITF and SOX10, which were both downregulated in the PD-L1CON cells. Moreover, consistent with increased expression of Receptor Tyrosine Kinases (RTK), which are associated with drug resistance, we found a large number of RTK genes were highly expressed in the PD-L1CON cell lines including EGFR (log2FC= 4.2, FDR adjusted p value = 0.054), PDGFRB (log2FC= 3.6, FDR adjusted p value = 0.03), and AXL (log2FC= 2.5, FDR adjusted p value = 0.34) (Figure S6A,G–I). These data support the notion that PD-L1CON cell lines have a transcriptional profile corresponding to a drug resistant phenotype.
CD274 expression was found to be cell-intrinsically increased upon development of resistance to MAPK inhibitors, when accompanied by transcriptomic reprogramming, but not mediated by BRAF mutations (splicing or amplification) reactivating the MAPK pathway [5]. We utilised this dataset to assess our PD-L1CON gene expression signature in the context of development of resistance. In cell lines where resistance co-occurred with transcriptomic reprogramming, there was an increase of CD274 expression as well as IFN, TNF, and “up” PD-L1CON scores (Figure 5B), and a corresponding reduction in differentiation, oxidative phosphorylation, MPAS, and “down” PD-L1CON scores. Moreover, we found that TFs associated with PD-L1CON expression, and involved in IFN signalling and TNF signalling (IRF1, JUN, and FOSL2), were all upregulated (Figure 5E), except for IRF4, which was downregulated in these cell lines. Consistently, IRF4 was the only TF involved in IFN and TNF signalling that was also downregulated in PD-L1CON cell lines (Figure 5C). In addition, TFs involved in melanocyte differentiation (including SOX10, MITF, and LEF1) were downregulated in cell lines that the authors had associated with the development of drug resistance (Figure 5E). In contrast, in cell lines where resistance was mediated by reactivation of the MAPK signalling pathway (involving BRAF splicing or amplification), only minimal changes in the expression of CD274 and the PD-L1CON related biological processes (including IFN signalling, TNF signalling, differentiation, and oxidative phosphorylation) were observed (Figure 5A). Additionally, the expression of TFs involved in IFN or TNF signalling, or melanocyte differentiation were unaltered when the drug resistance was mediated by BRAF splicing or amplification mutations (Figure 5D).

3.7. Transcriptomic Changes, including PD-L1CON Expression, Occur at Defined Stages during the Development of Drug Resistance

To investigate the timing of when PD-L1 expression changes and biological processes may occur during drug treatment, we investigated RNA-Seq data derived from two melanoma cell lines (M229 and M238) [5], which had been treated with an MAPK inhibitor, and for which RNA-Seq analysis was carried out at different time points between the initial resistance, generation of drug tolerance, and finally permanent acquired resistance. After two days of BRAFi treatment (BRAFi2D), the time points included; at the Drug Tolerant Persisters (DTP) stage where small subpopulation of persisting cells remained; at the Drug Tolerant Proliferative Persisters (DTPP) stage where proliferation was regained; and finally at the Single Drug Resistant (SDR) stage, which is a permanent resistant state to BRAFi (months to years of drug treatment). There was little or no PD-L1CON “up” score after 2 days of BRAFi treatment (Figure 5F), suggesting that the PD-L1CON signature is not up-regulated due to immediate on-target effects of treatment with BRAFi. CD274 expression was slightly decreased at the DTP stage, however increased at the proliferative phase (DTPP) and was further increased at the SDR phase (Figure 5F). TNF, IFN, and PD-L1CON “up” scores were largely increased at the DTPP stage, whereas dedifferentiation and PD-L1CON “down” scores were largely decreased at the SDR stage. This suggested that innate IFN and TNF signalling responses precede dedifferentiation. Further analysis revealed that JUN and FOSL2, which encode proteins that are part of the AP-1 complex, and are part of the TNF signalling pathway, were increased relatively early after two days of BRAFi treatment and were maintained until acquired resistance developed, whereas IRF1, which can drive PD-L1 expression, increased at the DTPP stage when CD274 was also overexpressed (Figure 5G). Furthermore, the dedifferentiation TFs, including MITF and LEF1, were downregulated early after two days of BRAFi treatment, whereas SOX10 and MEF2C were reduced upon the development of acquired resistance (SDR), which was also when the differentiation score was greatly reduced. This suggested that SOX10 and MEF2C downregulation are required or are the main contributors to the dedifferentiation gene expression signature. Overall, these data support the notion that PD-L1CON expression occurs at the proliferative stage of drug resistance and that this is mediated by a transcriptome reprogramming, and is accompanied by increases in TNF and IFN signalling along with the IRF1 expression. Moreover, these expression changes further increase when resistance is stabilised, at which point the differentiation expression signature is downregulated.

3.8. The PD-L1CON Expression Signature Is Associated with Transcriptomic Reprogramming of Melanomas following MAPK Inhibitor Resistance in Patients

We next investigated whether the PD-L1CON expression signature occurs in melanomas from patients who develop acquired resistance to MAPK inhibitor treatment. To address this, we analysed RNA-Seq data from matched pre- and post-treatment tumour biopsies, following progression on treatment [19]. We compared melanomas with BRAF mutational resistance mechanisms (BRAF splicing and amplification) to melanomas lacking known mutations that drive resistance. In melanomas where no mutational resistance mechanism was identified, we observed the same direction of up and down gene expression changes to the PD-L1CON cell lines, and to the transcriptomic changes of Song and colleagues’ dataset, where resistance was mediated by transcriptomic reprogramming. This consisted of elevated levels of CD274 expression, and elevated IFN, TNF, and the “up” PD-L1CON scores (normalised to matched pre-treatment tumours from the same patient), while oxidative phosphorylation, down PD-L1CON, and differentiation scores were reduced (Figure S7A). In addition, there were lower levels of the differentiation related TFs and higher levels of expression of JUN, IRF1, JUNB, and FOSL2 observed in melanomas with no BRAF mutational resistance mechanisms (normalised to patient matched pre-treatment tumours), compared to melanomas containing BRAF mutational resistance mechanisms (Figure S7B). These results support the notion that the PD-L1CON expression signature accompanies transcriptomic reprogramming and treatment resistance in melanomas lacking a BRAF mutational resistance mechanism.

4. Discussion

The use of targeted therapies based on MAPK pathway inhibition, or of anti-PD1 immunotherapy, has improved melanoma patient survival, although drug resistance greatly limits long-term survival for most patients with advanced melanoma. A large number of mutational and non-mutational (transcriptomic or epigenetic) resistance mechanisms have been identified in melanoma [5,19,51,52]. Moreover, these have been associated with a cell-intrinsic upregulation of the immunosuppressive PD-L1 protein [4,5]. In the present study, to explore how PD-L1 expression is associated with transcriptomic reprogramming and drug resistance, we have assessed the transcriptomic landscape of melanoma cell lines with constitutively high levels of PD-L1 expression (PD-L1CON). Overall, our study has shown that PD-L1CON cells have many similarities to melanoma cell lines and patient tumours that exhibit intrinsically upregulated PD-L1 expression following the development of resistance to MAPK pathway inhibitors via non-mutational and/or epigenetic mechanisms. Our research reveals that this similarity is inclusive of a gene expression profile corresponding to dedifferentiation, and active cytokine (TNF and IFN) signalling pathways. Dedifferentiation is a non-mutational mechanism of drug resistance, which melanomas commonly exploit. There did not appear to be an association of transcriptomic reprogramming with the patient’s sex, or with the tumour stage, although there were a limited number of samples in our analysis. Consistent with other studies showing that treatment of melanoma and melanocytes cells with TNF-α and IFN-γ cytokines can induce a dedifferentiated state, we have found that the dedifferentiated state of PD-L1CON cells is associated with constitutively active IFN and TNF signalling pathways [17,37,46,53]. Furthermore, given that the PD-L1CON and PD-L1IND cell line groups used in this study were classified into these groups independent of BRAF or NRAS mutation status (four cell lines were BRAF mutant, two NRAS mutant cell lines, and one BRAF/NRAS wild-type), these data suggest that the presence of the BRAF or NRAS genetic mutations did not influence the occurrence of a PD-L1CON transcriptomic expression pattern, and lend support to the notion that transcriptional or epigenetic factors can over-ride mutation-driven resistance in the absence of a targeted drug treatment.
It is important to note that almost all of the PD-L1CON melanoma cell lines (six out of seven) in this study had not been previously exposed to targeted inhibitors either in the patient (in vivo) or in vitro. Nevertheless, each PD-L1CON cell line contained the “resistant” gene expression signature. This suggests that stress factors, other than targeted inhibitor drugs, may have induced the resistant transcriptional profile. This could include stresses, such as hypoxia, or treatment using chemotherapeutic drugs in the patients before the cell lines were initiated. The clinical relevance of these findings is that, in some patients, PD-L1CON cells may occupy a subpopulation of cells in the tumour prior to targeted inhibitor therapy, which could then result in innate resistance. Moreover, one of the PD-L1CON melanoma cell lines, which was previously exposed to MAPK pathway inhibitors (i.e., CM143-post), had a related cell line from the same patient, which also exhibited PD-L1CON expression prior to treatment exposure (CM143-pre), suggesting that, in this patient, a PD-L1CON transcriptomic pattern may have represented a pre-existing subpopulation of cells that was present prior to treatment, and which gave rise to the drug-resistant phenotype.

5. Conclusions

Our study offers new insights into transcriptomic patterns associated with PD-L1CON melanoma cell lines, the expression of relevant TFs, and reprogramming of melanoma cells towards a drug resistant phenotype. For instance, SOX10 and MITF are master TF regulators of melanocyte differentiation, whereby their downregulation can promote dedifferentiation and confer resistance to MAPK pathway inhibitors. This then results in activation of oncogenic signalling pathways, providing alternatives to the MAPK pathway, such as pathways driven by RTKs [16,52,54,55,56,57]. Both SOX10 and MITF were significantly downregulated in the PD-L1CON cells, along with a dedifferentiation gene expression profile and increased expression of RTK genes which highlights the importance of PD-L1CON expression during adaptive drug resistance. Moreover, TFs that were significantly upregulated consisted of AP-1 (JUN and FOSL2), which have recently been shown to not only mediate the TNF signalling pathways, but also to alter the transcriptional and enhancer chromatin landscape by serving as pioneer TFs [58,59,60]. Transcriptomic reprogramming and differentiation often involves pioneer TFs, which facilitate and alter access to distinct regulatory elements, and consequently alter the chromatin landscape [61]. Moreover, SOX9 has also been shown to be a pioneer TF that mediates stem cell plasticity [62,63,64], and we found this was significantly upregulated in the PD-L1CON cells. Finally, this study also provides foundations for further research to investigate mechanistic pathways (such as chromatin remodeling) leading to transcriptional reprogramming and the involvement of gene expression signatures, such as PD-L1CON, in cancer drug-resistant phenotypes.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/cancers13174250/s1, Table S1: Information on the mutation status and the PD-L1IND and PD-L1CON group, Table S2: The PD-L1CON signature genes for all protein-coding genes and lncRNAs, Figure S1: (A,B) PD-L1 protein and mRNA expression in the PD-L1IND samples (n = 10) in comparison to the PD-L1CON (n = 7) melanoma cell lines. The bars represent Median Fluorescence Intensity (MFI) for PD-L1 staining normalised to isotype control by subtraction followed by log2 transformation: log2(MFIantibody-MFIisotype). The mRNA expression levels are log2 transformed TMM normalized values. (C) Correlation between the PD-L1 protein and mRNA expression. (D) PD-L1 protein levels prior to (control group) and following one day of IFNγ treatment for six out of ten PD-L1IND samples. In four out of the ten cell lines PDL1 had previously shown to be inducible with IFNg [6]. The paired t test p value is shown, Figure S2: (A,B) Volcano plot showing the differentially expressed protein-coding genes (A) and lncRNAs (B). The top 10 differentially expressed gene names are shown and moreover the CD274 (PD-L1) and PDCD1LG2 genes (PD-L2) are shown. The x-axis represents log2 fold change and the y-axis represents the −log10 FDR adjusted p value. (C) The CD274 (log2 fold change = 7.1, FDR adjusted p value = 0.001) and (D) PDCD1LG2 mRNA expression (log2 fold change = 5.0, FDR adjusted p value = 0.03) in the PD-L1CON in comparison to the PD-L1IND melanoma cell lines, Figure S3: (A) Barcode plots showing eight gene sets with differential expression changes in the PD-L1CON melanoma cell lines (red or right side) compared to the PD-L1IND cell lines (blue or left side). The x-axis shows the ranked list of all protein-coding genes according to log2 fold-change. The vertical bars on the x-axis represent the genes that are included in the genesets. The y-axis shows the enrichment score. The FDR adjusted p values using the CAMERA test are shown. (B) Correlation of CD274 expression with the oxidation phosphorylation score in the four external melanoma cell line datasets. The samples (x-axis) are ranked according to CD274 expression (log2 transformed and corresponds to the right side y-axis) and the CD274 expression level is shown by the black line. Grey vertical bars indicate the oxidative phosphorylation scores and corresponds to the left side y-axis. The moving average of the oxidative phosphorylation score is represented by the blue line. The Pearson correlation between the moving average and CD274 expression is shown. (C) The x-axis shows the significantly upregulated (n = 462), downregulated (n = 298) protein-coding genes and the upregulated (n = 106) and downregulated (n = 71) lncRNAs in the PD-L1CON samples compared to the PD-L1IND samples. The y-axis shows the correlation (Pearson) of the differentially expressed genes with the oxidative phosphorylation score. The x-axis is ranked according the lowest to highest correlation. The percentage value of genes higher than 0.25 or lower than −0.25 Pearson correlation value out of all genes in the respective gene set are shown, Figure S4: Hierarchical clustering of the PD-L1CON and PD-L1IND melanoma cell lines using melanoma invasive expression signature genes. The gene signatures were derived from two independent studies: (A) Widmer et al. [45] and (B) Jeff et al. [44], Figure S5: (A,B) Z-scores normalised to baseline (before stimulation) for CD274 expression, and the six scores related to PD-L1CON expression are shown for eight melanoma cell lines following IFN-γ stimulation for 2–5 weeks or TNF stimulation for 3 days (GSE152755). Four samples with a differentiated (pink dots) and four with a dedifferentiated (blue dots) expression signature at baseline are shown, Figure S6: (A) Hierarchical clustering of the PD-L1CON and the PDL1IND samples using receptor tyrosine kinase genes. Receptor tyrosine kinase genes that were differentially expressed in the PD-L1CON cell lines with an FDR adjusted p value < 0.05 are indicated by * whereas between 0.05 and 0.1 are indicated by + (B-I) mRNA expression levels of in the PD-L1CON samples compared to the PD-L1IND samples of various genes that have been found to confer resistance to BRAFi in BRAF mutant melanomas. FDR adjusted p values are shown, Figure S7: (A) Z-scores normalised to baseline (pre-treatment tumour) for CD274 expression and the seven scores related to PD-L1CON expression in the tumours that gained BRAF mutations upon acquired resistance (labelled “BRAF_mutation”, green) in comparison to tumours that did not gain any known mutational resistance (labelled “No_mutation”, purple). (B) Expression levels of 12 TFs in the tumours that gained BRAF mutations upon acquired resistance (labelled “BRAF_mutation”) in comparison to tumours that did not gain any known mutational resistance (labelled “No_mutation”). All values were normalised to the baseline (pre-treatment) tumour from the same patient. The t-test p value is shown.

Author Contributions

Conceptualization, A.A., A.C., E.J.R., P.H. and M.R.E.; methodology, A.A.; software, G.G. and P.A.S.; formal analysis, A.A., J.M.; resources, P.H. and M.R.E.; data curation, A.A., J.M., G.G. and M.P.; writing—original draft preparation, A.A.; writing—review and editing, A.A., E.J.R., J.M., G.G., P.A.S., M.P., P.H., A.C. and M.R.E.; funding acquisition, A.C. and M.R.E. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by grants from the Royal Society of New Zealand Rutherford Discovery Fellowship Program (to AC), the Maurice Wilkins Centre for Molecular Biodiscovery, the Health Research Council of New Zealand grant number 18-144, the Cancer Society of New Zealand grant number 16.06, and Cancer Research Trust NZ.

Institutional Review Board Statement

Ethical review and approval were waived for this study because the study used human cell lines, which had already been established, and human data from online data repositories, and all ethical approvals for use in this context had already been obtained.

Informed Consent Statement

Patient consent was waived due to the use of established human cell lines, and/or use of data from online data repositories.

Data Availability Statement

Transcriptomic data for PD-L1CON and PDL1IND cell lines are available at Database: NCBI GEO, accession number GSE107622 and GSE153595.

Acknowledgments

The authors thank Sarah Diermeier and Nicole Cloonan for comments and suggestions, and Glen Boyle for providing the MM127, MM595 and COLOF239F melanoma cell lines.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Domingues, B.; Lopes, J.M.; Soares, P.; Populo, H. Melanoma treatment in review. Immunotargets Ther. 2018, 7, 35–49. [Google Scholar] [CrossRef] [Green Version]
  2. Sun, C.; Mezzadra, R.; Schumacher, T.N. Regulation and Function of the PD-L1 Checkpoint. Immunity 2018, 48, 434–452. [Google Scholar] [CrossRef] [Green Version]
  3. Robert, C.; Schachter, J.; Long, G.V.; Arance, A.; Grob, J.J.; Mortier, L.; Daud, A.; Carlino, M.S.; McNeil, C.; Lotem, M.; et al. Pembrolizumab versus Ipilimumab in Advanced Melanoma. N. Engl. J. Med. 2015, 372, 2521–2532. [Google Scholar] [CrossRef]
  4. Audrito, V.; Serra, S.; Stingi, A.; Orso, F.; Gaudino, F.; Bologna, C.; Neri, F.; Garaffo, G.; Nassini, R.; Baroni, G.; et al. PD-L1 up-regulation in melanoma increases disease aggressiveness and is mediated through miR-17-5p. Oncotarget 2017, 8, 15894–15911. [Google Scholar] [CrossRef] [Green Version]
  5. Song, C.; Piva, M.; Sun, L.; Hong, A.; Moriceau, G.; Kong, X.; Zhang, H.; Lomeli, S.; Qian, J.; Yu, C.C.; et al. Recurrent Tumor Cell-Intrinsic and -Extrinsic Alterations during MAPKi-Induced Melanoma Regression and Early Adaptation. Cancer Discov. 2017, 7, 1248–1265. [Google Scholar] [CrossRef] [Green Version]
  6. Chatterjee, A.; Rodger, E.J.; Ahn, A.; Stockwell, P.A.; Parry, M.; Motwani, J.; Gallagher, S.J.; Shklovskaya, E.; Tiffen, J.; Eccles, M.R.; et al. Marked Global DNA Hypomethylation Is Associated with Constitutive PD-L1 Expression in Melanoma. iScience 2018, 4, 312–325. [Google Scholar] [CrossRef] [PubMed]
  7. Emran, A.A.; Chatterjee, A.; Rodger, E.J.; Tiffen, J.C.; Gallagher, S.J.; Eccles, M.R.; Hersey, P. Targeting DNA Methylation and EZH2 Activity to Overcome Melanoma Resistance to Immunotherapy. Trends Immunol. 2019, 40, 328–344. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Lee, S.J.; Jang, B.C.; Lee, S.W.; Yang, Y.I.; Suh, S.I.; Park, Y.M.; Oh, S.; Shin, J.G.; Yao, S.; Chen, L.; et al. Interferon regulatory factor-1 is prerequisite to the constitutive expression and IFN-gamma-induced upregulation of B7-H1 (CD274). FEBS Lett. 2006, 580, 755–762. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Garcia-Diaz, A.; Shin, D.S.; Moreno, B.H.; Saco, J.; Escuin-Ordinas, H.; Rodriguez, G.A.; Zaretsky, J.M.; Sun, L.; Hugo, W.; Wang, X.; et al. Interferon Receptor Signaling Pathways Regulating PD-L1 and PD-L2 Expression. Cell Rep. 2017, 19, 1189–1201. [Google Scholar] [CrossRef] [Green Version]
  10. Spranger, S.; Spaapen, R.M.; Zha, Y.; Williams, J.; Meng, Y.; Ha, T.T.; Gajewski, T.F. Up-regulation of PD-L1, IDO, and T(regs) in the melanoma tumor microenvironment is driven by CD8(+) T cells. Sci. Transl. Med. 2013, 5, 200ra116. [Google Scholar] [CrossRef] [Green Version]
  11. Donia, M.; Andersen, R.; Kjeldsen, J.W.; Fagone, P.; Munir, S.; Nicoletti, F.; Andersen, M.H.; Thor Straten, P.; Svane, I.M. Aberrant Expression of MHC Class II in Melanoma Attracts Inflammatory Tumor-Specific CD4+ T- Cells, Which Dampen CD8+ T-cell Antitumor Reactivity. Cancer Res. 2015, 75, 3747–3759. [Google Scholar] [CrossRef] [Green Version]
  12. Shen, X.; Zhang, L.; Li, J.; Li, Y.; Wang, Y.; Xu, Z.X. Recent Findings in the Regulation of Programmed Death Ligand 1 Expression. Front. Immunol. 2019, 10, 1337. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Chatterjee, A.; Stockwell, P.A.; Ahn, A.; Rodger, E.J.; Leichter, A.L.; Eccles, M.R. Genome-wide methylation sequencing of paired primary and metastatic cell lines identifies common DNA methylation changes and a role for EBF3 as a candidate epigenetic driver of melanoma metastasis. Oncotarget 2017, 8, 6085–6101. [Google Scholar] [CrossRef] [Green Version]
  14. Johansson, P.; Pavey, S.; Hayward, N. Confirmation of a BRAF mutation-associated gene expression signature in melanoma. Pigment. Cell Res. 2007, 20, 216–221. [Google Scholar] [CrossRef] [PubMed]
  15. Hoek, K.S.; Schlegel, N.C.; Brafford, P.; Sucker, A.; Ugurel, S.; Kumar, R.; Weber, B.L.; Nathanson, K.L.; Phillips, D.J.; Herlyn, M.; et al. Metastatic potential of melanomas defined by specific gene expression profiles with no BRAF signature. Pigment. Cell Res. 2006, 19, 290–302. [Google Scholar] [CrossRef] [PubMed]
  16. Muller, J.; Krijgsman, O.; Tsoi, J.; Robert, L.; Hugo, W.; Song, C.; Kong, X.; Possik, P.A.; Cornelissen-Steijger, P.D.; Geukes Foppen, M.H.; et al. Low MITF/AXL ratio predicts early resistance to multiple targeted drugs in melanoma. Nat. Commun. 2014, 5, 5712. [Google Scholar] [CrossRef]
  17. Tsoi, J.; Robert, L.; Paraiso, K.; Galvan, C.; Sheu, K.M.; Lay, J.; Wong, D.J.L.; Atefi, M.; Shirazi, R.; Wang, X.; et al. Multi-stage Differentiation Defines Melanoma Subtypes with Differential Vulnerability to Drug-Induced Iron-Dependent Oxidative Stress. Cancer Cell. 2018, 33, 890–904.e5. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [Green Version]
  19. Hugo, W.; Shi, H.; Sun, L.; Piva, M.; Song, C.; Kong, X.; Moriceau, G.; Hong, A.; Dahlman, K.B.; Johnson, D.B.; et al. Non-genomic and Immune Evolution of Melanoma Acquiring MAPKi Resistance. Cell 2015, 162, 1271–1285. [Google Scholar] [CrossRef] [Green Version]
  20. Kim, Y.J.; Sheu, K.M.; Tsoi, J.; Abril-Rodriguez, G.; Medina, E.; Grasso, C.S.; Torrejon, D.Y.; Champhekar, A.S.; Litchfield, K.; Swanton, C.; et al. Melanoma dedifferentiation induced by IFN-gamma epigenetic remodeling in response to anti-PD-1 therapy. J. Clin. Invest. 2021, 131. [Google Scholar] [CrossRef]
  21. Colaprico, A.; Silva, T.C.; Olsen, C.; Garofano, L.; Cava, C.; Garolini, D.; Sabedot, T.S.; Malta, T.M.; Pagnotta, S.M.; Castiglioni, I.; et al. TCGAbiolinks: An R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res. 2016, 44, e71. [Google Scholar] [CrossRef] [PubMed]
  22. Stockwell, P.A.; Chatterjee, A.; Rodger, E.J.; Morison, I.M. DMAP: Differential methylation analysis package for RRBS and WGBS data. Bioinformatics 2014, 30, 1814–1822. [Google Scholar] [CrossRef]
  23. Chatterjee, A.; Stockwell, P.A.; Rodger, E.J.; Morison, I.M. Comparison of alignment software for genome-wide bisulphite sequence data. Nucleic Acids Res. 2012, 40, e79. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Chatterjee, A.; Ahn, A.; Rodger, E.J.; Stockwell, P.A.; Eccles, M.R. A Guide for Designing and Analyzing RNA-Seq Data. Methods Mol. Biol. 2018, 1783, 35–80. [Google Scholar] [CrossRef] [PubMed]
  25. Zheng, H.; Brennan, K.; Hernaez, M.; Gevaert, O. Benchmark of long non-coding RNA quantification for RNA sequencing of cancer samples. Gigascience 2019, 8. [Google Scholar] [CrossRef]
  26. Bray, N.L.; Pimentel, H.; Melsted, P.; Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 2016, 34, 525–527. [Google Scholar] [CrossRef] [PubMed]
  27. Soneson, C.; Love, M.I.; Robinson, M.D. Differential analyses for RNA-seq: Transcript-level estimates improve gene-level inferences. F1000Research 2015, 4, 1521. [Google Scholar] [CrossRef]
  28. Lun, A.T.; Chen, Y.; Smyth, G.K. It’s DE-licious: A Recipe for Differential Expression Analyses of RNA-seq Experiments Using Quasi-Likelihood Methods in edgeR. Methods Mol. Biol. 2016, 1418, 391–416. [Google Scholar] [CrossRef]
  29. Liberzon, A.; Birger, C.; Thorvaldsdottir, H.; Ghandi, M.; Mesirov, J.P.; Tamayo, P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015, 1, 417–425. [Google Scholar] [CrossRef] [Green Version]
  30. Wu, D.; Smyth, G.K. Camera: A competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res. 2012, 40, e133. [Google Scholar] [CrossRef]
  31. Barbie, D.A.; Tamayo, P.; Boehm, J.S.; Kim, S.Y.; Moody, S.E.; Dunn, I.F.; Schinzel, A.C.; Sandy, P.; Meylan, E.; Scholl, C.; et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 2009, 462, 108–112. [Google Scholar] [CrossRef]
  32. Hanzelmann, S.; Castelo, R.; Guinney, J. GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinform. 2013, 14, 7. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Chiappinelli, K.B.; Strissel, P.L.; Desrichard, A.; Li, H.; Henke, C.; Akman, B.; Hein, A.; Rote, N.S.; Cope, L.M.; Snyder, A.; et al. Inhibiting DNA Methylation Causes an Interferon Response in Cancer via dsRNA Including Endogenous Retroviruses. Cell 2015, 162, 974–986. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Roulois, D.; Loo Yau, H.; Singhania, R.; Wang, Y.; Danesh, A.; Shen, S.Y.; Han, H.; Liang, G.; Jones, P.A.; Pugh, T.J.; et al. DNA-Demethylating Agents Target Colorectal Cancer Cells by Inducing Viral Mimicry by Endogenous Transcripts. Cell 2015, 162, 961–973. [Google Scholar] [CrossRef] [Green Version]
  35. Rooney, M.S.; Shukla, S.A.; Wu, C.J.; Getz, G.; Hacohen, N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 2015, 160, 48–61. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Becht, E.; Giraldo, N.A.; Lacroix, L.; Buttard, B.; Elarouci, N.; Petitprez, F.; Selves, J.; Laurent-Puig, P.; Sautes-Fridman, C.; Fridman, W.H.; et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 2016, 17, 218. [Google Scholar] [CrossRef]
  37. Riesenberg, S.; Groetchen, A.; Siddaway, R.; Bald, T.; Reinhardt, J.; Smorra, D.; Kohlmeyer, J.; Renn, M.; Phung, B.; Aymans, P.; et al. MITF and c-Jun antagonism interconnects melanoma dedifferentiation with pro-inflammatory cytokine responsiveness and myeloid cell recruitment. Nat. Commun. 2015, 6, 8755. [Google Scholar] [CrossRef]
  38. Gowrishankar, K.; Gunatilake, D.; Gallagher, S.J.; Tiffen, J.; Rizos, H.; Hersey, P. Inducible but not constitutive expression of PD-L1 in human melanoma cells is dependent on activation of NF-kappaB. PLoS ONE 2015, 10, e0123410. [Google Scholar] [CrossRef] [Green Version]
  39. Gao, J.; Shi, L.Z.; Zhao, H.; Chen, J.; Xiong, L.; He, Q.; Chen, T.; Roszik, J.; Bernatchez, C.; Woodman, S.E.; et al. Loss of IFN-gamma Pathway Genes in Tumor Cells as a Mechanism of Resistance to Anti-CTLA-4 Therapy. Cell 2016, 167, 397–404.e9. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. Zaretsky, J.M.; Garcia-Diaz, A.; Shin, D.S.; Escuin-Ordinas, H.; Hugo, W.; Hu-Lieskovan, S.; Torrejon, D.Y.; Abril-Rodriguez, G.; Sandoval, S.; Barthly, L.; et al. Mutations Associated with Acquired Resistance to PD-1 Blockade in Melanoma. N. Engl. J. Med. 2016, 375, 819–829. [Google Scholar] [CrossRef]
  41. Alavi, S.; Stewart, A.J.; Kefford, R.F.; Lim, S.Y.; Shklovskaya, E.; Rizos, H. Interferon Signaling Is Frequently Downregulated in Melanoma. Front. Immunol. 2018, 9, 1414. [Google Scholar] [CrossRef]
  42. Sucker, A.; Zhao, F.; Pieper, N.; Heeke, C.; Maltaner, R.; Stadtler, N.; Real, B.; Bielefeld, N.; Howe, S.; Weide, B.; et al. Acquired IFNgamma resistance impairs anti-tumor immunity and gives rise to T-cell-resistant melanoma lesions. Nat Commun. 2017, 8, 15440. [Google Scholar] [CrossRef]
  43. Shin, D.S.; Zaretsky, J.M.; Escuin-Ordinas, H.; Garcia-Diaz, A.; Hu-Lieskovan, S.; Kalbasi, A.; Grasso, C.S.; Hugo, W.; Sandoval, S.; Torrejon, D.Y.; et al. Primary Resistance to PD-1 Blockade Mediated by JAK1/2 Mutations. Cancer Discov. 2017, 7, 188–201. [Google Scholar] [CrossRef] [Green Version]
  44. Jeffs, A.R.; Glover, A.C.; Slobbe, L.J.; Wang, L.; He, S.; Hazlett, J.A.; Awasthi, A.; Woolley, A.G.; Marshall, E.S.; Joseph, W.R.; et al. A gene expression signature of invasive potential in metastatic melanoma cells. PLoS ONE 2009, 4, e8461. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Widmer, D.S.; Cheng, P.F.; Eichhoff, O.M.; Belloni, B.C.; Zipser, M.C.; Schlegel, N.C.; Javelaud, D.; Mauviel, A.; Dummer, R.; Hoek, K.S. Systematic classification of melanoma cells by phenotype-specific gene expression mapping. Pigment. Cell Melanoma Res. 2012, 25, 343–353. [Google Scholar] [CrossRef]
  46. Landsberg, J.; Kohlmeyer, J.; Renn, M.; Bald, T.; Rogava, M.; Cron, M.; Fatho, M.; Lennerz, V.; Wolfel, T.; Holzel, M.; et al. Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation. Nature 2012, 490, 412–416. [Google Scholar] [CrossRef] [PubMed]
  47. Mehta, A.; Kim, Y.J.; Robert, L.; Tsoi, J.; Comin-Anduix, B.; Berent-Maoz, B.; Cochran, A.J.; Economou, J.S.; Tumeh, P.C.; Puig-Saus, C.; et al. Immunotherapy Resistance by Inflammation-Induced Dedifferentiation. Cancer Discov. 2018, 8, 935–943. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Lambert, S.A.; Jolma, A.; Campitelli, L.F.; Das, P.K.; Yin, Y.; Albu, M.; Chen, X.; Taipale, J.; Hughes, T.R.; Weirauch, M.T. The Human Transcription Factors. Cell 2018, 172, 650–665. [Google Scholar] [CrossRef] [Green Version]
  49. Wang, Y.; Wang, H.; Yao, H.; Li, C.; Fang, J.Y.; Xu, J. Regulation of PD-L1: Emerging Routes for Targeting Tumor Immune Evasion. Front. Pharmacol. 2018, 9, 536. [Google Scholar] [CrossRef]
  50. Zerdes, I.; Matikas, A.; Bergh, J.; Rassidakis, G.Z.; Foukakis, T. Genetic, transcriptional and post-translational regulation of the programmed death protein ligand 1 in cancer: Biology and clinical correlations. Oncogene 2018, 37, 4639–4661. [Google Scholar] [CrossRef] [Green Version]
  51. Hugo, W.; Zaretsky, J.M.; Sun, L.; Song, C.; Moreno, B.H.; Hu-Lieskovan, S.; Berent-Maoz, B.; Pang, J.; Chmielowski, B.; Cherry, G.; et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell 2016, 165, 35–44. [Google Scholar] [CrossRef] [Green Version]
  52. Shaffer, S.M.; Dunagin, M.C.; Torborg, S.R.; Torre, E.A.; Emert, B.; Krepler, C.; Beqiri, M.; Sproesser, K.; Brafford, P.A.; Xiao, M.; et al. Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature 2017, 546, 431–435. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Natarajan, V.T.; Ganju, P.; Singh, A.; Vijayan, V.; Kirty, K.; Yadav, S.; Puntambekar, S.; Bajaj, S.; Dani, P.P.; Kar, H.K.; et al. IFN-gamma signaling maintains skin pigmentation homeostasis through regulation of melanosome maturation. Proc. Natl. Acad. Sci. USA 2014, 111, 2301–2306. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Sun, C.; Wang, L.; Huang, S.; Heynen, G.J.; Prahallad, A.; Robert, C.; Haanen, J.; Blank, C.; Wesseling, J.; Willems, S.M.; et al. Reversible and adaptive resistance to BRAF(V600E) inhibition in melanoma. Nature 2014, 508, 118–122. [Google Scholar] [CrossRef] [PubMed]
  55. Konieczkowski, D.J.; Johannessen, C.M.; Abudayyeh, O.; Kim, J.W.; Cooper, Z.A.; Piris, A.; Frederick, D.T.; Barzily-Rokni, M.; Straussman, R.; Haq, R.; et al. A melanoma cell state distinction influences sensitivity to MAPK pathway inhibitors. Cancer Discov. 2014, 4, 816–827. [Google Scholar] [CrossRef] [Green Version]
  56. Verfaillie, A.; Imrichova, H.; Atak, Z.K.; Dewaele, M.; Rambow, F.; Hulselmans, G.; Christiaens, V.; Svetlichnyy, D.; Luciani, F.; Van den Mooter, L.; et al. Decoding the regulatory landscape of melanoma reveals TEADS as regulators of the invasive cell state. Nat. Commun. 2015, 6, 6683. [Google Scholar] [CrossRef] [Green Version]
  57. Ji, Z.; Erin Chen, Y.; Kumar, R.; Taylor, M.; Jenny Njauw, C.N.; Miao, B.; Frederick, D.T.; Wargo, J.A.; Flaherty, K.T.; Jonsson, G.; et al. MITF Modulates Therapeutic Resistance through EGFR Signaling. J. Invest. Dermatol. 2015, 135, 1863–1872. [Google Scholar] [CrossRef] [Green Version]
  58. Martinez-Zamudio, R.I.; Roux, P.F.; de Freitas, J.; Robinson, L.; Dore, G.; Sun, B.; Belenki, D.; Milanovic, M.; Herbig, U.; Schmitt, C.A.; et al. AP-1 imprints a reversible transcriptional programme of senescent cells. Nat. Cell Biol. 2020, 22, 842–855. [Google Scholar] [CrossRef]
  59. Bi, M.; Zhang, Z.; Jiang, Y.Z.; Xue, P.; Wang, H.; Lai, Z.; Fu, X.; De Angelis, C.; Gong, Y.; Gao, Z.; et al. Enhancer reprogramming driven by high-order assemblies of transcription factors promotes phenotypic plasticity and breast cancer endocrine resistance. Nat. Cell Biol. 2020, 22, 701–715. [Google Scholar] [CrossRef]
  60. Watt, A.C.; Cejas, P.; DeCristo, M.J.; Metzger-Filho, O.; Lam, E.Y.N.; Qiu, X.; BrinJones, H.; Kesten, N.; Coulson, R.; Font-Tello, A.; et al. CDK4/6 inhibition reprograms the breast cancer enhancer landscape by stimulating AP-1 transcriptional activity. Nat. Cancer 2020. [Google Scholar] [CrossRef]
  61. Choukrallah, M.A.; Matthias, P. The Interplay between Chromatin and Transcription Factor Networks during B Cell Development: Who Pulls the Trigger First? Front. Immunol. 2014, 5, 156. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Sharma, A.; Cao, E.Y.; Kumar, V.; Zhang, X.; Leong, H.S.; Wong, A.M.L.; Ramakrishnan, N.; Hakimullah, M.; Teo, H.M.V.; Chong, F.T.; et al. Longitudinal single-cell RNA sequencing of patient-derived primary cells reveals drug-induced infidelity in stem cell hierarchy. Nat. Commun. 2018, 9, 4931. [Google Scholar] [CrossRef] [Green Version]
  63. Adam, R.C.; Yang, H.; Rockowitz, S.; Larsen, S.B.; Nikolova, M.; Oristian, D.S.; Polak, L.; Kadaja, M.; Asare, A.; Zheng, D.; et al. Pioneer factors govern super-enhancer dynamics in stem cell plasticity and lineage choice. Nature 2015, 521, 366–370. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Ge, Y.; Gomez, N.C.; Adam, R.C.; Nikolova, M.; Yang, H.; Verma, A.; Lu, C.P.; Polak, L.; Yuan, S.; Elemento, O.; et al. Stem Cell Lineage Infidelity Drives Wound Repair and Cancer. Cell 2017, 169, 636–650.e14. [Google Scholar] [CrossRef] [PubMed]
Cancers 13 04250 g001 550
Figure 1. The PD-L1CON group harbours distinct global mRNA and long non-coding RNA expression profiles compared to the PD-L1IND group. (A) A schematic diagram showing the four subgroups of melanoma tumours based on the levels of tumour infiltrative lymphocytes and PD-L1 protein expression [6]. We used melanoma cell lines and compared transcriptomes between cell lines corresponding to the low PD-L1 (group 2: TIL-/PDL1- or PD-L1IND) and high PD-L1 expression groups (group 4: TIL-/PDL1+ or PD-L1CON). Representative flow cytometry figure of the PD-L1 expression levels in each of the cell lines, representing the two groups are shown below. The flow cytometry figure shows isotype control (red), PD-L1 stained (blue), and PD-L1 stained following 24-hour IFNγ induction (green). The x-axis represents PD-L1 expression levels (log2) and y-axis represents cell counts. (B,C) Principal component analysis (PCA). Shows PC1 (x-axis) and PC2 (y-axis) for the PD-L1IND samples and the PD-L1CON samples using the top 500 protein-coding genes (B) or lncRNAs (C) with the highest variance. (D,E) Unsupervised hierarchical clustering of PD-L1CON and PD-L1IND melanoma cell lines using either expression of all protein-coding genes, or of lncRNAs. (F) The significantly differentially expressed mRNA–lncRNA gene pairs that are located within 100 kilobases to one another are shown. The x-axis and the y-axis show the protein-coding gene and lncRNA expression fold change (log2), respectively. Only the protein-coding gene names are shown. The distance between the mRNA–lncRNA gene pair is represented by the size of the circle, and the lncRNA subtype is shown by the colours.
Figure 1. The PD-L1CON group harbours distinct global mRNA and long non-coding RNA expression profiles compared to the PD-L1IND group. (A) A schematic diagram showing the four subgroups of melanoma tumours based on the levels of tumour infiltrative lymphocytes and PD-L1 protein expression [6]. We used melanoma cell lines and compared transcriptomes between cell lines corresponding to the low PD-L1 (group 2: TIL-/PDL1- or PD-L1IND) and high PD-L1 expression groups (group 4: TIL-/PDL1+ or PD-L1CON). Representative flow cytometry figure of the PD-L1 expression levels in each of the cell lines, representing the two groups are shown below. The flow cytometry figure shows isotype control (red), PD-L1 stained (blue), and PD-L1 stained following 24-hour IFNγ induction (green). The x-axis represents PD-L1 expression levels (log2) and y-axis represents cell counts. (B,C) Principal component analysis (PCA). Shows PC1 (x-axis) and PC2 (y-axis) for the PD-L1IND samples and the PD-L1CON samples using the top 500 protein-coding genes (B) or lncRNAs (C) with the highest variance. (D,E) Unsupervised hierarchical clustering of PD-L1CON and PD-L1IND melanoma cell lines using either expression of all protein-coding genes, or of lncRNAs. (F) The significantly differentially expressed mRNA–lncRNA gene pairs that are located within 100 kilobases to one another are shown. The x-axis and the y-axis show the protein-coding gene and lncRNA expression fold change (log2), respectively. Only the protein-coding gene names are shown. The distance between the mRNA–lncRNA gene pair is represented by the size of the circle, and the lncRNA subtype is shown by the colours.
Cancers 13 04250 g001
Cancers 13 04250 g002 550
Figure 2. PD-L1CON samples have a reprogrammed transcriptome inclusive of dedifferentiation, reduced oxidative phosphorylation, and an enhanced IFN signalling pathway. (A) The top 20 gene sets that were found to be differentially enriched using the CAMERA test are shown [30]. The 5 × 10−5 FDR adjusted p value significance threshold is shown by the dotted red line. Gene sets are categorised into those involved in the innate immune response (or cytokine signalling), melanocyte differentiation, oxidative phosphorylation, and others. (B) Heatmap showing the unsupervised clustering of the PD-L1CON and PD-L1IND melanoma cell lines according to the differentiation signature genes from Tsoi and colleagues [17]. Mutation status according to BRAF (V600), NRAS (Q61 and G13), and double wild-type is shown. (C) Correlation of CD274 expression with the IFN score, TNF score and differentiation score in the four external melanoma cell line datasets. The samples (x-axis) are ranked according to CD274 expression (log2 transformed and corresponds to the right-side y-axis) and the CD274 expression level is shown by the black line. Grey vertical bars indicate the respective scores and corresponds to the left side y-axis. The moving average is represented by the blue line. The Pearson correlation between the moving average and CD274 expression is shown. (D, first and second figure) The y-axis represents the correlation value (Pearson) of the CYT score with either the upregulated (D, first figure) or downregulated protein-coding genes (D, second figure) in the PD-L1CON samples. The x-axis represents all the up (462 genes) and downregulated (298 genes) protein-coding genes. The x-axis is ranked according to lowest to highest correlation value (Pearson). (D, third and fourth figure) The y-axis represents the correlation value (Pearson) of the up (D, third figure) and down PD-L1CON scores (D, fourth figure) (calculated with ssGSEA using the differentially expressed protein-coding genes) with the estimated absolute abundance of various immune subtypes (as measured using the MCPcounter bioinformatic tool).
Figure 2. PD-L1CON samples have a reprogrammed transcriptome inclusive of dedifferentiation, reduced oxidative phosphorylation, and an enhanced IFN signalling pathway. (A) The top 20 gene sets that were found to be differentially enriched using the CAMERA test are shown [30]. The 5 × 10−5 FDR adjusted p value significance threshold is shown by the dotted red line. Gene sets are categorised into those involved in the innate immune response (or cytokine signalling), melanocyte differentiation, oxidative phosphorylation, and others. (B) Heatmap showing the unsupervised clustering of the PD-L1CON and PD-L1IND melanoma cell lines according to the differentiation signature genes from Tsoi and colleagues [17]. Mutation status according to BRAF (V600), NRAS (Q61 and G13), and double wild-type is shown. (C) Correlation of CD274 expression with the IFN score, TNF score and differentiation score in the four external melanoma cell line datasets. The samples (x-axis) are ranked according to CD274 expression (log2 transformed and corresponds to the right-side y-axis) and the CD274 expression level is shown by the black line. Grey vertical bars indicate the respective scores and corresponds to the left side y-axis. The moving average is represented by the blue line. The Pearson correlation between the moving average and CD274 expression is shown. (D, first and second figure) The y-axis represents the correlation value (Pearson) of the CYT score with either the upregulated (D, first figure) or downregulated protein-coding genes (D, second figure) in the PD-L1CON samples. The x-axis represents all the up (462 genes) and downregulated (298 genes) protein-coding genes. The x-axis is ranked according to lowest to highest correlation value (Pearson). (D, third and fourth figure) The y-axis represents the correlation value (Pearson) of the up (D, third figure) and down PD-L1CON scores (D, fourth figure) (calculated with ssGSEA using the differentially expressed protein-coding genes) with the estimated absolute abundance of various immune subtypes (as measured using the MCPcounter bioinformatic tool).
Cancers 13 04250 g002
Cancers 13 04250 g003 550
Figure 3. Protein-coding and lncRNA expression signatures of PD-L1CON samples are associated with an increased IFN score and reduced differentiation and oxidative phosphorylation scores in the TCGA SKCM database. (A) The x-axis shows the upregulated (n = 462) (A, top row) and downregulated (n = 298) (A, bottom row) protein-coding genes (A) in the PD-L1CON samples. The y-axis shows the correlation (Pearson) of the protein-coding genes with CD274 expression, IFN score, differentiation score and oxidative phosphorylation score. The x-axis is ranked according the lowest to highest correlation. The percentage of genes higher than 0.25 or lower than −0.25 Pearson correlation value is shown. (B) The x-axis shows the upregulated (n = 106) and downregulated (n = 71) lncRNAs in the PD-L1CON samples. The y-axis shows the correlation (Pearson) of the lncRNAs with CD274 expression, IFN score, differentiation score and oxidative phosphorylation score. The x-axis is ranked according the lowest to highest correlation. The percentage of genes higher than 0.25 or lower than −0.25 Pearson correlation value is shown.
Figure 3. Protein-coding and lncRNA expression signatures of PD-L1CON samples are associated with an increased IFN score and reduced differentiation and oxidative phosphorylation scores in the TCGA SKCM database. (A) The x-axis shows the upregulated (n = 462) (A, top row) and downregulated (n = 298) (A, bottom row) protein-coding genes (A) in the PD-L1CON samples. The y-axis shows the correlation (Pearson) of the protein-coding genes with CD274 expression, IFN score, differentiation score and oxidative phosphorylation score. The x-axis is ranked according the lowest to highest correlation. The percentage of genes higher than 0.25 or lower than −0.25 Pearson correlation value is shown. (B) The x-axis shows the upregulated (n = 106) and downregulated (n = 71) lncRNAs in the PD-L1CON samples. The y-axis shows the correlation (Pearson) of the lncRNAs with CD274 expression, IFN score, differentiation score and oxidative phosphorylation score. The x-axis is ranked according the lowest to highest correlation. The percentage of genes higher than 0.25 or lower than −0.25 Pearson correlation value is shown.
Cancers 13 04250 g003
Cancers 13 04250 g004 550
Figure 4. Correlation of lineage-specific and immunity-associated transcription factor mRNAs with PD-L1 expression in melanoma and association with long non-coding RNAs. (A) Volcano plot showing all the differentially expressed TFs in the PD-L1CON group. The 25 and 19 significantly upregulated and downregulated TF mRNAs in the PD-L1CON samples are shown, respectively. The x-axis represents log2 fold change and the y-axis represents −log10 FDR adjusted p value. The red dashed line represent the FDR adjusted 0.05 p value (y-axis) and the 1 log2 fold change (x-axis). (B) Heatmap showing the correlation value (Pearsons) of the 25 up and 19 down regulated TFs with CD274 expression. Rows represent the TFs and are ordered according to the lowest to highest mean correlation with CD274 expression across the four melanoma cell line datasets and the TCGA SKCM RNA-seq dataset. (C) Correlation was calculated between CD274 expression and all TF (n = 1107) across the five gene expression datasets and the TF were ordered from highest to lowest mean average correlation across the five datasets. The direction of the 25 upregulated and 19 downregulated TFs are shown by the direction sidebar. (D) The significantly differentially expressed TFs to lncRNA pairs that are located within 100 kilobases to each other are shown. The x-axis and the y-axis show the TF mRNA and lncRNA expression fold-change (log2), respectively. The distance between the TF mRNA and lncRNA gene are shown by the size of the circle.
Figure 4. Correlation of lineage-specific and immunity-associated transcription factor mRNAs with PD-L1 expression in melanoma and association with long non-coding RNAs. (A) Volcano plot showing all the differentially expressed TFs in the PD-L1CON group. The 25 and 19 significantly upregulated and downregulated TF mRNAs in the PD-L1CON samples are shown, respectively. The x-axis represents log2 fold change and the y-axis represents −log10 FDR adjusted p value. The red dashed line represent the FDR adjusted 0.05 p value (y-axis) and the 1 log2 fold change (x-axis). (B) Heatmap showing the correlation value (Pearsons) of the 25 up and 19 down regulated TFs with CD274 expression. Rows represent the TFs and are ordered according to the lowest to highest mean correlation with CD274 expression across the four melanoma cell line datasets and the TCGA SKCM RNA-seq dataset. (C) Correlation was calculated between CD274 expression and all TF (n = 1107) across the five gene expression datasets and the TF were ordered from highest to lowest mean average correlation across the five datasets. The direction of the 25 upregulated and 19 downregulated TFs are shown by the direction sidebar. (D) The significantly differentially expressed TFs to lncRNA pairs that are located within 100 kilobases to each other are shown. The x-axis and the y-axis show the TF mRNA and lncRNA expression fold-change (log2), respectively. The distance between the TF mRNA and lncRNA gene are shown by the size of the circle.
Cancers 13 04250 g004
Cancers 13 04250 g005 550
Figure 5. PD-L1CON expression is associated with treatment resistance to MAPK inhibitors when associated with transcriptome reprogramming. (A,B) Five melanoma cell lines that acquired BRAF mutations (splicing and copy number amplification) upon development of resistance and six melanoma cell lines that did not acquire mutations but exhibited a transcriptome reprogrammed state upon development of resistance are shown. CD274 expression and seven scores which includes the IFN score, TNF score, differentiation score, oxidative phosphorylation score, the MPAS score, and the up and down PD-L1CON scores are shown for the parental (before treatment) and after BRAFi resistance (SDR = single drug resistance, DDR = double drug resistance). The bar plots show values normalised to the corresponding baseline sample. Bar plot values are averages across all samples within either the BRAF mutation resistance group or the transcriptome reprogrammed resistance group. FDR adjusted p values for paired t test are shown for each bar. Error bars represent standard error. (C) Volcano plot shows the TFs that were differentially expressed in the PD-L1CON samples and in addition, involved in IFN signalling, TNF signalling, both IFN and TNF signalling (labelled as “IFN_TNF”), and the melanocyte differentiation pathway. The CD274 gene is shown in red. (D,E) Bar plots shows the expression changes of TFs (involved in the IFN, TNF signalling, and melanocyte differentiation) following acquired resistance with either BRAF mutations or transcriptome reprogramming in melanoma cell lines. Error bars represents standard error. (F) Cumulative barplot shows the CD274 expression and the seven scores at different stages of BRAFi resistance for the M229 and M238 melanoma cell lines. Values are normalised to the parental cell line, and resistance stages include (1) two days of BRAFi treatment (BRAFi2D), (2) Drug Tolerant Persisters (DTP) where a small subpopulation of persisting cells remain, (3) Drug Tolerant Proliferative Persisters (DTPP) where proliferation was regained, and (4) single drug resistance (SDR), a permanent resistant state to BRAFi. (G,H) Cumulative barplot shows the mRNA expression changes of TFs involved in the IFN signalling, TNF signalling, and melanocyte differentiation at different stages of BRAFi resistance for the M229 and M238 melanoma cell lines.
Figure 5. PD-L1CON expression is associated with treatment resistance to MAPK inhibitors when associated with transcriptome reprogramming. (A,B) Five melanoma cell lines that acquired BRAF mutations (splicing and copy number amplification) upon development of resistance and six melanoma cell lines that did not acquire mutations but exhibited a transcriptome reprogrammed state upon development of resistance are shown. CD274 expression and seven scores which includes the IFN score, TNF score, differentiation score, oxidative phosphorylation score, the MPAS score, and the up and down PD-L1CON scores are shown for the parental (before treatment) and after BRAFi resistance (SDR = single drug resistance, DDR = double drug resistance). The bar plots show values normalised to the corresponding baseline sample. Bar plot values are averages across all samples within either the BRAF mutation resistance group or the transcriptome reprogrammed resistance group. FDR adjusted p values for paired t test are shown for each bar. Error bars represent standard error. (C) Volcano plot shows the TFs that were differentially expressed in the PD-L1CON samples and in addition, involved in IFN signalling, TNF signalling, both IFN and TNF signalling (labelled as “IFN_TNF”), and the melanocyte differentiation pathway. The CD274 gene is shown in red. (D,E) Bar plots shows the expression changes of TFs (involved in the IFN, TNF signalling, and melanocyte differentiation) following acquired resistance with either BRAF mutations or transcriptome reprogramming in melanoma cell lines. Error bars represents standard error. (F) Cumulative barplot shows the CD274 expression and the seven scores at different stages of BRAFi resistance for the M229 and M238 melanoma cell lines. Values are normalised to the parental cell line, and resistance stages include (1) two days of BRAFi treatment (BRAFi2D), (2) Drug Tolerant Persisters (DTP) where a small subpopulation of persisting cells remain, (3) Drug Tolerant Proliferative Persisters (DTPP) where proliferation was regained, and (4) single drug resistance (SDR), a permanent resistant state to BRAFi. (G,H) Cumulative barplot shows the mRNA expression changes of TFs involved in the IFN signalling, TNF signalling, and melanocyte differentiation at different stages of BRAFi resistance for the M229 and M238 melanoma cell lines.
Cancers 13 04250 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ahn, A.; Rodger, E.J.; Motwani, J.; Gimenez, G.; Stockwell, P.A.; Parry, M.; Hersey, P.; Chatterjee, A.; Eccles, M.R. Transcriptional Reprogramming and Constitutive PD-L1 Expression in Melanoma Are Associated with Dedifferentiation and Activation of Interferon and Tumour Necrosis Factor Signalling Pathways. Cancers 2021, 13, 4250. https://0-doi-org.brum.beds.ac.uk/10.3390/cancers13174250

AMA Style

Ahn A, Rodger EJ, Motwani J, Gimenez G, Stockwell PA, Parry M, Hersey P, Chatterjee A, Eccles MR. Transcriptional Reprogramming and Constitutive PD-L1 Expression in Melanoma Are Associated with Dedifferentiation and Activation of Interferon and Tumour Necrosis Factor Signalling Pathways. Cancers. 2021; 13(17):4250. https://0-doi-org.brum.beds.ac.uk/10.3390/cancers13174250

Chicago/Turabian Style

Ahn, Antonio, Euan J. Rodger, Jyoti Motwani, Gregory Gimenez, Peter A. Stockwell, Matthew Parry, Peter Hersey, Aniruddha Chatterjee, and Michael R. Eccles. 2021. "Transcriptional Reprogramming and Constitutive PD-L1 Expression in Melanoma Are Associated with Dedifferentiation and Activation of Interferon and Tumour Necrosis Factor Signalling Pathways" Cancers 13, no. 17: 4250. https://0-doi-org.brum.beds.ac.uk/10.3390/cancers13174250

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop