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

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.


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-L1 IND . 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-L1 CON [2,12].
In this study, we describe transcriptomic features of melanoma cell lines with constitutive high expression of PD-L1 (PD-L1 CON ) compared to melanoma cell lines with low inducible levels of PD-L1 expression (PD-L1 IND ). We found that PD-L1 CON 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-L1 CON mRNA expression, and with factors associated with resistance to MAPK pathway inhibitors.

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-L1 IND group and >10,000 for PD-L1 CON group. Gene expression values confirmed that there were significantly higher levels of CD274 expression in the PD-L1 CON group, with an increased log 2 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.

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).

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].

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.

Generation and Processing of Transcriptome Data
RNA-seq processing for PD-L1 IND and PD-L1 CON 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].

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 arti-cles [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 crestlike, (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]. 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 cellintrinsic PD-L1 expression (i.e., PD-L1 constitutive expression, or PD-L1 CON ) 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 PDL1 IND ), or at a constitutively high level (PDL1 CON ). 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-L1 CON cell lines and ten PD-L1 IND 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-L1 IND cell lines in our panel were inducible, which confirmed them as being able to respond to IFNγ induction ( Figure S1D), while PD-L1 CON weren't inducible [38]. , 6 of 23 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-L1 CON and PD-L1 IND melanoma subgroups. Differential expression analysis identified 462 upregulated and 298 downregulated protein coding genes in the PD-L1 CON group, compared to the PD-L1 IND group (FDR adjusted p value < 0.05) ( Table S2). As expected, CD274 expression was significantly increased in the PD-L1 CON group (FDR adjusted p value = 0.001, log 2 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-L1 CON group (FDR adjusted p value = 0.03, log 2 FC = 5.0, Figure S2A,D). For the long non-coding RNA (lncRNA) analysis, we identified 106 upregulated and 71 downregulated lncRNA in PD-L1 CON lines compared to the PD-L1 IND 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-L1 CON lines compared to the PD-L1 IND group. These results suggest coordinated expression changes for the mRNA and lncRNA expression. Overall, PD-L1 CON 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.

PD-L1 CON 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-L1 CON cell lines compared to PD-L1 IND cell lines, we performed gene-set enrichment analysis using CAM-ERA [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-L1 CON 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-L1 CON 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-L1 CON group (FDR adjusted p value = 7.5 × 10 −4 and 2.1 × 10 −4 , respectively) ( Figure S3A). These findings suggest that PD-L1 CON 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-L1 CON melanoma cells [6], and were a characteristic feature associated with PD-L1 expression. The downregulated genes in the PD-L1 CON 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-L1 CON cell lines were enriched for genes characteristic of the undifferentiated state ( Figure 2B). The PD-L1 IND 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-L1 CON 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-L1 CON samples were enriched for an invasive gene expression profile. Using two independent invasive gene expression gene sets [44,45], we found PD-L1 CON cell lines had higher levels of expression of genes involved in melanoma invasion ( Figure S4A,B) although three PD-L1 IND cell lines (NZM22, WM115, CM145.pre) had weak expression levels of invasive genes, whereas CM138 clustered more closely to the PD-L1 CON cell lines.

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-L1 CON 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-L1 CON samples (referred to as the up and down PD-L1 CON 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-L1 CON samples did not alter the PD-L1 CON expression in two of the PD-L1 CON samples that were also used in this study [38], which demonstrates that PD-L1 CON 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-L1 CON expression is similar to that in tumours with an immune cell dependent PD-L1 expression, we asked whether the PD-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON samples, these signatures are not dependent on any external stimulation.

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-L1 CON 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-L1 CON cell lines as a surrogate for PD-L1 CON expression. A large proportion of the upregulated proteincoding genes (n = 462) found in the PD-L1 CON 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 lncR-NAs ( Figure 3B). Collectively, these results support the finding that PD-L1 CON expression is associated with an enhanced IFN and TNF signalling, and reduced differentiation and oxidative phosphorylation.  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-L1 CON 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.

Lineage Specific, TNF, and IFN Associated Transcription Factors Are Differentially Expressed in the PD-L1 CON 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-L1 CON 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-L1 CON 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 (log 2 FC = −10.7, FDR adjusted p value = 0.0003) and MITF expression was significantly reduced in expression (log 2 FC= −3.8, FDR adjusted p value = 0.04) in the PD-L1 CON cell lines ( Figure S6E,F). The upregulation of these TF expression further supports our observation that the PD-L1 CON 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).

The PD-L1 CON 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-L1 CON cell lines had increased expression of cMET (log 2 FC = 2.6, FDR adjusted p value = 0.07), reduced expression of LEF1 (log 2 FC = −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-L1 CON 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-L1 CON   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-L1 CON 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-L1 CON scores (Figure 5B), and a corresponding reduction in differentiation, oxidative phosphorylation, MPAS, and "down" PD-L1 CON scores. Moreover, we found that TFs associated with PD-L1 CON 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-L1 CON 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-L1 CON 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).

Transcriptomic Changes, including PD-L1 CON 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-L1 CON "up" score after 2 days of BRAFi treatment ( Figure 5F), suggesting that the PD-L1 CON 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-L1 CON "up" scores were largely increased at the DTPP stage, whereas dedifferentiation and PD-L1 CON "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-L1 CON 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.  . (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 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-L1 CON 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-L1 CON 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

The PD-L1 CON Expression Signature Is Associated with Transcriptomic Reprogramming of Melanomas following MAPK Inhibitor Resistance in Patients
We next investigated whether the PD-L1 CON 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-L1 CON 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-L1 CON scores (normalised to matched pre-treatment tumours from the same patient), while oxidative phosphorylation, down PD-L1 CON , 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-L1 CON expression signature accompanies transcriptomic reprogramming and treatment resistance in melanomas lacking a BRAF mutational resistance mechanism.

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-L1 CON ). Overall, our study has shown that PD-L1 CON 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-L1 CON cells is associated with constitutively active IFN and TNF signalling pathways [17,37,46,53]. Furthermore, given that the PD-L1 CON and PD-L1 IND 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-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON 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-L1 CON expression prior to treatment exposure (CM143-pre), suggesting that, in this patient, a PD-L1 CON 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.

Conclusions
Our study offers new insights into transcriptomic patterns associated with PD-L1 CON 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-L1 CON cells, along with a dedifferentiation gene expression profile and increased expression of RTK genes which highlights the importance of PD-L1 CON 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-L1 CON 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-L1 CON , in cancer drug-resistant phenotypes.  Table S1: Information on the mutation status and the PD-L1 IND and  PD-L1 CON group, Table S2: The PD-L1 CON signature genes for all protein-coding genes and lncRNAs, Figure S1: (A,B) PD-L1 protein and mRNA expression in the PD-L1 IND samples (n = 10) in comparison to the PD-L1 CON (n = 7) melanoma cell lines. The bars represent Median Fluorescence Intensity (MFI) for PD-L1 staining normalised to isotype control by subtraction followed by log 2 transformation: log 2 (MFI antibody -MFI isotype ). The mRNA expression levels are log 2 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-L1 IND 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 . The x-axis shows the ranked list of all protein-coding genes according to log 2 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 (log 2 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-L1 CON samples compared to the PD-L1 IND 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-L1 CON and PD-L1 IND 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-L1 CON 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-L1 CON and the PDL1 IND samples using receptor tyrosine kinase genes. Receptor tyrosine kinase genes that were differentially expressed in the PD-L1 CON 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-L1 CON samples compared to the PD-L1 IND 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-L1 CON 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. 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-L1 CON and PDL1 IND cell lines are available at Database: NCBI GEO, accession number GSE107622 and GSE153595.