Changes in the Transcriptome and Chromatin Landscape in BRAFi-Resistant Melanoma Cells

Metastatic and drug-resistant melanoma are leading causes of skin cancer–associated death. Mitogen-associated protein kinase (MAPK) pathway inhibitors and immunotherapies have provided substantial benefits to patients with melanoma. However, long-term therapeutic efficacy has been limited due to emergence of treatment resistance. Despite the identification of several molecular mechanisms underlying the development of resistant phenotypes, significant progress has still not been made toward the effective treatment of drug-resistant melanoma. Therefore, the identification of new targets and mechanisms driving drug resistance in melanoma represents an unmet medical need. In this study, we performed unbiased RNA-sequencing (RNA-seq) and assay for transposase-accessible chromatin with sequencing (ATAC-seq) to identify new targets and mechanisms that drive resistance to MAPK pathway inhibitors targeting BRAF and MAPK kinase (MEK) in BRAF-mutant melanoma cells. An integrative analysis of ATAC-seq combined with RNA-seq showed that global changes in chromatin accessibility affected the mRNA expression levels of several known and novel genes, which consequently modulated multiple oncogenic signaling pathways to promote resistance to MAPK pathway inhibitors in melanoma cells. Many of these genes were also associated with prognosis predictions in melanoma patients. This study resulted in the identification of new genes and signaling pathways that might be targeted to treat MEK or BRAF inhibitors resistant melanoma patients. The present study applied new and advanced approaches to identify unique changes in chromatin accessibility regions that modulate gene expression associated with pathways to promote the development of resistance to MAPK pathway inhibitors.

Introduction

Melanoma is the deadliest form of skin cancer, and aggressive or drug-resistant forms of melanoma can metastasize to various distal organs (1, 2). Genomic sequencing has identified oncogenic BRAF mutations in greater than 50% of melanoma tumors (3, 4). Acquired oncogenic BRAF mutations result in the constitutive activation of BRAF→MEK→ERK (MAPK) pathway, which is necessary for melanoma growth and progression (5–7). These findings have led to the development of several BRAF and MEK inhibitors (BRAFi and MEKi, respectively) that have received approval by the US Food and Drug Administration (FDA) for use in the treatment of unresectable metastatic melanoma (8–10). However, despite an initial robust response to BRAF- and MEK-targeted therapies, patients with melanoma typically acquire treatment resistance within a few months, resulting in disease progression (11, 12). Due to the high prevalence of acquired resistance to BRAFi and MEKi, intensive efforts have focused on identifying the underlying mechanisms (4, 12, 13), contributing to substantial progress in the treatment of advanced and drug-resistant melanoma (14, 15). However, in a subset of cases, the mechanisms underlying acquired BRAFi and MEKi resistance remain unknown, and continued efforts toward identifying the drivers of treatment resistance in these cases remain necessary to identify more efficient, durable, and potentially personalized treatment options.

Mechanism of resistance to BRAF/MEK inhibitors can be MAP kinase pathway dependent and independent. The goal of our study was to investigate the role of MAP kinase pathway independent- reprogramming of chromatin landscape in MAPK pathway inhibitor resistant melanoma cells. Here, we focused on identifying alterations in genomic distribution of accessible chromatin site and its impact on transcriptional network to identify the functional genes and pathways that can be targeted for the treatment of MAPK pathway inhibitor resistant melanoma. To do so, in this study we performed deep-sequencing analyses in BRAFi resistant and sensitive melanoma cells, after confirming their growth phenotype. Our results revealed that the combined use of MAPK pathway inhibitors significantly inhibited the growth of BRAFi-sensitive cells but had only modest effects on BRAFi-resistant cells. RNA-sequencing (RNA-seq) analysis identified altered expression patterns for several genes involved in various functional pathways in BRAFi-resistant cells compared with BRAFi-sensitive cells. Assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) was performed to measure alterations in chromatin accessibility, which revealed that chromatin accessibility was altered in both the promotor, in between and downstream regions of several genes. These changes resulted in alterations in the expression of several functional pathways, promoting resistance. An integrated analysis combining the RNA-seq and ATAC-seq datasets identified a group of genes for which changes in chromatin accessibility aligned with the changes in mRNA expression levels. Pathway analysis of these differentially expressed genes revealed the involvement of several oncogenic signaling pathways in the development of resistance to MAPK pathway inhibitors in melanoma cells. The present study describes the application of new and advanced approaches for identifying changes in chromatin accessibility that modulate gene expression and signaling pathway activities to promote MAPK pathway inhibitor resistance in BRAF-mutant melanoma cells.

Materials and Methods

Cell Culture

M229 parent and M229 BRAFi-resistant lines were gifts from Roger Lo, University of California, Los Angeles, and Neil Rosen, Memorial Sloan Kettering Cancer Center, New York. They were maintained in a humidified atmosphere of 5% CO₂ at 37°C in Dulbecco’s modified Eagle medium (Life Technologies, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Life Technologies) and 1% penicillin/streptomycin (Life Technologies).

Cell Viability Assay (MTT Assay)

For MTT assays, 3 × 10³ melanoma cells (M229-Par and M229-Res) were plated in triplicate in a volume of 100 µL in 96-well plates. After 24 h, the cells were treated with different concentrations of inhibitors. Cell viability was evaluated 5 days after treatment. To measure cell viability, 20 µl 5 mg/mL MTT solution dissolved in 1× phosphate-buffered saline was added to each well of the 96-well plate and incubated for 2 h at 37°C. The MTT solution was then removed, and 100 µL dimethyl sulfoxide (DMSO) was added to each well. After the contents were mixed by pipetting, absorbance was measured at 590 nm and 630 nm using the Biotek Synergy MX Multi Format Microplate Reader (Biotek, Winooski, VT, USA). The average absorbance at 630 nm was subtracted from the average absorbance at 590 nm, and the growth rate was plotted relative to the growth rate of vehicle-treated cells.

Soft-Agar Assay

Soft-agar assays were performed by seeding 1 × 10⁴ melanoma cells (M229-Par and M229-Res) onto 0.4% low melting–point agarose (Sigma-Aldrich, Burlington, MA, USA) layered on top of 0.8% agarose. After 24 h, the cells were treated with different concentrations of inhibitors, as shown in the Figure 1, or DMSO (control). After 3–4 weeks of incubation, colonies were stained with 0.005% crystal violet and imaged under a microscope. Colony sizes were measured using ImageJ software (https://imagej.nih.gov/ij/) and plotted as percent relative colony size compared with control colonies. Statistical analysis was performed by Student’s t-test in GraphPad Prism 8.0 (GraphPad, San Diego, CA, USA).

FIGURE 1

Figure 1 Effects of MAPK pathway inhibitors on M229-Par and M229-Res melanoma cell growth. (A–C) The indicated melanoma cell lines were treated with the indicated concentrations of dabrafenib, vemurafenib, or trametinib for 5 days and subjected to MTT assays. Cell viability is plotted relative to DMSO-treated cells. (D) The indicated melanoma cell lines were treated with either 10 nM dabrafenib or 1 nM trametinib, alone or in combination, for 5 days and subjected to MTT assays. Cell viability is plotted relative to DMSO-treated cells. (E) The indicated melanoma cell lines were treated with either 0.5 µM vemurafenib or 1 nM trametinib, alone or in combination, for 5 days and subjected to MTT assays. Cell viability is plotted relative to DMSO-treated cells. (F) The indicated melanoma cell lines were treated with either 10 nM dabrafenib or 1 nM trametinib, alone or in combination, and subjected to clonogenic assays. Representative images are shown. (G) The indicated melanoma cell lines were treated with either 0.5 µM vemurafenib or 1 nM trametinib, alone or in combination, and subjected to clonogenic assays. Representative images are shown. (H) The indicated melanoma cell lines were treated with either 10 nM dabrafenib or 1 nM trametinib, alone or in combination, and subjected to soft-agar assays. Representative images are shown. (I) The indicated melanoma cell lines were treated with either 0.5 µM vemurafenib or 1 nM trametinib, alone or in combination, and subjected to soft-agar assays. Representative images are shown. (J–K) Relative colony sizes in the images shown in panels (H, I), respectively. Data represent the mean ± standard error of three biological replicates. ns = not significant, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Clonogenic Assay

The clonogenic abilities of vehicle-treated (control) and inhibitor-treated melanoma cells (M229-Par and M229-Res) were measured in triplicate using clonogenic assays in which 0.5 × 10³ cells were seeded in a 6-well plate. After 24 h, the cells were treated with different concentrations of inhibitors, as shown in the Figure 1, or DMSO (control). After 1–2 weeks of treatment, colonies were fixed with a fixing solution containing 50% methanol and 10% acetic acid, followed by staining with 0.05% Coomassie blue (Sigma-Aldrich). The relative number of colonies was calculated by counting the number of colonies in each sample and plotting the average number of colonies.

RNA-Sequencing

Total RNA was extracted from frozen cell pellet samples using the Qiagen RNeasy Plus Universal mini kit, according to the manufacturer’s instructions (Qiagen, Hilden, Germany). RNA samples were quantified using a Qubit 2.0 Fluorometer (Life Technologies), and RNA integrity was verified using an Agilent TapeStation 4200 (Agilent Technologies, Palo Alto, CA, USA). RNA-seq libraries were prepared using the NEBNext Ultra RNA Library Prep Kit for Illumina, according to the manufacturer’s instructions (NEB, Ipswich, MA, USA). Briefly, mRNAs were initially enriched with Oligo d(T) beads. Enriched mRNAs were fragmented for 15 minutes at 94°C. First-strand and second-strand cDNA were subsequently synthesized. cDNA fragments were end-repaired and adenylated at the 3’ends, and universal adapters were ligated to cDNA fragments, followed by index addition and library enrichment by PCR with limited cycles. The sequencing library was validated on the Agilent TapeStation and quantified using a Qubit 2.0 Fluorometer and by quantitative PCR (KAPA Biosystems, Wilmington, MA, USA).

The sequencing libraries were clustered on two lanes of a flowcell. After clustering, the flowcell was loaded on the Illumina HiSeq instrument (4000 or equivalent), according to the manufacturer’s instructions. The samples were sequenced using a 2 × 150-bp paired-end configuration. Image analysis and base calling were conducted using HiSeq Control Software. Raw sequence data (.bcl files) generated from Illumina HiSeq were converted into fastq files and de-multiplexed using Illumina bcl2fastq 2.17 software. One mismatch was allowed for index sequence identification.

RNA-Sequencing Analysis

After investigating the quality of the raw data, sequence reads were trimmed to remove possible adapter sequences and nucleotides with poor quality using Trimmomatic v.0.36. The trimmed reads were mapped to the reference genome available on ENSEMBL using the STAR aligner v.2.5.2b. The STAR aligner uses a splice aligner that detects splice junctions and incorporates them to help align entire read sequences. BAM files were generated during this step. Unique gene hit counts were calculated using feature counts from the Subread package v.1.5.2. Only unique reads that fell within exon regions were counted.

Differentially expressed genes were identified using the DESeq2 program (17). Genes showing altered expression with p < 0.05 and fold change > 1.5 were considered to be differentially expressed. Goseq (18) was used to perform the gene ontology (GO) enrichment analysis, and Kobas was used to perform the KEGG pathway analysis (19).

ATAC-Sequencing and Data Analysis

M229-Par and M229-Res cells were washed and treated with DNAse I (Life Tech, Cat. #EN0521) to remove genomic DNA contamination. Live cell samples were quantified and assessed for viability using a Countess Automated Cell Counter (ThermoFisher Scientific, Waltham, MA, USA). After cell lysis and cytosol removal, nuclei were treated with Tn5 enzyme (Illumina, Cat. #20034197) for 30 min at 37°C and purified with a MinElute PCR Purification Kit (Qiagen, Cat. #28004) to produce tagmented DNA samples. Tagmented DNA was barcoded with a Nextera Index Kit v2 (Illumina, Cat. #FC-131-2001) and amplified via PCR prior to an SPRI Bead cleanup to yield purified DNA libraries.

The reads were first mapped to the latest UCSC genome set using Bowtie2 version 2.1.0 (16). Mitochondrial reads, duplicate reads, and non-unique reads were removed before peak calling. MACS2 was used for peak calling using BAMPE mode (20). Differentially expressed peaks were identified using the DEseq2 program (17). Peaks showing altered expression with p < 0.05 and fold change > 1.5 were considered differentially expressed. Downstream genes of the differential peaks were used for GO and pathway enrichment analysis. Goseq (18) was used to perform the GO enrichment analysis, and Kobas was used to perform the KEGG pathway analysis (19).

Integrated Analysis of RNA-Seq and ATAC-Seq Data

RNA-seq and ATAC-seq data were analyzed to identify same-direction changes in mRNA expression and chromatin accessibility. This integration was used to assess pathway enrichment using KEGG pathway analysis (KEGG; www.genome.jp/kegg/pathway.html).

Statistical Analysis

All experiments were conducted in three biological replicates. The results for individual experiments were expressed as the mean ± SEM. P-values were calculated by t-test using GraphPad Prism version 8.0h for Macintosh, GraphPad Software, San Diego, California, USA (www.graphpad.com).

Results

Effects of MAPK Pathway Inhibitors on the Growth of BRAFi-Resistant and BRAFi-Sensitive Melanoma Cells

BRAFi-resistant M229 cells (M229-Res) were generated through the continuous exposure of BRAFi-sensitive M229 cells (M229-Par) to increasing concentrations of vemurafenib (PLX4032) in vitro, which allowed these cells to acquire a resistant phenotype (21). To determine whether the M229-Res cells generated by continuous vemurafenib exposure were also resistant to more potent and stronger BRAFi, such as dabrafenib, or to MEKi, such as trametinib, we performed short-term 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) survival assays. Vemurafenib is a potent inhibitor of B-Raf^V600E, with a half-maximal inhibitory concentration (IC₅₀) of 31 nM in a cell-free assay (22). Vemurafenib shows a 10-fold increase in selectivity for B-Raf^V600E compared with wild-type B-Raf in enzymatic assays, and cellular selectivity can exceed 100-fold for B-Raf^V600E compared with wild-type B-Raf (23). Dabrafenib is also a specific inhibitor of BRAFV600 mutants, with an IC₅₀ of 0.7 nM in cell-free assays, and presents with 7- and 9-fold less potency against wild-type B-Raf and c-Raf, respectively (24). Dabrafenib is approximately 50 times more effective than vemurafenib for inhibiting B-Raf^V600E (25). Trametinib is a highly specific and potent MEK1/2 inhibitor, with an IC₅₀ ranging from 0.92 to 1.8 nM in cell-free assays and displaying no inhibitory effects against the kinase activities of c-Raf, B-Raf, or ERK1/2 (26–28). Using the MTT assay, we found that the inhibition of the BRAF→MEK→ERK pathway by vemurafenib, dabrafenib, or trametinib resulted in stronger effects in M229-Par cells than in M229-Res cells (Figures 1A–C). We also observed that BRAFi and MEKi when used in combination, strongly inhibited the growth of M229-Par cells with modest effects on M229-Res cells (Figures 1D, E). Combined BRAFi and MEKi treatment also strongly inhibited the growth of M229-Par cells compared with M229-Res cells in clonogenic (Figures 1F, G) and soft-agar assays (Figures 1H–K). These results indicate that M229-Res cells generated by continuous exposure to vemurafenib also developed resistance against the more potent BRAFV600 inhibitor dabrafenib and the MEK1/2 inhibitor trametinib, both alone and in combination.

RNA-Sequencing Analysis Identifies Alterations in the Transcriptome of BRAFi-Resistant Relative to BRAFi-Sensitive Melanoma Cells

To identify transcriptional changes in BRAF-mutant M229 BRAFi-resistant cells as compared to BRAFi-sensitive cells, we performed RNA-seq comparing M229-Res cells with M229-Par cells. RNA-seq identified 12,314 differentially expressed genes (adjusted p-values ≤ 0.05), including 6,139 downregulated and 6,175 upregulated genes (Figure 2A; Supplementary Table S1) in M229-Res cells compared with M229-Par cells. A heat map showing the top 100 affected genes (50 upregulated and 50 downregulated) and a volcano plot showing the top 30 genes (15 upregulated and 15 downregulated), based on the log₂ fold change values, were plotted (Figures 2B, C).

FIGURE 2

Figure 2 RNA-sequencing identified differentially expressed mRNAs between M229-Par and M229-Res cells. (A) Total number of upregulated or downregulated genes with adjusted p-values ≤ 0.05 in M229-Res compared with M229-Par samples. (B) Heat map showing differentially expressed genes (up- or downregulated) in the indicated comparisons. The top 50 upregulated and the top 50 downregulated genes based on p-values are shown. (C) Volcano plot showing differentially expressed genes (up- or downregulated) in the indicated comparisons. The top 15 upregulated and the top 15 downregulated genes based on p-values are also labeled. (D, E) KEGG pathway analysis showing key upregulated (D) and downregulated (E) biological pathways associated with differentially expressed mRNAs in M229-Res cells compared with M229-Par cells.

To explore the functional pathways likely to be activated or repressed by changes in gene expression, we performed functional pathway analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG). Our results showed that genes upregulated in M229-Res cells compared with M229-Par cells were associated with the activation of many oncogenic pathways, including phosphoinositide 3-kinase (PI3K)–protein kinase B (AKT), MAPK, focal adhesion, and proteoglycan-based signaling cascades (Figure 2D). Additionally, the identified downregulated genes were associated with the inhibition of anti-proliferative pathways, such as cellular senescence, autophagy, and glycosaminoglycan degradation (Figure 2E). Our results demonstrate that transcriptional changes in M229-Res cells promote the activation of pro-oncogenic signaling pathways and the inhibition of anti-cancer pathways.

ATAC-Sequencing Analysis Identifies Alterations in Chromatin Accessibility Regions in BRAFi-Resistant Compared With BRAFi-Sensitive Melanoma Cells

Open chromatin accessible regions contain cis-regulatory elements that might modulate gene expression and activity (29, 30). We, therefore, performed chromatin accessibility profiling analyses in both M229-Par and M229-Res cells using ATAC-seq. We identified a total of 71,542 (adjusted p-values ≤ 0.05) accessible regions in M229-Res cells compared with M229-Par cells, including 37,814 regions with negative chromatin accessibility and 33,728 regions with positive chromatin accessibility (Figures 3A, B; Supplementary Table S2). Altered chromatin accessibility regions were distributed across upstream, intergenic, and downstream regions (Figure 3C). A heat map showing the top 100 genes (50 upregulated and 50 downregulated) due to changes in chromatin accessibility was plotted (Figure 3D). To understand how changes in chromatin accessibility impact the regulation of functional pathways involved in the development of resistance to MAPK pathway inhibitors in M229 cells, we performed functional pathways analysis using KEGG. Our results showed that changes in chromatin accessibility upregulated the expression of genes involved in several oncogenic pathways, such as the Rap1, Hippo, and extracellular matrix receptor interaction–dependent signaling cascades (Figure 3E). Apoptosis, among other pathways, was downregulated in M229-Res cells compared with M229-Par cells (Figure 3F). These results demonstrate that M229-Res cells display a unique chromatin accessibility profile that supports the development of MAPK pathway inhibitor resistance.

FIGURE 3

Figure 3 ATAC-sequencing identified differentially expressed mRNAs between M229-Par and M229-Res cells. (A) Heatmaps showing differential genomic regions with increased or decreased chromatin accessibility based on ATAC-seq in M229-Res compared with M229-Par samples. (B) Chromatin-immunoprecipitation (ChiP) peaks over chromosomes analyzed by ATAC-seq in M229-Res and M229-Par samples. (C) Pie-chart for the indicated samples mapping the locations of annotated peaks identified by ATAC-seq. (D) Heat map showing the top 50 upregulated and the top 50 downregulated genes with increased and decreased chromatin accessibility based on p-values. (E, F) KEGG pathway analysis showing key upregulated (E) and downregulated (F) biological pathways associated with genes located in regions with increased or decreased chromatin accessibility in M229-Res cells compared with M229-Par cells.

Integrated ATAC-Sequencing and RNA-Sequencing Data Analysis Identifies Alterations in Chromatin Accessibility Regions That Align With Changes in mRNA Expression Level in BRAFi-Resistant Compared With BRAFi-Sensitive Melanoma Cells

To identify correlations between accessible chromatin regions and mRNA expression levels, we integrated the data obtained from ATAC-seq and RNA-seq analyses, which revealed 5,646 significant changes in chromatin accessibility (adjusted p-values ≤ 0.05) aligned with significant changes in mRNA expression levels (adjusted p-values ≤ 0.05). Chromatin accessibility increased for 2,038 regions, associated with the upregulation of mRNA expression, whereas 3,608 chromatin regions became less accessible, inhibiting mRNA expression (Supplementary Table S3). A heat map was plotted to display the top 100 genes (50 upregulated and 50 downregulated) based on the log₂ fold changes in mRNA expression among those aligned with changes in chromatin accessibility (Figure 4A). These included transcriptional regulators, ion channels, enzymes, kinases, phosphatases, growth factors, G-protein coupled receptors and transmembrane receptors (Figure 4B). A subset of top upregulated and downregulated genes owing to their association with the pathways involved in promoting MAPK pathway inhibitor resistance and treatment failure, predicted 3-year survival outcomes among melanoma patients (Figure 4C).

FIGURE 4

Figure 4 Integrated analysis of ATAC-sequencing and RNA-sequencing to identify differentially expressed genes between M229-Par and M229-Res cells for which changes in mRNA expression aligned with changes in chromatin accessibility. (A) Heatmap for the top 100 (50 upregulated or with increased chromatin accessibility and 50 downregulated or with reduced chromatin accessibility) genes showing similar patterns in both the ATAC-seq and RNA-seq analyses in M229-Res cells compared with M229-Par cells. (B) Genes obtained integrated from analysis categorized based on function. (C) Survival analysis (3-year) for patients with melanoma according to high and low expression levels of genes showing similar patterns in both the ATAC-seq and RNA-seq analyses in M229-Res cells compared with M229-Par cells using the Human Protein Atlas dataset. (D, E) KEGG pathway analysis showing key upregulated (D) and downregulated (E) biological pathways associated with genes that display coherence between changes in chromatin accessibility and changes in mRNA expression levels in M229-Res cells compared with M229-Par cells. (F) Model: ATAC-seq integrated with RNA-seq is a new and advanced approach for identifying unique changes in chromatin accessibility regions that modulate gene expression and signaling pathway activities to promote the MAPK pathway inhibitor resistance.

To identify which functional pathways were affected by these gene expression changes, we performed KEGG pathway analysis (Supplementary Table S3) and identified several key biological processes enriched in M229-Res cells compared with M229-Par cells (Figures 4D, E), which may be involved in promoting a treatment-resistant phenotype. In particular, we observed the significant upregulation of PI3K–AKT, MAPK, Rap1, Ras, and proteoglycan, actin cytoskeleton and ECM receptor interaction dependent signaling pathways in M229-Res cells compared with M229-Par cells (Figure 4D). We also observed the inhibition of tumor growth inhibitory pathways, such as autophagy and the programmed death-ligand 1 (PD-L1) and programmed cell death protein 1 (PD-1) checkpoint pathway (Figure 4E). These results suggest that changes in the global chromatin state in M229-Res cells lead to altered gene expression associated with biological pathways that promote MAPK pathway inhibitor resistance in BRAF-mutant melanoma cells. Additionally, the expression of many genes predicted survival among melanoma patients.

Discussion

BRAF is a member of the Raf kinase family, and the oncogenic V600E mutation in BRAF has been identified in 90% of melanoma cases, leading to the activation of the MAPK pathway (31–33). Several oncogenic BRAF-targeting inhibitors have been approved by the US FDA, including vemurafenib and dabrafenib, for the clinical treatment of metastatic melanoma (8, 34, 35). Although BRAFi therapy results in an impressive initial clinical response against BRAF-mutant metastatic melanoma, the durability of this response is limited by the rapid emergence of acquired BRAFi resistance, which often occurs within a few months of treatment initiation (36–39). In the clinic, BRAFi therapy is often combined with other MAPK pathway inhibitors, such as MEKi, to obtain durable effects for the suppression of melanoma growth and the avoidance of drug resistance (6, 10, 12). However, acquired resistance to these agents remains a major hurdle preventing the success of targeted therapies and limiting their benefits. One approach to overcoming this limitation is to understand the mechanisms underlying acquired resistance (40–42), which can contribute to modifying therapeutic regimens or developing combination therapies to prevent the emergence of drug resistance.

In this study, we performed a large-scale, deep-sequencing analysis to investigate reprogramming of chromatin landscape, a MAP kinase pathway independent mechanism in acquired MAPK pathway inhibitor resistance in melanoma. Our results are summarized in Figure 4F. We found that melanoma cells that are resistant to the BRAFi vemurafenib are also resistant to the more potent BRAFi dabrafenib and to combination treatment including both BRAFi (vemurafenib and dabrafenib) and MEKi (trametinib). Previous studies have shown that cancer cells have distinct genetic, epigenetic and transcriptional states, which allows them to exist in heterogeneous functional populations and that poses a major obstacle to cancer diagnosis and treatment (43). To identify the differential chromatin state, MAPK pathway inhibitor–resistant and sensitive melanoma cells were analyzed using high-throughput sequencing methods. Analysis of the results let us to discover global changes in chromatin accessibility regions located upstream, in between, and downstream of numerous genes, resulting in changes in mRNA expression. These genes that were enriched in MAPK pathway inhibitor–resistant melanoma cells belonged to distinct functional classes which includes kinases, phosphatases, transcription regulator, transporters, growth factor, enzyme, g-protein coupled receptors. Previous studies have provided the evidence that kinases (44, 45), phosphatases (46), transcription regulator (47), transporters (48), growth factor (49, 50), enzymes (51), g-protein coupled receptors (52) drive drug resistance in cancer and provides the opportunity for targeting them to combat drug resistance. Additionally, the expression levels of many of these altered genes were able to predict survival in melanoma patients. Thus, in the future, these candidate genes may serve as biomarkers that can predict the subgroup of patients able to benefit from MAPK pathway inhibitor therapy, identify patients who will develop resistance to MAPK pathway inhibitors early on and select the optimal therapeutic approaches for treating these patients.

Pathway analysis performed with the genes identified as potentially upregulated due to changes in chromatin accessibility were associated with several known oncogenic signaling pathways involved in tumor growth, metastasis, and cancer drug resistance, such as Ras signaling (21, 53), the MAPK pathway (12, 54), and PI3–AKT signaling (55, 56). In addition to known tumor-promoting signaling cascades, our study identified pathways that never been investigated for their role in promoting resistance to MAPK pathway inhibitors in melanoma, such as relaxin, calcium, and Rap1 signaling. Additionally, our study identified many genes that were downregulated due to closed chromatin, associated with the suppression of anti-growth signaling pathways, such as senescence and apoptosis.

Apart from genes that changed at both chromatin accessibility level and transcription level, there were candidates that only changed at the mRNA level. This was because of the limitation of these two methods. Fundamentally, the transcriptome measured via RNA seq is the result of transcription, posttranscriptional regulation and RNA degradation, while the chromatin accessibility changes measured via ATAC seq provides information on chromatin accessibility across the genome that effect transcriptional initiation site availability at that particular time point. Additionally, the probability of ATAC-seq to not accurately predict all the changes in chromatin accessibility also exist. These are few limitations of using ATAC seq and RNA seq together to identify the genes that are altered in drug resistant state as compared to the sensitive state and has to be considered while concluding the results. Although important, to resolve this issue is beyond the current scope of our manuscript.

In sum, our results suggest that the development of drug resistance in cancer cells is a complex process. Hence, a deeper understanding of these newly identified mechanisms will provide better insights into the development of MAPK pathway inhibitor resistance in melanoma and potentially lead to more efficient treatment options. These studies also lay the foundation for further examinations of newly identified genes and pathways involved in the development of resistance to MAPK pathway inhibitors in our study in preclinical mouse models of melanoma. Collectively, our study results provide insight into the comprehensive changes in chromatin accessibility changes that regulate the transcriptional outputs and signaling cascades to promote resistance to MAPK pathway inhibitors in melanoma. Our study also identified new biomarkers, targets, and signaling pathways that can be investigated to formulate new melanoma treatments, particularly for patients who have developed resistance to MAPK pathway inhibitors (BRAFi + MEKi).

Data Availability Statement

The datasets presented in this study - the RNA-seq and ATAC-seq are deposited with accession number GSE202122 (GEO accession# GSE202118 for RNA-seq and GSE202117 for ATAC-seq) and can be found in online repository Z2VvQG5jYmkubmxtLm5vaC5nb3Y=.

Author Contributions

RG conceived the research and designed the experiments. KR and PG performed the experiments. YE helped with bioinformatics analysis. RG interpreted the data and wrote the manuscript. All authors reviewed, made suggestions, approved the paper, and provided comments.

Funding

We gratefully acknowledge grants from the National Institutes of Health (NIH) [R03CA221926 (to RG), R03CA230815 (to RG), R03CA248913 (to RG), and R01CA233481 (to RG)].

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.937831/full#supplementary-material

Supplementary Table 1 | All genes differentially expressed genes (adjusted p-values ≤ 0.05) in M229-Par and M229-Res cells identified via RNA sequencing.

Supplementary Table 2 | All significantly chromatin accessible regions altered with (adjusted p-values ≤ 0.05) in M229-Par and M229-Res cells identified via ATAC sequencing.

Supplementary Table 3 | Integration of ATAC sequencing and RNA sequencing data with significantly altered chromatin accessible regions aligning with gene expression changes in in M229-Par and M229-Res cells.

References

1. Eddy K, Shah R, Chen S. Decoding Melanoma Development And Progression: Identification Of Therapeutic Vulnerabilities. Front Oncol (2020) 10:626129. doi: 10.3389/fonc.2020.626129

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Davis LE, Shalin SC, Tackett AJ. Current State Of Melanoma Diagnosis And Treatment. Cancer Biol Ther (2019) 20:1366–79. doi: 10.1080/15384047.2019.1640032

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Scatena C, Murtas D, Tomei S. Cutaneous Melanoma Classification: The Importance Of High-Throughput Genomic Technologies. Front Oncol (2021) 11:635488. doi: 10.3389/fonc.2021.635488

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Janostiak R, Malvi P, Wajapeyee N. Anaplastic Lymphoma Kinase Confers Resistance To Braf Kinase Inhibitors In Melanoma. Iscience (2019) 16:453–67. doi: 10.1016/j.isci.2019.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Klein O, Clements A, Menzies AM, O'toole S, Kefford RF, Long GV. Braf Inhibitor Activity In V600r Metastatic Melanoma. Eur J Cancer (2013) 49:1073–9. doi: 10.1016/j.ejca.2012.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Luebker SA, Koepsell SA. Diverse Mechanisms Of Braf Inhibitor Resistance In Melanoma Identified In Clinical And Preclinical Studies. Front Oncol (2019) 9:268. doi: 10.3389/fonc.2019.00268

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu D, Liu Z, Condouris S, Xing M. Braf V600e Maintains Proliferation, Transformation, And Tumorigenicity Of Braf-Mutant Papillary Thyroid Cancer Cells. J Clin Endocrinol Metab (2007) 92:2264–71. doi: 10.1210/jc.2006-1613

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Subbiah V, Baik C, Kirkwood JM. Clinical Development Of Braf Plus Mek Inhibitor Combinations. Trends Cancer (2020) 6:797–810. doi: 10.1016/j.trecan.2020.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Tanda ET, Vanni I, Boutros A, Andreotti V, Bruno W, Ghiorzo P, et al. Current State Of Target Treatment In Braf Mutated Melanoma. Front Mol Biosci (2020) 7:154. doi: 10.3389/fmolb.2020.00154

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Gupta R, Bugide S, Wang B, Green MR, Johnson DB, Wajapeyee N. Loss Of Bop1 Confers Resistance To Braf Kinase Inhibitors In Melanoma By Activating Map Kinase Pathway. Proc Natl Acad Sci USA (2019) 116:4583–91. doi: 10.1073/pnas.1821889116

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Lelliott EJ, Mcarthur GA, Oliaro J, Sheppard KE. Immunomodulatory Effects Of Braf, Mek, And Cdk4/6 Inhibitors: Implications For Combining Targeted Therapy And Immune Checkpoint Blockade For The Treatment Of Melanoma. Front Immunol (2021) 12:661737. doi: 10.3389/fimmu.2021.661737

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Czarnecka AM, Bartnik E, Fiedorowicz M, Rutkowski P. Targeted Therapy In Melanoma And Mechanisms Of Resistance. Int J Mol Sci (2020) 21:4576. doi: 10.3390/ijms21134576

CrossRef Full Text | Google Scholar

13. Kakadia S, Yarlagadda N, Awad R, Kundranda M, Niu J, Naraev B, et al. Mechanisms Of Resistance To Braf And Mek Inhibitors And Clinical Update Of Us Food And Drug Administration-Approved Targeted Therapy In Advanced Melanoma. Onco Targets Ther (2018) 11:7095–107. doi: 10.2147/OTT.S182721

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Polkowska M, Czepielewska E, Kozlowska-Wojciechowska M. Drug Combinations As The New Standard For Melanoma Treatment. Curr Treat Options Oncol (2016) 17:61. doi: 10.1007/s11864-016-0436-y

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Abe R, Fujita Y, Yamagishi S. Angiogenesis And Metastasis Inhibitors For The Treatment Of Malignant Melanoma. Mini Rev Med Chem (2007) 7:649–61. doi: 10.2174/138955707780859440

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Langmead B, Salzberg SL. Fast Gapped-Read Alignment With Bowtie 2. Nat Methods (2012) 9:357–9. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Love MI, Huber W, Anders S. Moderated Estimation Of Fold Change And Dispersion For Rna-Seq Data With Deseq2. Genome Biol (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene Ontology Analysis For Rna-Seq: Accounting For Selection Bias. Genome Biol (2010) 11:R14. doi: 10.1186/gb-2010-11-2-r14

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, et al. Kobas 2.0: A Web Server For Annotation And Identification Of Enriched Pathways And Diseases. Nucleic Acids Res (2011) 39:W316–22. doi: 10.1093/nar/gkr483

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-Based Analysis Of Chip-Seq (Macs). Genome Biol (2008) 9:R137. doi: 10.1186/gb-2008-9-9-r137

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Nazarian R, Shi H, Wang Q, Kong X, Koya RC, Lee H, et al. Melanomas Acquire Resistance To B-Raf(V600e) Inhibition By Rtk Or N-Ras Upregulation. Nature (2010) 468:973–7. doi: 10.1038/nature09626

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Fisher R, Larkin J. Vemurafenib: A New Treatment For Braf-V600 Mutated Advanced Melanoma. Cancer Manag Res (2012) 4:243–52. doi: 10.2147/CMAR.S25284

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Bollag G, Hirth P, Tsai J, Zhang J, Ibrahim PN, Cho H, et al. Clinical Efficacy Of A Raf Inhibitor Needs Broad Target Blockade In Braf-Mutant Melanoma. Nature (2010) 467:596–9. doi: 10.1038/nature09454

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Menzies AM, Long GV, Murali R. Dabrafenib And Its Potential For The Treatment Of Metastatic Melanoma. Drug Des Devel Ther (2012) 6:391–405. doi: 10.2147/DDDT.S38998

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cebollero A, Puertolas T, Pajares I, Calera L, Anton A. Comparative Safety Of Braf And Mek Inhibitors (Vemurafenib, Dabrafenib And Trametinib) In First-Line Therapy For Braf-Mutated Metastatic Melanoma. Mol Clin Oncol (2016) 5:458–62. doi: 10.3892/mco.2016.978

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wu PK, Park JI. Mek1/2 Inhibitors: Molecular Activity And Resistance Mechanisms. Semin Oncol (2015) 42:849–62. doi: 10.1053/j.seminoncol.2015.09.023

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Thota R, Johnson DB, Sosman JA. Trametinib In The Treatment Of Melanoma. Expert Opin Biol Ther (2015) 15:735–47. doi: 10.1517/14712598.2015.1026323

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kostanyan T, Shazly T, Kaplowitz KB, Wang SZ, Kola S, Brown EN, et al. Longer-Term Baerveldt To Trabectome Glaucoma Surgery Comparison Using Propensity Score Matching. Graefes Arch Clin Exp Ophthalmol (2017) 255:2423–8. doi: 10.1007/s00417-017-3804-9

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang B, Lin J, Zhang J, Wang X, Deng X. Integrated Chromatin Accessibility And Transcriptome Landscapes Of 5-Fluorouracil-Resistant Colon Cancer Cells. Front Cell Dev Biol (2022) 10:838332. doi: 10.3389/fcell.2022.838332

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wang X, Yan J, Shen B, Wei G. Integrated Chromatin Accessibility And Transcriptome Landscapes Of Doxorubicin-Resistant Breast Cancer Cells. Front Cell Dev Biol (2021) 9:708066. doi: 10.3389/fcell.2021.708066

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Holderfield M, Deuker MM, Mccormick F, Mcmahon M. Targeting Raf Kinases For Cancer Therapy: Braf-Mutated Melanoma And Beyond. Nat Rev Cancer (2014) 14:455–67. doi: 10.1038/nrc3760

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ascierto PA, Kirkwood JM, Grob JJ, Simeone E, Grimaldi AM, Maio M, et al. The Role Of Braf V600 Mutation In Melanoma. J Transl Med (2012) 10:85. doi: 10.1186/1479-5876-10-85

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yaeger R, Corcoran RB. Targeting Alterations In The Raf-Mek Pathway. Cancer Discovery (2019) 9:329–41. doi: 10.1158/2159-8290.CD-18-1321

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhu Z, Liu W, Gotlieb V. The Rapidly Evolving Therapies For Advanced Melanoma–Towards Immunotherapy, Molecular Targeted Therapy, And Beyond. Crit Rev Oncol Hematol (2016) 99:91–9. doi: 10.1016/j.critrevonc.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Livingstone E, Zimmer L, Vaubel J, Schadendorf D. Braf, Mek And Kit Inhibitors For Melanoma: Adverse Events And Their Management. Chin Clin Oncol (2014) 3:29. doi: 10.3978/j.issn.2304-3865.2014.03.03

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Bugide S, Parajuli KR, Chava S, Pattanayak R, Manna DLD, Shrestha D, et al. Loss Of Hat1 Expression Confers Brafv600e Inhibitor Resistance To Melanoma Cells By Activating Mapk Signaling Via Igf1r. Oncogenesis (2020) 9:44. doi: 10.1038/s41389-020-0228-x

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Spagnolo F, Ghiorzo P, Queirolo P. Overcoming Resistance To Braf Inhibition In Braf-Mutated Metastatic Melanoma. Oncotarget (2014) 5:10206–21. doi: 10.18632/oncotarget.2602

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Savoia P, Zavattaro E, Cremona O. Clinical Implications Of Acquired Braf Inhibitors Resistance In Melanoma. Int J Mol Sci (2020) 21:9730. doi: 10.3390/ijms21249730

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Cohen JV, Sullivan RJ. Developments In The Space Of New Mapk Pathway Inhibitors For Braf-Mutant Melanoma. Clin Cancer Res (2019) 25:5735–42. doi: 10.1158/1078-0432.CCR-18-0836

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Guo W, Wang H, Li C. Signal Pathways Of Melanoma And Targeted Therapy. Signal Transduct Target Ther (2021) 6:424. doi: 10.1038/s41392-021-00827-6

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Kozar I, Margue C, Rothengatter S, Haan C, Kreis S. Many Ways To Resistance: How Melanoma Cells Evade Targeted Therapies. Biochim Biophys Acta Rev Cancer (2019) 1871:313–22. doi: 10.1016/j.bbcan.2019.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Villanueva J, Vultur A, Herlyn M. Resistance To Braf Inhibitors: Unraveling Mechanisms And Future Treatment Options. Cancer Res (2011) 71:7137–40. doi: 10.1158/0008-5472.CAN-11-1243

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yuan H, Yan M, Zhang G, Liu W, Deng C, Liao G, et al. Cancersea: A Cancer Single-Cell State Atlas. Nucleic Acids Res (2019) 47:D900–8. doi: 10.1093/nar/gky939

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Borgo C, Ruzzene M. Role Of Protein Kinase Ck2 In Antitumor Drug Resistance. J Exp Clin Cancer Res (2019) 38:287. doi: 10.1186/s13046-019-1292-y

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Petti C, Picco G, Martelli ML, Trisolini E, Bucci E, Perera T, et al. Truncated Raf Kinases Drive Resistance To Met Inhibition In Met-Addicted Cancer Cells. Oncotarget (2015) 6:221–33. doi: 10.18632/oncotarget.2771

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Stebbing J, Lit LC, Zhang H, Darrington RS, Melaiu O, Rudraraju B, et al. The Regulatory Roles Of Phosphatases In Cancer. Oncogene (2014) 33:939–53. doi: 10.1038/onc.2013.80

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Martinez-Balibrea E, Ciribilli Y. Editorial: Transcriptional Regulation As A Key Player In Cancer Cells Drug Resistance. Front Oncol (2021) 11:764506. doi: 10.3389/fonc.2021.764506

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Xiao H, Zheng Y, Ma L, Tian L, Sun Q. Clinically-Relevant Abc Transporter For Anti-Cancer Drug Resistance. Front Pharmacol (2021) 12:648407. doi: 10.3389/fphar.2021.648407

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Wilson TR, Fridlyand J, Yan Y, Penuel E, Burton L, Chan E, et al. Widespread Potential For Growth-Factor-Driven Resistance To Anticancer Kinase Inhibitors. Nature (2012) 487:505–9. doi: 10.1038/nature11249

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Dai Z, Huang Y, Sadee W. Growth Factor Signaling And Resistance To Cancer Chemotherapy. Curr Top Med Chem (2004) 4:1347–56. doi: 10.2174/1568026043387746

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Morad SAF, Cabot MC. The Onus Of Sphingolipid Enzymes In Cancer Drug Resistance. Adv Cancer Res (2018) 140:235–63. doi: 10.1016/bs.acr.2018.04.013

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lappano R, Maggiolini M. G Protein-Coupled Receptors: Novel Targets For Drug Discovery In Cancer. Nat Rev Drug Discovery (2011) 10:47–60. doi: 10.1038/nrd3320

CrossRef Full Text | Google Scholar

53. Yadav V, Zhang X, Liu J, Estrem S, Li S, Gong XQ, et al. Reactivation Of Mitogen-Activated Protein Kinase (Mapk) Pathway By Fgf Receptor 3 (Fgfr3)/Ras Mediates Resistance To Vemurafenib In Human B-Raf V600e Mutant Melanoma. J Biol Chem (2012) 287:28087–98. doi: 10.1074/jbc.M112.377218

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Wang L, Leite De Oliveira R, Huijberts S, Bosdriesz E, Pencheva N, Brunen D, et al. An Acquired Vulnerability Of Drug-Resistant Melanoma With Therapeutic Potential. Cell (2018) 173:1413–1425 E14. doi: 10.1016/j.cell.2018.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Chapman PB. Mechanisms Of Resistance To Raf Inhibition In Melanomas Harboring A Braf Mutation. Am Soc Clin Oncol Educ Book (2013) 33. doi: 10.1200/EdBook_AM.2013.33.80

CrossRef Full Text | Google Scholar

56. Haarberg HE, Smalley KS. Resistance To Raf Inhibition In Cancer. Drug Discovery Today Technol (2014) 11:27–32. doi: 10.1016/j.ddtec.2013.12.004

CrossRef Full Text | Google Scholar