Transcriptomic Alterations Induced by Vemurafenib After Treatment of Melanoma: A Comprehensive Bioinformatics Analysis

Backgrounds Melanoma is a highly aggressive kind of cancer with very poor prognosis. B-raf inhibitor vemurafenib has indeed harvested substantial clinical benets. Nevertheless, its drug resistance has also hampered scientists effort towards successful melanoma treatment. In this study, we used data derived from the GEO database to analyze the effect on vemurafenib sensitive cell lines after vemurafenib treatment. GEO datasets GSE42872 (cohort1), GSE127988 (cohort2), GSE110054 (cohort3) were included in the analysis. Results We found 25 common Differentially Expressed Genes(DEGs) in 3 datasets, including 10 upregulated genes and 15 downregulated genes after vemurefenib application. Analysis using web tool TIMER showed signicant correlation of the upregulated genes with immune inltration level in skin cell melanoma. GO enrichment analysis showed that after vemurafenib treatment, all datasets showed downregulation in DNA replication and cell cycle arrest. Meanwhile, genes related to neuro-generation, extracellular matrix and cell-cell adhesion were signicantly enriched in all three datasets. KEGG analysis showed that pathways like P53, PI3K-Akt, and Rap signaling pathways were enriched in DEGs after vemurafenib administration. give


Background
Melanoma, a highly malignant tumor, has puzzled many scholars with its treatment [1].Approximately 1.7% of newly diagnosed primary malignancies worldwide are cases of cutaneous malignant melanoma which cause about 55,500 deaths annually [2]. Advanced melanoma is often associated with a lifethreatening tumor [3]. Before the advent of immune checkpoint inhibitors and gene-targeted drugs, there were few effective plans to treat advanced melanoma, and clinical trials did not yield satisfactory results [4]. PD1/PDL1 immunocheckpoint inhibitors have achieved good results in the treatment of melanoma, but there is no decrease in enthusiasm for the research of gene targeted therapy to improve the clinical effect [5]. B-RAF mutation, a genetic mutation found in 50% to 60% of melanomas, often occurs in the 600th codon, replacing valine by glutamic acid (so called B-RAF V600E ). This mutation leads to the activation of the BRAF/MEK/ERK pathway, causing continuous cell growth [6].
Vemurafenib, a novel drug rst synthesized in 2005, inhibits BRAF primarily. This drug was rst approved by the FDA in 2011 to treat melanoma with BRAF mutations [7]. However, patients who receive targeted therapies initially respond well to treatment, but there is resistance in the long term [8]. At present, researchers speculate that drug resistance may be related to secondary mutations, bypass signaling pathways, tumor microenvironment changes, and immune in ltration [9]. In order to provide supporting material for current drug resistance studies, we conducted a detailed bioinformatics analysis to explore what genetic changes were caused by vemurafenib in BRAF mutant melanoma cells.

Results
The ow chart The process of our study is summarized in Figure S1.
Characteristic BRAF overexpression and mutation in human melanoma tissues We enquired into the BRAF mutation and expression level within TCGA database using TCGA, cBioportal, and GEPIA. Firstly, a pan-cancer analysis showed that BRAF mRNA level in multiple melanomas was signi cantly higher than that in various other tumor tissues demonstrating that BRAF maybe a characteristic marker within multiple kinds of melanomas(Figure1A). Secondly, using TCGA database, we found that BRAF ranks third among the most frequently mutated genes in skin cutaneous melanoma (SKCM), with almost half of the total cases involved (Figure1B). Further analysis using GEPIA exhibited that BRAF expression level was signi cantly upregulated as melanoma stage upgrades ( Figure 1C).
Finally using TCGA datasets, we plotted the survival curve of the patients with and without BRAF mutation. The result was quite intriguing: the mutated group had a signi cant better outcome than the not mutated group, implying that treatment targeting BRAF mutation had signi cant effect, which is also the key point that we are going to address later on.
Expression Pro le before and after vemurafenib administration in vemurafenib sensitive cell lines In order to understand the changes in RNA expression after the application of vemurafenib cell lines, we used the GEO datasets GSE42872 (cohort1), GSE127988 (cohort2), GSE110054 (cohort3) to conduct the differential analysis using R. Among them, GSE127988 contains expression pro le data of 3 different vemurafenib concentrations: 0.1M, 0.316M, and1M. Therefore, 5 parallel experiments were included in the study. Using Venn diagram, we identi ed 25 common differential genes that co-existed in all those 5 parallel experiments (we deleted one invalid output from the 26 genes that served as an outcome Figure2A). The scatter plots of genetic expression were made ( Figure 3A,3C,3E). Meanwhile, a heatmap was drawn to illustrate the expression deviance of the 25 common differential genes before and after the treatment. These genes are potentially helpful in identifying a melanoma patient's response to vemurafenib treatment (Fig 3B,3D 3F). And the relationship between 10 highly expressed genes and immune in ltration was studied in TIMER, as shown in the Figure S4.
Differential genes were enriched in cellular metabolism, replication, extracellular adhesion, and certain tumorigenesis pathways.
We then used the R package 'ClusterPro ler' to further investigate where those differentially genes are signi cantly enriched. Not surprisingly, GO enrichment showed that the differential genes were predominantly enriched within biological process such as "DNA replication" "Cell Cycle" "DNA packaging", implying that after vemurafenib administration, the cellular replication activity was signi cantly reduced ( Figure 1A, C, E, FigureS1A). Moreover, cellular component enrichment analysis showed that those enriched genes were most frequently distributed within the chromosomal regions, as well as those regions that are extremely vital to cellular replication and metabolism, such as replication fork, spindle, and kinetochore. ( Figure S2C).Besides, molecular function GO analysis showed that enzymes that are essential in DNA replication and formation were also signi cantly enriched (FigureS2E).Interestingly, GO enrichment in all three GEO datasets showed that mRNAs concerned with extracellular matrix or cell-cell adhesion were signi cantly enriched, or to put it more precisely, upregulated ( Figure S2B,D,F). This nding is important in that it implies vemurafenib may enhance cellular adhesion via acting on extracellular matrix, thus reducing migration in melanoma cells, which also requires further investigation. Besides, we also found that molecular functions related to neurogenesis was signi cantly enriched in cohort 1 and 2( Figure S2B).KEGG analysis showed that multiple cellular signaling pathways were enriched after using the Braf inhibitor, such as the P53 signaling pathway was deactivated after using vemurafenib, while PI3K-AKT signaling pathways, Rap signaling pathways, as well as pathways related to focal adhesion were somehow activated after vemurafenib application( Figure 4B,4D,4F, FigureS3).

Protein-Protein interaction networks combined with TCGA survival analysis reveals certain transcriptionally regulated hub genes within melanoma cells
We rst drew a Venn diagram of the differentially expressed genes within the GEO datasets GSE42872, GSE 127988 (which contains 3 different vemurafenib concentration 0.1M,0.316M, and 1M), and GSE234231 ( Figure 2A). Not surprisingly, the differentially expressed genes have multiple genes in common. However, due to lack of su cient samples, the datasets of GSE110054 only had 98 differential genes which met our previous criteria. Therefore we selected the genes that are included in at least 4 parallel experiments, which is illustrated as samples encircled in red lines(Fig2B) .
By using the STRING database and the cytoscape software, we conducted a protein-protein interaction analysis of those common genes, (Figure5A). the result of our PPI analysis showed multiple hub genes that are at the core of the PPI network, which was further con rmed by survival analysis. Using GEPIA web tool, we examined the hub genes (meaning genes that have the maximal degree of freedom-interaction with other proteins) in our PPI network. We eventually found that the overexpression of most of the hub genes were related with poor survival outcomes. And most of those genes were genes that are essential in cell cycle and replication, for instance, gene CDK1 that encodes cyclin dependent kinase 1 (Figure5B). Apparently, most of the hub genes were signi cantly downregulated (Table1), showing the drug's core effect relies on melanoma cell cycle arrest.

Discussion
Targeted therapies for melanoma are evolving rapidly with the advance of basic research [12]. In the past, advanced melanoma therapy has ba ed many researchers, with an estimated 5-year survival rate of only 16 percent for patients with advanced malignant melanoma [13]. With the development of precision medicine, B-Raf mutations and immune checkpoint changes have been found in the pathogenesis of melanoma, which is an exciting discovery and has greatly promoted the treatment of melanoma [14]. This was followed by the emergence of B-Raf inhibitors (vemurafenib, dabrafenib) and immune checkpoint inhibitors (Ipilimumab, Pembrolizumab) that gave patients with advanced melanoma a ray of life [15]. A growing number of clinical trials have shown that these two type of drugs can signi cantly improve survival in patients with advanced melanoma [16].
In 2002, researchers found that melanoma cells often harbor BRAF mutations [14].This opens the door to targeted therapies. The activation process of this pathway is: growth factor binding receptor tyrosin kinases (RTKs) leads to automatic phosphorylation of the receptor, thus phosphorylation of RAS-GDP into RAS-GTP(active state). Activated RAS binds to RAF to activate downstream sites such as MEK/ERK [17,18,19].This cascade is highly conserved, with many protein-protein interaction networks that mediate cell proliferation, differentiation, and apoptosis [20].Normally, BRAF proteins are highly regulated by the body and are often bound to heat shock proteins(HSP) in a low active state [21]. Once the BRAF mutation occurs, it leads to the continuous activation of functional domains, causing a cascade of downstream sites, thereby promoting cell growth and inhibiting apoptosis [22,23].
Vemurafenib, which was primarily synthesized in 2005, targets BRAF effectively [24,25]. It was found to be selective for melanoma cells with BRAF mutation compared to non-BRAF mutant melanoma cells and has good clinical effect and strong compliance [24]. Vemurafenib had a combined response rate of about 50% and signi cantly improved progression-free survival compared with the previous chemotherapy regimen-dacarbazine [26]. Because BRAF mutations occur in 50 to 60 percent of melanomas, the drug could be a promising treatment for BRAF mutated melanomas [25].In addition, vemurafenib has shown promising results in the treatment of lung, colorectal, breast and ovarian cancer with BRAF mutations [27,28]. It mainly acts on the "DFG motif", which is continuously activated in BRAF mutations, to achieve inhibitory pathway activity [29,30]. Because of these bene ts, it was approved by the FDA in 2011 as the rst drug to treat cancers harboring BRAF mutations [31,32].
Despite all its virtues, it did not escape the fate of drug resistance [33]. Like other kinase inhibitors, vemurafenib shows striking clinical effects at rst, but resistance develops after 6-8 months [34]. Resurrection of mutations, contradictory activation of drugs, tumor microenvironment, and immune in ltration are speculated to be the reasons for the development of drug resistance and tumor progression [35,36]. Studies have found that MAPK pathway is highly expressed in drug-resistant tumors, so scholars speculated that BRAF inhibitor and MEK inhibitor could be used in combination to reduce drug resistance [37]. In order to prove the viewpoints above, a number of clinical studies have been carried out, which indeed received good clinical effects, showing its potential practicability [38].
In this study, using the datasets derived from the vemurafenib cell lines, we investigated the core genes Moreover, Enrichment results showed a signi cant upregulation in genes that are concerned with extracellular matrix and cellular adhesion. Whether vemurafenib is directly correlated with reduced cancer cell migration and invasion requires more validation. Despite all these ndings, we have seen enrichment in focal adhesion, which has been proved to participate in drug resistance by activating Bcl-2 through ERK signaling pathways. Moreover, we also found multiple tumorigenic genes downregulated that are can signi cantly reduce patient survival, including CDK1, CDK2, CDC45, PRC1, PLK1, CCNA2-most of them being key player in cell cycles , were at the core in the PPI network after vemurafenib application.
We found that the function of downregulated genes after vemurafenib treatment was enriched in cell replication (particularly in cell cycle), which was closely related to the expression of cyclin and CDK, revealing its signi cant role in inhibiting the growth of tumor cells. Tumorigenesis is the process of loss of normal regulation resulting in continuous growth. Cell cycle refers to the time that a cell goes through mitosis, including the G1, S, G2, and M phases [39]. The G1, S and G2 phases are collectively referred to as the intermitotic phase, which is the preparation period for the M phase [40]. G1/S checkpoint is an important site of cell regulation, which is completed by dephosphorylation and phosphorylation of cyclin, CDK and CKI [41]. Cyclin and CDK positively promote cell cycle, while CKI is involved in negative regulation. High expression of cyclin and CDK was detected in many tumors [42]. Cyclin D1 may be related to the activation of MAPK signaling pathway, and vemurafenib can inhibit this pathway. In our analysis, we also found that vemurafenib can reduce the expression of cyclin and CDK, suggesting the rationality of vemurafenib treatment. In conclusion, vemurafenib inhibits genes highly related to cyclin and CDK expression, thereby reducing tumor cell replication and slowing tumor growth.
Tumor microenvironment is one of the hotspots of current research, among which the extracellular matrix(ECM) is widely mentioned in the resistance of melanoma.
Fibronectin, a class of extracellular matrix, is often found to be expressed in drug-resistant cells with BRAF mutations .Its mechanism of resistance is activation of the α5β1/ Akt pathway [43]. In addition, Eishu Hirata et al. found that the abnormal activation of vemurafenib on melanoma-associated broblasts is associated with stromal generation and induction of integrin β1 (ITGB1)/FAK/Src signaling in melanoma cells, which provides a mechanism of drug resistance [44]. Our analysis found that the upregulated genes after vemurafenib administration are associated with ECM function, which may be a potential mechanism for drug resistance after targeted therapy.
There are many adhesion molecules on the surface of melanoma cells, such as cadherin, integrin and immunoglobulin superfamily [45]. The cell adhesion molecule MCAM/ MUC18 enables melanoma cells to have metastatic potential and enhanced tumorigenicity. MCAM/MUC18 mediated homomorphic and heteromorphic adhesion between melanoma cells and endothelial cells, respectively [46]. These two interactions may promote transfer at different stages of the cascade [47]. The loss of E-cadherin and the increase of N-cadherin in melanoma cells often disrupt the normal homeostasis, promote the occurrence and development of tumor, and may mediate the occurrence of epithelial-mesenchymal transformation (EMT) [48]. Therefore, some scholars speculate that the increase of intercellular adhesion may be the mechanism of drug resistance in targeted therapy. In our analysis, there was also an increase in cellular adhesion after vemurafenib treatment, which may be related to the development of drug resistance.
Neurogenesis often occurs after repairing nerve damage, such as a stroke. It is often accompanied by changes in the local microenvironment, causing the regeneration of blood vessels and nerves. Such local microenvironmental alterations are also present in melanoma cells. Roshini Prakash et al. found that neurogenesis could increase the brain metastasis of melanoma [49]. According to their study, neurogenesis can promote the action of melanoma on components such as ECM, potentially leading to the activation of some signaling pathways and increasing the metastatic potential of melanoma. In our study, vemurafenib led to the overexpression of genes whose functions were enriched in neurogenesis. This may also be related to the development of drug resistance after vemurafenib treatment, resulting in adverse outcomes.
Consistent with current clinical observations, vemurafenib has a dual effect in the treatment of melanoma. On the one hand, it targets BRAF mutations, causing downstream cascades (such as the MAPK pathway) to be suppressed, resulting in low expression of the associated proto-oncogenes. The results of our analysis showed that, the suppressed genes were involved in promoting cell replication, suggesting that vemurafenib inhibited tumor cell proliferation. However, on the other hand, vemurafenib can also cause some adverse events, such as the promotion of cell adhesion, changes in ECM and neurogenesis, which may be the result of the contradictory activation of vemurafenib and may be related to the acquisition of drug resistance, but it still needs to be con rmed by further studies.
Currently, the combination of BRAF inhibitor and MEK inhibitor has been found to reduce contradictory activation, thus avoiding resistance to vemurafenib, improving clinical outcomes and promoting patient survival. Our ndings can provide a direction for future exploration of the mechanism of vemurafenib drug resistance.But there are limitations to our study. First, we lack in vivo and in vitro trials to verify the results. And then, there was heterogeneity in the different experimental data we included.

Conclusion
Vemurafenib is one of the drugs commonly used in targeted therapies that inhibit BRAF mutations. In clinical trials, vemurafenib is found to have excellent clinical results at rst, but resistance emerges after a few months. Researches on drug resistance continue to mount. Our data analysis also found that vemurafenib had a dual effect: inhibiting the expression of harmful genes and promoting patient survival, but also promoting adverse events leading to drug resistance and tumor progression. Therefore, it is necessary to further study the drug resistance mechanism of vemurafenib and how to reduce its drug resistance, so as to shed light on the treatment of melanoma.

Materials And Methods
Oncogenic Properties and clinical Manifestations of BRAF UALCAN is an user-friendly interactive web-portal to perform to in-depth analyses of TCGA gene expression data.
We used the UALCAN(http://ualcan.path.uab.edu) databases to verify the differential mRNA expression of BRAF in human. We investigated relative BRAF mRNA expression in variant types of tumors and different stages SKCM patients using UACLAN databases.

TCGA datasets
The Cancer Genome Atlas (TCGA), which is a joint effort between NCI and the National Human Genome Research Institute beginning in 2006, molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types, bringing together researchers from diverse disciplines and multiple institutions. We investigated the TCGA datasets in order to gure out the genetic mutations in BRAF genes within skin cutaneous melanoma patients, and to analyze the survival curve between the mutated and not mutated cases.

GEO datasets
The Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo) is a public repository at the National Center of Biotechnology Information (NCBI) for storing high-throughput gene expression datasets. We selected potential GEO datasets using search terms "melanoma"[MeSH Terms] AND vemurafenib [All Fields] AND "Homo sapiens" [Organism] "Expression pro ling by high throughput sequencing" [ Dataset Type] in the GEO datasets to identify potential datasets. Then, we further screened these datasets according to the inclusion criteria above, deleting invalid outcomes and samples that are too few for statistical analysis (with sample volume less than 3), as well as single-cell sequencing analysis. Finally, datasets GSE42872 (cohort1), GSE127988 (cohort2), GSE110054 (cohort3) were included in our study, all of which were experiments conducted using melanoma cell lines M297.
R software and data processing R software (version4.0.3) was used to download and process the raw data. after the expression data has been processed, we used log2fold change to quantify the changes in expression between the control group and the vemurafenib group. R packages from Bioconductor( www.bioconductor.org ) were used in the analysis of the data; the volcano plots were drawn using the package "ggscatter"; the heatmaps were drawn using the packages "pheatmap". Pvalue<0.01 and logFC>2 was considered statistically signi cant.
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were conducted using the package "ClusterPro ler" to investigate whether those signi cantly upregulated or downregulated genes were enriched in certain pathways.

TIMER
Timer is an open web tool that can be used to explore the association of genes in immune cell in ltration in different tumors. With TIMER, B cells, CD4+T cells, CD8+T cells, macrophages, neutrophils, and dendritic cells can all be studied. In this paper, we used TIMER to explore the role of 10 highly expressed genes in immune in ltration in melanoma.

String database and Cytoscape software
The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, www.string-db.org) is a userfriendly, interactive database to analyze protein-protein interaction. In this study, we used the string database to investigate common differential genes in order to understand the protein interaction after using the Cytoscape, an open-source software platform for visualizing complex networks was used in our study to further analyze the genes that are analyzed using STRING database.
The genes were ranked according to their degrees (meaning the number of a protein's interaction with other proteins), and those genes with the highest degree scores (with scores ranking from 1 to 15) were recognized as hub genes.

Survival analysis
Gene Expression Pro ling Interactive Analysis(GEPIA, gepia.cancer-pku.cn ) is a interactive web tool for analyzing gene expression within different types of cancers and normal samples, based on the data acquired from TCGA database [11].
We conducted the survival analysis and drew the survival curve of between expression of the hubgenes from the PPI network, and the patient's survival outcome using the GEPIA web server. Therefore, the genes analyzed from GEO datasets are ultimately validated by TCGA data.

Declarations
Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and material All the data for this article can be found on TCGA, GEO, GEPIA.

Competing interests
All authors declare that No con ict of interest exists.