Differential Gene Expression Common to Acquired and Intrinsic Resistance to BRAF Inhibitor Revealed by RNA-Seq Analysis.

Melanoma cells have been shown to respond to BRAF inhibitors; however, intrinsic and acquired resistance limits their clinical application. In this study, we performed RNA-Seq analysis with BRAF inhibitor-sensitive (A375P) and -resistant (A375P/Mdr with acquired resistance and SK-MEL-2 with intrinsic resistance) melanoma cell lines, to reveal the genes and pathways potentially involved in intrinsic and acquired resistance to BRAF inhibitors. A total of 546 differentially expressed genes (DEGs), including 239 up-regulated and 307 down-regulated genes, were identified in both intrinsic and acquired resistant cells. Gene ontology (GO) analysis revealed that the top 10 biological processes associated with these genes included angiogenesis, immune response, cell adhesion, antigen processing and presentation, extracellular matrix organization, osteoblast differentiation, collagen catabolic process, viral entry into host cell, cell migration, and positive regulation of protein kinase B signaling. In addition, using the PANTHER GO classification system, we showed that the highest enriched GOs targeted by the 546 DEGs were responses to cellular processes (ontology: biological process), binding (ontology: molecular function), and cell subcellular localization (ontology: cellular component). Ingenuity pathway analysis (IPA) network analysis showed a network that was common to two BRAF inhibitor-resistant cells. Taken together, the present study may provide a useful platform to further reveal biological processes associated with BRAF inhibitor resistance, and present areas for therapeutic tool development to overcome BRAF inhibitor resistance.

INTRODUCTION

Activating mutations in the BRAF gene are observed in around 50% of melanomas (Davies ; Cancer Genome Atlas Network, 2015). The detection of BRAF mutations initiated great scientific efforts to develop drugs for the treatment of metastatic melanoma harboring BRAF mutations (Holderfield ). Vemurafenib (PLX4032) was the first drug developed to specifically target the BRAF V600E mutation (Bollag ). We also previously reported UAI-201 (also known as UI-152) as a potent ATP-competitive inhibitor of RAF proteins (Kim ). UAI-201 is more than 1,000- fold more selective at inhibiting the proliferation of tumor cell lines bearing the BRAF V600E mutation when compared with that of cells carrying wild-type (WT) BRAF (Kim ).

In most instances, tumor cells gradually develop a loss of response to anticancer drugs, acquiring resistance during the course of treatment (Dobbelstein and Moll, 2014). BRAF inhibitors also show poor clinical responses due to acquired resistance of melanoma cells (Rizos ). Therefore, understanding the molecular mechanisms of drug resistance is necessary to improve the efficacy of BRAF inhibitors (Arozarena and Wellbrock, 2017). Although multiple BRAF inhibitor resistance mechanisms have been described, most studies have shown that the majority of resistances are centered on reactivation of the MAPK-pathway (Johannessen ; Villanueva ; Shi ). Interestingly, a recent study reported that U34 enzymes, which catalyze modifications of wobble uridine 34 tRNA, promote the survival and resistance to therapy of melanoma cells by regulating specific mRNA translation (Rapino ). Our previous study indicated that induction of resistance to BRAF inhibitors is associated with the inability of SPRY2 to inhibit BRAF-V600E activity in cells expressing mutant BRAF (Ahn ). In addition, our previous miRNA microarray profiling study revealed that ectopically expressed miR-1246 can confer resistance to BRAF inhibitors (Kim ). On the contrary, melanoma expressing WT BRAF present intrinsic resistance and do not respond to selective BRAF inhibitors (Flaherty ; Kim ), which also limits the effectiveness of this therapeutic approach. Rather, secondary cancers have subsequently found to be induced by the paradoxical activation of MAPK signaling by selective BRAF inhibitors in cells harboring WT BRAF, but activated RAS (Halaban ; Hatzivassiliou ; Heidorn ). Several authors have recently reported the development of a second generation of BRAF inhibitors, so-called paradox breakers that do not exhibit paradoxical ERK activation (Zhang ; Agianian and Gavathiotis, 2018).

Several published studies have identified gene markers of BRAF inhibitor resistance (Konieczkowski ; Ji ; Verfaillie ); however, most of these studies have focused on acquired resistance. In this study, we performed RNA sequencing with Illumina HiSeq 2500 technology on three melanoma cell lines: BRAF inhibitor-sensitive A375P BRAF V600E cells, their BRAF inhibitor-resistant counterparts (A375P/Mdr), and SK-MEL-2 cells with intrinsic resistance to BRAF inhibition, to reveal the potential genes and pathways involved in intrinsic and acquired resistance to BRAF inhibitors. We show that 546 DEGs are deregulated in two BRAF inhibitor-resistant cells compared with BRAF inhibitor-sensitive cells. Differential gene expression and pathway analysis confirmed that resistance to BRAF inhibitors is not triggered by one common mechanism but involves several functional pathways.

MATERIALS AND METHODS

Cell lines and culture

We used three melanoma cell lines: BRAF inhibitor-sensitive A375P BRAF V600E cells, their BRAF inhibitor-resistant counterparts (A375P/Mdr), and SK-MEL-2 BRAF-WT cells with intrinsic resistance to BRAF inhibitors. A375P and SK-MEL-2 cells were acquired from either the Korean Cell Line Bank (KCLB; Seoul, Korea) or YOUAI Co., Ltd (Suwon, Korea). A375P/Mdr cells were previously established by chronic selection with increasing doses of an inhibitor of oncogenic BRAF (Ahn and Lee, 2013). The SK-MEL-2 cell line expressing WT BRAF has intrinsic resistance to BRAF inhibition because the BRAF inhibitor lacks activity against cell lines that express WT BRAF. All cell lines were maintained at 37°C in DMEM supplemented with 10% FBS (Thermo Fisher Scientific, Carlsbad, CA, USA), penicillin–streptomycin, and glutamine (Thermo Fisher Scientific). A375P/Mdr cells were further propagated in a growth medium containing the BRAF inhibitor PLX4720 (1 μM, Selleck Chemicals, Houston, TX, USA). Before use in experiments, A375P/Mdr cells were maintained in a PLX4720-free culture medium and subcultured at least three times. For experimental purposes, cells were cultured in 60-mm tissue culture dishes until they reached ∼80% confluence. PLX4720 was dissolved in dimethyl sulfoxide (DMSO) and was diluted immediately prior to each experiment. The final DMSO concentration was less than 0.1% in all experiments.

RNA preparation

Total RNA was prepared using the RNeasy Midi Kit (Qiagene, Valencia, CA, USA) from A375P, A375P/Mdr, and SK-MEL-2 cells grown under preferred culture conditions as described above. Sample integrity was verified using a Nano-drop 1000 and Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).

Illumina sequencing

The cDNA libraries were created using RNA isolated from three melanoma cell lines using the Illumina TruSeq RNA Library Preparation Kit (Illumina Inc., San Diego, CA, USA) following the manufacturers protocol. Briefly, PolyA selection was performed using RNA purification beads on 0.1–4 μg of total RNA per sample to obtain mRNA. For each RNA sample, 10–400 ng of purified mRNA was fragmented and used to synthesize first-strand cDNA using reverse transcriptase and random hexamer primers. Second-strand cDNA synthesis was performed using dNTPs including dUTP, DNA polymerase, and RNase. Adding dATP to all free 3′ ends and adaptor ligation added unique index sequences to the fragments performed end repair. The products were then amplified by PCR and purified. The cDNA libraries were quantified using a Qubit 2.0 fluorometer (Invitrogen) and assessed for correct fragment size (∼260 bp) using an Agilent Bioanalyzer 2100 (Agilent Technologies). Libraries were normalized to 10 nM, pooled, denatured, and sequenced on the Hiseq2500 Sequencing platform (Illumina Inc.) using a 101 bp paired-end approach targeting approximately 5000 million reads per sample.

RNA-Seq analysis

Paired end sequencing (100 bp) was performed using an Illumina HiSeq 2500 instrument. All the 100 bp paired end raw reads were quality checked for low quality bases and adapter sequences (processing of raw reads). Quality check was performed using FastQC, followed by trimming at Q20. Overall alignment rate was higher than 97% for all samples. Low quality bases and reads were filtered from the datasets prior to read mapping. Cutadapt (v1.9.1) (Martin and Wang, 2011) and Sickle (v1.33) were used to remove low quality reads and adaptor contamination. We used HISAT2 (v2.0.5) to map sequence reads to UCSC (the University of California Santa) reference human genome hg19.

Analysis of differentially expressed genes

The quantification of gene expression was estimated Fragments Per Kilobase Million (FPKM). Statistically significant expression changes between A375P and A375P/Mdr cells were estimated using Cufflinks v2.2.1 (cuffnorm). We determined significantly expressed genes in each sample as having a log2 fold change at ≥1.

Functional gene ontology

DEGs with a log2 fold change ≥1, as identified by Cuffdiff, were analyzed to identify enriched biological processes based on gene ontology (GO). This was performed using ClueGO v2.3.3, a plugin of the software platform Cytoscape v3.4.0-RC1 (Bindea ). Processes that had p-value≤.01 were considered significant and displayed.

Pathway analysis

Identification of biological pathways involving significant DEGs, was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 (Huang , 2009b). Pathways containing ≥2 significant DEGs were categorized and displayed.

Network analysis by Ingenuity Pathway Analysis software

Network analysis was performed with Ingenuity Pathway Analysis (IPA) software (http://www.ingenuity.com/) to increase confidence in the observations of differentially expressed genes (DEGs) by correlation with biological pathways. The list of DEGs in two BRAF inhibitor-resistant cell lines, containing gene identifiers and corresponding expression values, was uploaded to the IPA software (Qiagen). IPA computes a score for each network according to the fit of the set of supplied focus genes. The network interaction of the focused genes in the network is based on their connectivity in Ingenuity Knowledge Base. These scores indicated the likelihood of focus genes belonging to a network, versus those obtained by chance.

RESULTS

Comparative transcriptome profiling

We used a Next Generation Sequencing (NGS) approach with Illumina HiSeq 2500 technology as a means of examining and comparing the transcriptomes of BRAF inhibitor-sensitive (A375P) and -resistant (A375P/Mdr and SK-MEL-2) melanoma cell lines. These cell lines have been previously characterized with respect to BRAF inhibitor resistance (Ahn and Lee, 2013). For the FPKM normalization method, the aligned RNA-Seq reads were normalized by FPKM using Cufflinks-Cuffnorm (v2.2.1). Differential gene expression analysis revealed several genes that were differentially regulated in BRAF inhibitor-resistant A375P/Mdr cells when compared with BRAF inhibitor-sensitive A375P cells with log2 fold change at ≥1 and q-value <0.05. From a total of 25,271 genes, 1,931 significant DEGs were identified in A375P/Mdr cells when compared with A375P cells, of which 887 were up-regulated and 1044 down-regulated (Supplementary Table 1). To provide a more detailed overview of the expression patterns of genes responsive to BRAF inhibitor resistance, a heat-map of the hierarchical clustering was generated based on transcriptome data from three melanoma cell lines for the 1,931 significant DEGs identified in A375P/Mdr cells. This analysis classified the DEGs into five groups with distinct expression profiles (Fig. 1). In particular, DEGs in group-I (n=486), which are mainly involved in tissue development, regulation of cell adhesion, anatomical structural morphogenesis, positive regulation of transcription, and the cell cycle, were repressed in two BRAF inhibitor-resistant cell lines. Expression of group V genes (n=123) was significantly induced in two BRAF inhibitor-resistant cell lines. This group includes genes encoding cellular response to the lipopolysaccharide, ovulation cycle process, and regulation of interleukin-1 beta production, as well as those involved in regulation of epithelial cell apoptotic processes. Group II and IV contained down- and up-regulated unique DEGs, respectively, only in A375P/Mdr cells. Group II (n=357) included those genes related to the cell cycle, and cellular responses to DNA damage. The largest, group-IV DEGs (n=708) were mainly related to the inflammatory response and cell migration. DEGs for group III (n=257), induced only in SK-MEL-2 cells, included genes encoding chromosome organization mechanisms, and DNA repair. Interestingly, the pattern of gene expression changes in SK-MEL-2 with intrinsic resistance was similar to that of parental A375P rather than A375P/Mdr with acquired resistance.

Fig. 1.

Hierarchical clustering and heat map analysis of differentially expressed genes. Heat maps showing results of gene expression clustering between BRAF inhibitor-sensitive (A375P) and resistant cells (A375P/Mdr and SK-MEL-2). Clustering groups are represented by different colors above each heat map. Red bars indicate relatively high signal intensity for the gene product in question, with blue representing lower intensity and grey intermediate. Clustering of the samples is indicated by the dendrogram on top.

Functional GO analysis

We focused our attention on sequences that were differentially expressed in common (>2-fold change and p<0.05) between A375P/Mdr and SK-MEL-2 vs. A375P cells, regardless of whether the resistance was intrinsic or acquired. In SK-MEL-2 with intrinsic resistance, 5,727 significant DEGs with a log2 fold change of ≥2 were detected out of which 2,871 were up-regulated and 2,855 down-regulated when compared to A375P cells (Supplementary Table 2). Among the differentially regulated genes, 546 genes were found common across A375P/Mdr and SK-MEL-2 (Supplementary Table 3). Of these, 239 genes were up regulated and the remaining 307 down regulated (Fig. 2A). GO analysis of 546 DEGs common in A375P/Mdr and SK-MEL-2 identified ten processes including angiogenesis, immune response, cell adhesion, extracellular matrix organization, collagen catabolic process, cell migration and positive regulation of protein kinase B signaling (Fig. 2B). Table 1 lists the top ten genes commonly up- and down-regulated in two BRAF inhibitor-resistant cell lines compared to A375P cells and the genes of interest.

Fig. 2.

Pathway and Gene ontology (GO) analysis of DEGs common between two BRAF inhibitor-resistant cells. (A) Venn diagram showing unique as well as common DEGs in A375P/Mdr and SK-MEL-2 cells as compared to A375P cells. (B) Ten significant pathways associated with DEGs common in BRAF inhibitor-resistant cells. Significant DEGs with a log2 fold change ≥1 were analyzed for enriched biological processes using the ClueGO application for Cytoscape. The vertical axis represents the pathway category, and the horizontal axis represents the −log2 (p-value) of the pathways. Greater −log2 (p-value) scores correlate with increased statistical significance.

Table 1.

The genes of interest and top 10 DEGs common to two BRAF inhibitor-resistant melanoma cells based on log2 fold change values

Genes	Functions	Fold-change
		A375P/Mdr	SK-MEL-2
Up-regulated Genes in in BRAF inhibitor-resistant cells
SERPINE1	Serpin peptidase inhibitor, clade E	50.19	5.52
BASP1	Brain abundant, membrane attached signal protein 1	45.50	22.02
IGFN1	Immunoglobulin-like and fibronectin type III domain	35.72	12.66
C15orf48	Chromosome 15 open reading frame 48	33.01	11.66
CCL3	Chemokine (C-C motif) ligand 3	19.14	3.07
KRT75	Keratin 75, type II	18.39	3.59
ACSL5	Acyl-CoA synthetase long-chain family member 5	17.03	16.92
NUPR1	Nuclear protein, transcriptional regulator, 1	16.71	44.19
SERPINB2	Serpin peptidase inhibitor, clade B (ovalbumin)	15.97	11.12
KRT34	Keratin 34, type I	15.76	2.56
CASP1	Caspase 1, apoptosis-related cysteine peptidase	13.20	38.94
MAGEB2	Melanoma antigen family B2	8.03	61.23
JUN	Jun proto-oncogene	5.50	2.43
MCAM	Melanoma cell adhesion molecule	2.56	2.63
MOAP1	Modulator of apoptosis 1	2.16	3.10
Down-regulated Genes in in BRAF inhibitor-resistant cells
COL9A3	Collagen, type IX, alpha 3	−91.97	−2.52
TIMP3	TIMP metallopeptidase inhibitor 3	−77.28	−9.02
CD96	CD96 molecule	−52.94	−4.72
SPRR2D	Small proline-rich protein 2D	−33.19	−28.37
NES	Nestin	−29.38	−6.56
S100A4	S100 calcium binding protein A4	−25.79	−17.49
HLA-DRA	Major histocompatibility complex, class II, DR alpha	−23.68	−135.06
FST	Follistatin	−17.69	−75.48
SLITRK6	SLIT and NTRK-like family, member 6	−17.72	−30.46
DUSP9	Dual specificity phosphatase 9	−2.97	−3.17
SPRY2	Sprouty homolog 2 (Drosophila)	−2.89	−2.29

Pathway analysis of DEGs

To further understand their biological functions and pathways, the 546 DEGs common to the two BRAF inhibitor-resistant cell lines were classified according to GO terms at the biological process, molecular function, and cellular compartment levels, respectively, using the PANTHER GO classification system (Ashburner ). The DEGs were correlated with a wide range of biological processes and molecular functions. In respect to molecular function (Fig. 3A), the most common differential expression proteins were associated with binding, followed by catalytic activity, receptor activity, signal transducer activity, transporter activity, and structural molecule activity. With regard to biological processes (Fig. 3B), the majority of proteins were involved in cellular processes (27.1%), metabolic processes (19.6%), biological regulation (11.3%), response to stimulus (10.1%), developmental processes (8%), multicellular organismal processes (5.8%), cellular component organization or biogenesis (5.6%), and localization (4.7%). With respect to cellular component (Fig. 3C), DEGs were cell part proteins (36%), organelle proteins (24.3%), macromolecular complex proteins (10.3%), membrane proteins (13.1%), extracellular region proteins (3.5%), and extracellular matrix proteins (3.5%).

Fig. 3.

Gene enrichment and pathway analysis. A total of 546 common DEGs were subjected to gene ontology (GO) classification using the PANTHER gene analysis tool. The three main GO categories include biological process, cellular component, and molecular function.

Network analysis by IPA

Finally, gene network analysis was performed with IPA to look at the associations of the 546 DEGs common between BRAF inhibitor-resistant cell lines by performing direct interaction analysis, and with other genes by indirect interaction analysis. We identified 25 networks (p-scores>18) in BRAF inhibitorsensitive vs. -resistant cells (Supplementary Table 4). The two significant networks, which were considered to be associated with drug resistance, were ‘cellular movement, cancer, cellular response to therapeutics’ (Fig. 4A) and ‘cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation’ (Fig. 4B). Table 2 includes the molecules, score, and focus molecules in two networks. There were 20 genes involved in the down-regulated network of ‘cellular movement, cancer, cellular response to therapeutics’, and 16 genes in the network of ‘cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation’. Interestingly, the number of down-regulated genes was twice that of the up-regulated genes (down-regulated:up-regulated=13:7) in ‘cellular movement, cancer, cellular response to therapeutics’, while the opposite was true (down-regulated:up-regulated=5:11) for ‘cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation’.

Fig. 4.

IPA-identified two significant gene networks. (A) Cellular movement, cancer, cellular response to therapeutics. (B) Cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation. The lines between genes represent known interactions (solid lines represent direct interactions; dashed lines represent indirect interactions). Each gene is displayed using shapes that represent the functional class of the gene product, as indicated in the legend. A total of 25 networks (p-scores>18) are provided in Supplementary Table 2.

Table 2.

Networks predicted by IPA in BRF inhibitor-resistant cells

Top diseases and functions	Molecules in network	Score	Focus molecules
Cellular Movement, Cancer, Cellular Response to Therapeutics	14-3-3, Alpha catenin, BCAM, Cadherin, CaMKII, Creb, CXADR, DIO2, EGR3, EPHA3, FSCN1, GLI2, Gsk3, Hdac, IGF1R, Mlc, MME, PDGF BB, PI3K (family), Pka catalytic subunit, PP2A, Ppp2c, QPCT, Rap1, RTL1, RUNX2, S100A4, SLC26A2, SNAI1, SNAI2, SORCS1, STIM2, TCF, TP73, TSPAN1	24	20
Cell Death and Survival, Cell-To-Cell Signaling and Interaction, Cellular Growth and Proliferation	CASP1, chemokine, COL19A1, COL9A3, collagen, Collagen type ix, Complement, cytochrome C, FN1, Growth hormone, Hedgehog, Hsp27, IFIH1, IFIT3, Ifn, IFN Beta, Ifn gamma, IL1, IL18, Interferon-α Induced, Laminin (family), MAC, Mmp, OASL, PAEP, PRB1/PRB2, Pro-inflammatory Cytokine, PTGFR, SERPINB2, SLC5A3, Tlr, TLR4, Tnf (family), TNFAIP3, UBA7	18	16

Red and blue denote up-regulated and down-regulated genes among the differentially expressed genes.

DISCUSSION

Intrinsic and acquired resistance limits the therapeutic benefits of oncogenic BRAF inhibitors in melanoma. Although several studies have reported the specific DEGs involved in acquired resistance to BRAF inhibitor (Konieczkowski ; Ji ; Verfaillie ), little is known about the genes associated with intrinsic resistance. Thus, we investigated genes whose expression levels changed in common between intrinsic resistant cells and acquired resistance cells. In the present study, we used an RNA-Seq approach with Illumina HiSeq 2500 technology as a means of examining and comparing the transcriptomes of BRAF inhibitor-sensitive (A375P) and -resistant (A375P/Mdr with acquired resistance and SK-MEL-2 with intrinsic resistance) melanoma cell lines. A heat-map of hierarchical clustering of 1,931 significant DEGs with a log2 fold change of ≥2 in A375P/Mdr cells was generated to visualize differential gene expression status between BRAF inhibitor-sensitive cells and -resistant cells. The analysis classified 1,931 DEGs into five groups with distinct expression profiles. In particular, groups I and V contained down- and up-regulated DEGs, respectively, common to A375P/Mdr and SK-MEL-2 cells. Down-regulated DEGs in group-I (n=486) are mainly involved in tissue development, regulation of cell adhesion, anatomical structural morphogenesis, positive regulation of transcription, and cell cycle, while up-regulated group V DEGs (n=123), include genes encoding cellular response to lipopolysaccharide, the ovulation cycle process, and regulation of interleukin-1 beta production, as well as those involved in regulation of epithelial cell apoptotic processes.

Among differentially regulated genes, a total of 546 DEGs with change >2-fold were found common between A375P/Mdr and SK-MEL-2, regardless of whether resistance was intrinsic or acquired. Of these, 239 genes were up regulated and the remaining 307 down regulated. Some up-regulated genes, including caspase 1 and modulator of apoptosis 1, have been previously shown to be associated with apoptosis. Another up regulated gene, melanoma antigen family B2 (MAGEB2), has recently been implicated in carcinogenesis and identified a potential cancer biomarker (van Duin ; Pattani ). In particular, we observed that melanoma cell adhesion molecule (MCAM), which mediates chemoresistance in smallcell lung cancer via the PI3K/AKT/SOX2 signaling pathway (Tripathi ), was among the significantly up regulated genes. Therefore, we propose that the observed high levels of MCAM might be involved in increased resistance to BRAF inhibitors. A group of genes, which were significantly differentially down regulated in the BRAF inhibitor-resistant cells, included the dual specificity phosphatase 9 (DUSP9), and Sprouty 2 (SPRY2). It has been shown that DUSP9 contributes to chemotherapy-induced BCSC enrichment through inactivation of ERK and activation of p38 MAPK signaling pathways (Lu ). SPRY2 is a potent antagonist of the FGFR-induced Ras/MAPK signal pathway (Hacohen ). We previously reported that the enhancement of RAF-1 kinase activity by knockdown of SPRY2 is associated with high sensitivity to paclitaxel in v-Ha-ras-transformed NIH 3T3 fibroblasts (Ahn ), suggesting that up regulation of SPRY2 might significantly contribute to BRAF inhibitor resistance.

Using the PANTHER GO classification system, we found that the highest enriched Gos, targeted by the 546 DEGs, were responses to cellular process (ontology: biological process), binding (ontology: molecular function), and cellular organization (ontology: cellular component). GO classification by molecular function showed that DEGs were assigned 9 different functions (binding, structural molecule, antioxidant, signal transducer, receptor, channel regulator, catalytic activity, transregulator, and transporter). In addition, GO classification based on biological process included 13 biological groups; the 8 groups that contained the majority of the DEGs were metabolic process, multicellular organismal process, developmental process, response to stimulus, biological regulation, localization, cellular process, and cellular component organization or biogenesis. The GO classification by cellular component showed that DEGs were mapped to extracellular region, organelle, cell part, extracellular matrix, macromolecular complex, membrane, cell junction, and the synpase. On the contrary, KEGG pathway analysis of the common 546 DEGs enriched in ten pathways including angiogenesis, immune response, cell adhesion, extracellular matrix organization, collagen catabolic process, cell migration and positive regulation of protein kinase B signaling. In particular, angiogenesis is among the most affected.

The up- and down-regulated genes used for pathway analysis were also used for network analysis with IPA to identify associated biological networks. The two significant networks, considered to be associated with drug resistance, were ‘cellular movement, cancer, cellular response to therapeutics’ and ‘cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation’. Interestingly, the number of down-regulated genes was twice that of up regulated genes in the network of ‘cellular movement, cancer, cellular response to therapeutics’, while the opposite was true for the network ‘cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation’. These results imply that these differentially-expressed genes influence biological processes leading to the development of BRAF inhibitor resistance. In particular, IGF1R, RUNX2, GLI2, and FSCN1 act as vital crossroad regulators in the network of ‘cellular movement, cancer, cellular response to therapeutics’ network. Some of these genes along with various other genes like CXADR, SNAI2, TSPAN1, MME, RTL1, EGR3, FSCN1, and STIM3 also show indirect associations. Genes such as CASP1, TNFAIP3, TLR4, IL18, and FN1 might play an important role in the network of ‘cell death and survival, cell-to-cell signaling and interaction, cellular growth and proliferation’. This again allowed us to group subsets of these focus genes into functional networks associated with BRAF inhibitor resistance.

In conclusion, the present data shows, for the first time, 546 genes whose expression levels changed in common in intrinsic as well as acquired resistant cells. These DEGs are involved in various biological functions and provide a mechanism through which resistance to BRAF inhibitor may occur. Further studies are needed to identify the key pathways responsible for therapeutic escape from BRAF inhibitor treatment and to develop a promising strategy for preventing and abrogating the onset of both intrinsic and acquired drug resistance.