RNA-seq Characterization of Melanoma Phenotype Switch in 3D Collagen after p38 MAPK Inhibitor Treatment

Melanoma phenotype plasticity underlies tumour dissemination and resistance to therapy, yet its regulation is incompletely understood. In vivo switching between a more differentiated, proliferative phenotype and a dedifferentiated, invasive phenotype is directed by the tumour microenvironment. We found that treatment of partially dedifferentiated, invasive A375M2 cells with two structurally unrelated p38 MAPK inhibitors, SB2021920 and BIRB796, induces a phenotype switch in 3D collagen, as documented by increased expression of melanocyte differentiation markers and a loss of invasive phenotype markers. The phenotype is accompanied by morphological change corresponding to amoeboid–mesenchymal transition. We performed RNA sequencing with an Illumina HiSeq platform to fully characterise transcriptome changes underlying the switch. Gene expression results obtained with RNA-seq were validated by comparing them with RT-qPCR. Transcriptomic data generated in the study will extend the present understanding of phenotype plasticity in melanoma and its contribution to invasion and metastasis.


Introduction
Despite significant progress in therapy, melanoma remains a life-threatening disease with a high risk of early metastasis. Understanding the biological processes underlying melanoma metastasis is expected to generate new advances and therapeutic modalities in melanoma management [1]. Melanoma is a cancer arising from melanocytes, pigmentproducing cells located in the epidermis and a few other tissues. It shows unique biological features related to its developmental origin, namely specific differentiation plasticity and migratory behaviour. The central or master regulator of the melanocyte lineage is the transcription factor MITF [2,3]. MITF regulates expression of pigmentation-related genes, and coordinates differentiation and proliferation in developing melanocytes. In addition, it affects many other biological processes, including cell survival, invasion, senescence, metabolism, autophagy, lysosomal biogenesis, and DNA damage repair [2]. MITF is regulated at transcription and multiple post-transcriptional levels. Its transcription in the melanocyte lineage is driven by SOX10, WNT/β-catenin and CREB, and is directly or indirectly suppressed by NF-κB, JNK, hedgehog, TGFβ, activated RhoA, HIF1α, Brn2/POU3F2 and ATF4 [2][3][4][5][6].
Cancer cells, including melanoma, can individually invade a 3D environment in either a mesenchymal mode, which depends on integrin-mediated adhesion and extracellular protease activity, or in an amoeboid mode that is independent of integrins and localised proteolysis, and is also characterised by largely dissociated microtubules and high RhoA-ROCK-dependent contractility [36][37][38]. As RhoA-ROCK signalling is also required for mesenchymal invasion, it is an attractive target for antimetastatic therapy [39,40]. A375M2 metastatic melanoma cells have been a popular model for studying amoeboid invasion and its regulation, for they show prominent blebbing amoeboid morphology in, or on top of, soft collagen matrices, good long-term viability under such conditions, and, Cancer cells, including melanoma, can individually invade a 3D environment in either a mesenchymal mode, which depends on integrin-mediated adhesion and extracellular protease activity, or in an amoeboid mode that is independent of integrins and localised proteolysis, and is also characterised by largely dissociated microtubules and high RhoA-ROCK-dependent contractility [36][37][38]. As RhoA-ROCK signalling is also required for mesenchymal invasion, it is an attractive target for antimetastatic therapy [39,40]. A375M2 metastatic melanoma cells have been a popular model for studying amoeboid invasion and its regulation, for they show prominent blebbing amoeboid morphology in, or on top of, soft collagen matrices, good long-term viability under such conditions, and, importantly, display substantial plasticity in terms of induced migratory phenotype transitions between amoeboid and mesenchymal modes, referred to as invasion plasticity. So far, amoeboid-mesenchymal transition (AMT) in A375M2 has been achieved by chemical inhibition of ROCK kinases, or gene targeting or knockdowns of miscellaneous regula-tors involved in the underlying signalling network [29,30,[41][42][43][44][45]. Here, we show, for the first time, that treatment of A375M2 cells embedded in a 3D collagen matrix with compounds SB202190 and BIRB796 induces a phenotype switch, as documented by decreased expression of dedifferentiated/invasive phenotype-associated genes and upregulation of differentiated/proliferative phenotype-associated genes. The differentiation phenotype switch is accompanied by a morphological shift from a rounded to an elongated shape, which is indicative of AMT, suggesting that melanoma phenotype plasticity might be interconnected with invasion plasticity. SB202190 is more potent than BIRB796 in all observed effects, particularly in the induction of differentiation-related markers. While both compounds have been developed and used as p38 MAPK inhibitors, they can significantly inhibit several p38-unrelated kinases (see Table 1) [46,47]. Some of these kinases may affect melanoma cell phenotype, as exemplified by SB202190-specific inhibition of BRAF [48]. Recently, an important property distinguishing SB202190 from BIRB796 was identified: a SB202190-specific, p38-independent ability to robustly translocate lysosomal membranebound TFEB, a MITF-related protein, to the nucleus of melanoma cells [48]. The same effect was observed for TFEB and TFE3 (another MITF family member) in non-melanoma cells [49]. We confirmed that SB202190, but not BIRB796, elicits nuclear translocation of the MITF-A isoform in melanoma cells (see Figshare file, "Subcellular localisation of MITF" [50]). While MITF-M, the major isoform in melanoma/melanocytes, is constitutively localised in the nucleus, the localisation is not regulated by mTORC1, unlike the other MITF isoforms and TFEB, and does not change after SB202190 or BIRB796 treatment; it forms heterodimers with the other isoforms and with TFEB, TFEC and TFE3 [51,52]. Massive nuclear translocation of these proteins due to SB202190 treatment thus might potentiate transactivation by MITF-M or replace it in the transactivating complexes. Indeed, MITF-A was found to be a more potent activator of DCT than MITF-M, and TFE3 may replace MITF-M in the complex with LEF-1 to activate the DCT promoter [53,54]. x/+ p38delta x JNK2 x/+ xx/++ JNK3 x RSK2 x BRSK2 x

Cells Culture, Treatments and Morphological Analysis
A375M2 cells were maintained in DMEM (4.5 g/L glucose, pyruvate) supplemented with 10% foetal bovine serum and 50 µg/mL gentamicin (all from Sigma-Aldrich, St. Louis, MO, USA) at 37 • C in a humidified atmosphere with 5% CO 2 . The cultures were regularly tested for mycoplasma contamination. 3D cell culture experiments were performed with rat tail collagen (SERVA, Heidelberg, Germany) at a concentration of 1 mg/mL and DMEM supplemented with 1% foetal bovine serum, 15mM HEPES and 50 µg/mL gentamicin. In all experiments, cells were treated with 10 µM SB202190 (Santa Cruz Biotechnology, Santa Cruz, CA, USA), 10 µM BIRB796 alias doramapimod (LC Laboratories), or an equivalent volume of DMSO only (0.1% final concentration). To quantify cell morphology in the 3D collagen matrix, 100,000 cells were seeded in 250 µL of collagen in a 48-well plate. Hoffman modulation contrast microscopy images were taken 48 h later from approx. 20 planes along the z-axis. Cells were classified as elongated or rounded as described previously [55]. Please note that different A375-derived cell lines with almost identical names (A375M2, A375 M-2, A375M, etc.) were created in similar ways in a few completely independent studies [56,57]. A375M2 cells used in this study are those from Clark et al., and were kindly provided by Prof. Richard Hynes.

RNA Extraction and Sequencing
A total of 330,000 A375M2 cells were cultured in 500 µL of 3D collagen gel for 48 h in a 24-well plate. RNA was isolated as described in detail previously [55]. Polyadenylated RNA was enriched with Oligo d(T) 25

Reverse Transcription-Quantitative Polymerase Chain Reaction (RT-qPCR)
All the RT-qPCR experiments were performed according to MIQE guidelines [61]. Briefly, total RNA was extracted from cells embedded in 3D collagen matrix as described above. RNA reverse transcription was performed using M-MuLV Reverse Transcriptase (NEB) with 1 µM oligo(dT) and (dN) 15 primer and 1.5 µg of total RNA. SYBR green-based qPCR was performed in a LightCycler 480 Instrument (Roche). For primer details, see the file, PCR.primers.xlsx, available from Figshare [50]. Cq and relative expression values were calculated in R with ReadqPCR and NormqPCR packages, using TARDBP and THRAP3 as reference gene indices (selected based on geNORM analysis) [62,63]. Significance of differences was analysed with two-way ANOVA and pairwise paired t-tests, with p-value adjustment for multiple testing (Benjamini-Hochberg [64]).

Statistical Analysis
To estimate differential gene expression from RNA sequencing data, a workflow based on the STAR aligner and DESeq2 R package was used, as described previously [65]. We used rlog-transformed values for principal component analysis. Gene set enrichment analysis was performed with Enrichr [66]. We utilised core R functions to compare data produced in this work against two public datasets (Verfaillie et al. 2015 andFolberg et al. 2006) [11,67]. We calculated the numbers of shared genes (represented by gene symbols) for all possible pairs of gene sets, and the related Jaccard indices and hypergeometric right-tail cumulative p-values, to estimate the statistical significance of the intersections. As an approximation, we used the total number of all HGNC gene symbols assigned to protein-coding and non-coding RNA genes (27,371) as the size of the source population in hypergeometric tests. Venn diagrams representing intersections of the gene sets were generated using a Multiple List Comparator online tool [68].

SB202190 and BIRB796 Induce a Phenotype Switch in A375M2 Cells Cultured in 3D Collagen
Treatment of A375M2 melanoma cells with 10 µM compounds SB202190 and BIRB796 for 48 h in 3D collagen resulted in an elongated shape in the majority of the cell population (Figure 2A,C), also see Figshare item, "Bright field microscopy photographs of A375M2 3D cell cultures"). The compounds were well tolerated, with only a mild decline in cell numbers after SB202190 treatment ( Figure 2B). RT-qPCR analysis of transcripts associated with the dedifferentiated/invasive phenotype (INHBA, PODXL, IL1A, IL1B), and of MITFregulated genes (DCT, MLANA) associated with the differentiated/proliferative phenotype, confirmed a phenotype switch consistent with the characteristic gene expression signature described in the literature ( Figure 2D). SB202190 was more potent than BIRB796 in both morphological and gene expression level effects.

Transcriptomic Profiling of the Phenotype Switch-Associated Changes with RNA-seq
To uncover the transcriptomic changes underlying the dual phenotype switch, we performed RNA sequencing of mRNA samples from SB202190, BIRB796 or DMSO-treated cells prepared in four independent biological replicates. The cells were kept for 48 h in 3D

Transcriptomic Profiling of the Phenotype Switch-Associated Changes with RNA-seq
To uncover the transcriptomic changes underlying the dual phenotype switch, we performed RNA sequencing of mRNA samples from SB202190, BIRB796 or DMSO-treated cells prepared in four independent biological replicates. The cells were kept for 48 h in 3D collagen with the compounds, and then the whole samples, including the collagen, were homogenised and further processed for RNA sequencing. Polyadenylated RNA was enriched with oligo-dT magnetic beads, converted into a stranded cDNA library, and sequenced with an Illumina HiSeq sequencer (Figure 3). The sequencing yielded approx. 21-26.5 million paired-end reads per sample. Raw reads trimmed from adapter sequences were quality checked with FastQC software [69]. FastQC's plots of Phred scores by position showed typical profiles with decreasing quality towards the ends of reads ( Figure 4A). The sequencing yielded approx. 21-26.5 million paired-end reads per sample. Raw reads trimmed from adapter sequences were quality checked with FastQC software [69]. FastQC's plots of Phred scores by position showed typical profiles with decreasing quality towards the ends of reads ( Figure 4A). Highlighted in red are the genes whose higher expression is well known to be characteristic of either the differentiated/proliferative (E) or of the dedifferentiated/invasive phenotype (F).
The mapping metrics of the STAR aligner showed an average paired read length of 240 nucleotides. On average, the aligner uniquely mapped 93.6% of the fragments to the human genome; 3.1% of the fragments were multi-mapped and 3.1% of the fragments were excluded for being too short. For the complete mapping metrics, see Figshare file, "RNA sequencing raw data and mapping metrics" [50]. Principal component analysis of normalised gene expression profiles showed that the treatments were the dominant source of variation ( Figure 4B). Differential gene expression analysis with DESeq2 [70] identified that 2700 out of 16,612 Ensembl genes (16.3%) changed with adjusted p-value < 0.1 and fold change > 1.5 after treatment with SB202190. BIRB796 treatment changed the expression of 1638 out of 16,612 (9.9%) genes. The distribution of affected transcripts by gene expression level is depicted as an MA plot in Figure 4C. Complete results are The mapping metrics of the STAR aligner showed an average paired read length of 240 nucleotides. On average, the aligner uniquely mapped 93.6% of the fragments to the human genome; 3.1% of the fragments were multi-mapped and 3.1% of the fragments were excluded for being too short. For the complete mapping metrics, see Figshare file, "RNA sequencing raw data and mapping metrics" [50]. Principal component analysis of normalised gene expression profiles showed that the treatments were the dominant source of variation ( Figure 4B). Differential gene expression analysis with DESeq2 [70] identified that 2700 out of 16,612 Ensembl genes (16.3%) changed with adjusted p-value < 0.1 and fold change > 1.5 after treatment with SB202190. BIRB796 treatment changed the expression of 1638 out of 16,612 (9.9%) genes. The distribution of affected transcripts by gene expression level is depicted as an MA plot in Figure 4C. Complete results are available from Figshare ("Differential gene expression analysis") [50]. Selected differentially expressed genes from the gene lists, characteristic of either invasive or proliferative phenotypes (see Figure 1), are listed in Table 2. Next, we analysed the differentially expressed genes for over-representation of gene lists from published datasets using the Enrichr web server (Table 3). We found that the genes upregulated by SB202190 and BIRB796 are significantly enriched with genes with high expression in a subset of melanoma cell lines from the NCI-60 cancer cell line panel. Importantly, all the cell lines in this subset were previously classified as "proliferative" within a large panel of melanoma cell lines, see the HOPP database [71]. The genes downregulated by the compounds showed many more concordant enrichments; Table 3 displays the top enriched GO Biological Process terms and TRRUST transcriptional regulatory networks [72]. For complete results, see the Figshare file, "Gene set enrichments".

Validation and Reproducibility of the RNA-seq Results
To validate the biological reproducibility of the results of differential gene expression analysis, we performed RT-qPCR measurements of three selected transcripts-CEMIP, S100A16 and STC1-in RNA samples independent from those used for the RNA sequencing (three biological replicates). The corresponding genes were identified as significantly differentially expressed in both SB202190 and BIRB796 data, and are also present in published gene expression signatures of the invasive phenotype, e.g., S100A16 in Verfaillie et al. 2015, Tirosh et al. 2016and Jeffs et al. 2009. The results confirmed that the genes were similarly affected in the independent samples ( Figure 4D). We also wondered to what extent are the data generated in the presented study matching the published gene expression signatures characterising the phenotype switch. For comparison, we chose two published gene expression profiles of invasive and proliferative phenotypes with a total number of affected genes similar to our study- Verfaillie et al. 2015 andFolberg et al. 2006 [11,67]. The comparison revealed that the two published studies overlap, with the data from SB202190-treated A375M2 cells being similar. Data from BIRB796-treated cells overlap by approx. half of the numbers, but still significantly (see Figshare data, "Comparison of A375M2 data sets against public datasets" [50]). We also analysed all intersections of the gene lists in detail and found subsets shared by all the lists (11 genes for "proliferative" and 25 genes for "invasive" sets) containing known markers of the respective phenotypes ( Figure 4E,F).

Discussion
The signalling network controlling amoeboid phenotype in melanoma, as has been previously uncovered in A375M2 cells, is virtually identical to the one responsible for the dedifferentiated/invasive phenotype (i.e., NF-κB, TGFβ and non-canonical WNT signalling) [41,42,75]. The amoeboid phenotype is also known to be suppressed by βcatenin-stimulated MITF expression [76]. EPHA2, which is highly expressed in dedifferentiated/invasive melanoma cells, stimulates RhoA activity and promotes amoeboid migration [77]. It thus seems that amoeboid invasion is a manifestation of the dedifferentiated/invasive phenotype observable under suitable 3D conditions (in terms of cell viability and matrix composition). Our data bring essential additional support for such a hypothesis and encourage further large-scale 3D culture-based or intravital analyses required for definitive confirmation. We suggest that the presented data are of high relevance to melanoma biology and cancer metastasis research fields, as they bridge invasion and phenotype plasticity in melanoma cells.  Data Availability Statement: Supporting data referenced in the text, complete results of differential gene expression analysis, and published data comparisons have been deposited in the Figshare repository [50]. Adapter-trimmed RNA sequencing data have been deposited in the ArrayExpress database at EMBL-EBI (accession number E-MTAB-9273) [60]. The scripts used in data processing are available from Figshare [50] (see "Code used in data processing").