Abstract
Background/Aim: Approximately 50% of melanomas harbor the BRAF V600E mutation and targeted therapies using BRAF inhibitors improve patient outcomes. Nonetheless, resistance to BRAF inhibitors develops rapidly and remains a challenge in melanoma treatment. In this study, we attempted to isolate long noncoding RNAs (lncRNAs) involved in BRAF inhibitor resistance using a comprehensive screening method. Materials and Methods: We used a CRISPR-Cas9 synergistic activation mediator (SAM) protein complex in a genome-scale transcriptional activation assay to screen for candidate lncRNA genes related to BRAF inhibitor resistance. Correlation analysis was performed between expression levels of isolated lncRNA genes and IC50 of dabrafenib in a BRAF-mutated melanoma cell line. Next, online databases were used to construct the lncRNA–miRNA–mRNA regulatory network. Finally, we evaluated the significance of the expression levels of these lncRNAs and mRNAs as biomarkers using clinical specimens. Results: We isolated three BRAF inhibitor resistance–associated lncRNA genes, namely SNHG16, NDUFV2-AS1, and LINC01502. We constructed a lncRNA–miRNA–mRNA network of 13 nodes consisting of three lncRNAs, six miRNAs, and four mRNAs. The lncRNAs and target mRNAs from each regulatory axis significantly and positively correlated with each other. Finally, Kaplan–Meier analysis showed that higher expression levels of MITF, which was up-regulated by LINC01502, were significantly associated with worse prognosis in BRAF V600E-mutated melanoma. Conclusion: The identification of these BRAF inhibitor resistance–associated lncRNA genes at the genomic scale and the establishment of the lncRNA–miRNA–mRNA regulatory network provides new insights into the underlying mechanisms of BRAF inhibitor resistance in melanoma.
Approximately 50% of melanomas harbor the BRAF V600E mutation (1). Following the discovery of BRAF inhibitors, response rates, progression-free survival, and overall survival in patients diagnosed with BRAF-mutated metastatic melanoma have greatly improved (2, 3). Nonetheless, resistance to BRAF inhibitors develops rapidly and remains a challenge in treating melanoma (4, 5).
Drug resistance mainly results from genetic and epigenetic factors (6, 7). The main mechanism underlying the genetic factors in BRAF inhibitor resistance is the reactivation of the mitogen-activated protein kinase (MAPK) signaling pathway and activation of the phosphatidylinositol-3-kinase (PI3K)/AKT pathway (8). In recent years, several signaling pathways have been shown to be involved in the development of drug resistance in cancer cells. To date, pathways including the mammalian target of rapamycin (mTOR) pathway (9), Yes-associated protein 1/PDZ-binding motif (YAP/TAZ) pathway (10), c-Jun N-terminal kinase (JNK)/c-Jun pathway (11), and WNT5A/β-catenin pathway (12) highly correlated with BRAF inhibitor resistance.
In addition to genetic mechanisms, drug resistance can also emerge through epigenetic mechanisms (6). LncRNAs are transcripts longer than 200 nucleotides that are not translated into proteins. They are thought to be an essential factor involved in diverse cell functions and disease development through transcriptional interference, sense and antisense RNA hybridization, and interaction with proteins or small RNA precursors (13). Moreover, they can act as “miRNA sponges” that affect the competing endogenous RNA (ceRNA) network (14).
Recently, researchers have identified different lncRNAs that influence BRAF inhibitor resistance in melanoma (6). The overexpression of SAMMSON confers drug resistance in cancer cells (15), and down-regulation of TSLNC8 in BRAF inhibitor– sensitive cells reduce the cytotoxicity of the drug (16). However, few genome-scale screening studies have been conducted to identify potential regulators of BRAF inhibitor resistance. More drug resistance–associated lncRNAs are hypothesized to exist, that may play various roles in carcinogenesis and chemoresistance, but are yet to be discovered.
Herein, we screened lncRNA genes that may be related to BRAF inhibitor resistance, using CRISPR-Cas9 synergistic activation mediator (SAM) protein complex (17), and we then carried out genome-scale bioinformatics analyses to investigate their roles in carcinogenic processes and BRAF-inhibitor resistance regulation.
Materials and Methods
Cell lines and culture conditions. The human malignant melanoma cell line A375, harboring the BRAF V600E mutation, and HEK293 cell line were purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). A375 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Wako Bio, Tokyo, Japan) supplemented with 10% fetal bovine serum and 1% penicillin/streptomycin. HEK293 cells were cultured in DMEM supplemented with 10% fetal bovine serum. All cells were tested for mycoplasma contamination and were maintained at 37°C in 5% CO2.
Genome-scale transcriptional activation screen using CRISPR-Cas9 SAM system. To isolate lncRNAs related to BRAF inhibitor resistance via a genome-scale screening, we constructed a CRISPR-Cas9 SAM using a lentiviral transduction system and CRISPR lncRNA Activation Pooled Library (17) purchased from Addgene (#1000000106, Addgene, Watertown, MA, USA). The library provided three components: dCas9-VP64 fusion protein, MS2-P65-HSF1 plasmid, and single guide RNA (sgRNA). The sgRNA library consists of 96,458 independent sequences, which targeted 10,504 intergenic lncRNA transcriptional start sites (TSSs) for the transcriptional activation of lncRNA genes.
First, to produce lentiviruses, plasmids and packaging plasmids were co-transfected into HEK293T cells by using Polyethylenimine “Max” (PEI MAX) (MW 40,000), (#24765-100, Polysciences, Inc., Warington, PA, USA). After 48 h, the cell culture medium containing lentiviral particles was collected and filtered with a 0.45-μm Hawach Sterile PES Syringe Filter (#SLPES2545S, NIPPON Genetics, Tokyo, Japan). A Lenti-X™ Concentrator (#631231, TAKARA, Shiga, Japan) was used to determine the concentration of filtered infectious lentiviral vector particles.
Next, 5×105 cells of the A375 cell line were seeded into 6-well plates and cultured with the lentivirus-containing medium for two days to transduce the dCas9-VP64 and MS2-P65-HSF1 components at a multiplicity of infection (MOI) of less than 0.7. Stably expressing cells were selected using the antibiotic blasticidin (10 μg/ml) and hygromycin (1,000 μg/ml) (Wako Bio, Tokyo, Japan) for 7 days. Then, the sgRNA library was transduced into 1.8×108 established cells, and stable cell clones were selected using zeocin (800 μg/ml) (InvivoGen, San Diego, CA, USA) for 7 days.
Dabrafenib resistance screening. We divided the stably expressing cells into either the treatment or the control group. We first seeded 1.8×108 established cells into T-225 flasks for screening selection in different infection replicates. Culture medium with dabrafenib (97 nM) was added to the treatment group, while normal culture growth medium was added to the control group. Two infection replicates per screen were used to account for stochastic noise. Screening experiment was repeated two times.
gDNA sequencing and analysis. Total gDNA was extracted from 6×107 cells in different infection replicates using a Quick-DNA Midi Prep Plus Kit (ZYMO RESEARCH, Irvine, CA, USA), and PCR was performed on the sgRNA-containing region using EmeraldAmp PCR Master Mix (TAKARA) and analyzed using a 96-well T100 PCR thermal cycler (Bio-Rad, Hercules, CA, USA). The sequence of PCR primers is described in Supplemental Data 1.
Adaptors were added to each end of the PCR products via a second PCR for next-generation sequencing. The sgRNA sequences were analyzed on the next-generation sequencer NovaSeq 6000 platform (Illumina, San Diego, CA, USA). Adapter sequences from the DNA sequencing reads were trimmed using the program Cutadapt and regions with low-quality scores were removed using Trimmomatic (18). Comparisons between samples and between groups based on the output sgRNA count information were performed by using the CRISPR screening analysis program MAGeCK (19), in which the read counts of each sequence from BRAF inhibitor-selected cells were compared with matched read counts from the control cells.
Each target gene corresponded to multiple gRNAs of different sequences, and all sgRNAs targeting each gene were scored and ranked based on their aggregation values of negative/positive selection using a modified robust ranking aggregation (RRA) algorithm. RRA uses an underlying probabilistic model that makes the algorithm parameter free and robust to outliers, noise, and errors. Significance scores also provide a rigorous way to retain only the statistically relevant genes in the final list and the basis for reordering and identifying significant elements (20).
In addition, we obtained a spreadsheet of the lncRNA library with associated genomic locations in the hg19 reference genome from the Zhang laboratory (17). Conversion of TCONS/NR ID into Ensembl ID was carried out by confirming the location information in the Ensembl project GRCh37 (Supplemental Data 2, Gene ID converted).
Correlation analysis. In this study, the expression data of lncRNAs in cell lines was obtained from the TANRIC expression database (https://www.tanric.org). The IC50 values of the BRAF-mutated melanoma cell lines were obtained from the Cancer Dependency Map (DepMap) (https://depmap.org/portal/). Correlation analysis between the lncRNA expression levels and IC50 values of cell lines was carried out using EZR (ver 1.54; Saitama Medical Center, Jichi Medical University, Saitama, Japan), a graphical user interface for R (21).
Construction of lncRNA–miRNA–mRNA network. The mRNAs that exhibited a significantly positive correlation with the candidate lncRNAs were obtained from the TANRIC database. The target miRNAs of the mRNAs were predicted using miRWalk (22) and miRTarBase database (23). We selected common miRNA–mRNA relationships in both databases. DIANA-LncBase v3 is a reference repository with experimentally supported miRNA targets of noncoding transcripts, from which we obtained miRNA targets for the candidate lncRNAs (24). We then used the lncRNA–miRNA and mRNA–miRNA relationships with experimental evidence from these databases for further analysis. Visualization of the lncRNA– miRNA–mRNA network was performed using Cytoscape software (ver 3.91) (25).
Expression level of candidate genes and Kaplan–Meier survival analysis. We used the GEPIA database (http://gepia.cancer-pku.cn/) to analyze the expression levels of candidate genes in melanoma tumor tissues from The Cancer Genome Atlas (TCGA) database and in non-tumor tissues from the Genotype-Tissue Expression (GTEx) database.
We also conducted a Kaplan–Meier survival analysis on the candidate genes. TCGA data were obtained from cBioportal (https://www.cbioportal.org/). The subtypes and BRAF mutation information of the samples analyzed in this study were obtained from the downloaded clinical data.
Spearman correlation test and Kaplan–Meier survival of the lncRNA and mRNA expression levels between BRAF-mutated melanoma cohorts were performed using EZR (ver 1.54). The cutoff value for the expression level was determined at the median level. The p-values of the Kaplan–Meier survival curves were calculated using the log-rank test. The results were considered significant when the p-value was less than 0.05.
GO and KEGG pathway enrichment analysis. GO and KEGG pathway enrichment analyses were performed using David (ver 6.8) (26) and QIAGEN IPA software (QIAGEN Ingenuity Pathway Analysis). In a function annotation chart obtained from DAVID, GO and pathway terms with p-value less than 0.1 were considered a statistically significantly difference.
Results
BRAF inhibitor resistance–associated lncRNA genes identified in genome-scale CRISPR-Cas9 SAM screening. Study design and workflow are shown as a flowchart (Figure 1).
Design and workflow of the study.
To screen lncRNA genes related to BRAF inhibitor (dabrafenib) resistance, we used genome-scale CRISPR-Cas9 SAM system. In this study, an sgRNA showed a tendency to be targeting a BRAF inhibitor resistance–associated lncRNA when its count was larger in the treatment group than in the control group. From the next-generation sequencing results, all sgRNAs targeting each lncRNA were scored and ranked based on the modified robust-ranking RRA algorithm. All lncRNAs that ranked high in the screening results (indicating a high correlation with BRAF inhibitor resistance) were initially extracted, and then the top 50 lncRNAs with the greatest number of corresponding sgRNAs were selected from the initial list of highly ranked lncRNAs (Supplemental Data 2).
We analyzed the relationship between expression levels of the lncRNA genes and dabrafenib sensitivity in a melanoma cell line carrying the BRAF mutation. We found that the expression levels of SNHG16, NDUFV2-AS1, and LINC01502 positively correlated with the IC50 value of dabrafenib in the melanoma cell line (Figure 2), which suggests that these three lncRNAs are likely to be associated with BRAF inhibitor resistance.
Correlation between lncRNA expression levels and IC50 of dabrafenib in the BRAF-mutated melanoma cell lines. Correlation between the expression levels of candidate lncRNAs SNHG16 (A), NDUFV2-AS1 (B), LINC01502 (C), and dabrafenib IC50 values in the BRAF-mutated melanoma cell lines.
lncRNA–miRNA–mRNA network. lncRNAs affect the expression of other genes through various mechanisms. One of these mechanisms can be described by the ceRNA hypothesis. miRNAs regulate gene expression by interacting with their target genes. lncRNAs function as ceRNAs, competing for shared miRNAs to regulate the expression of mRNAs.
To further investigate the regulatory relationships of these candidate lncRNAs and predict potential mechanisms, we constructed a regulatory lncRNA–miRNA–mRNA network based on the ceRNA hypothesis using the interaction relationships of lncRNAs, miRNAs, and mRNAs obtained from databases (Supplemental Data 3).
The network included 13 nodes, consisting of three lncRNAs, six miRNAs, and four mRNAs (Figure 3). This network reflects the regulatory relationships in the expression of the candidate lncRNAs and provides a basis for predicting the functions and interactions of genes.
Interaction schema between miRNA, lncRNAs, and up-regulated mRNAs.
GO and pathway enrichment analysis. To study the potential canonical function of the three lncRNAs involved in the lncRNA–miRNA–mRNA network, we performed GO and KEGG analyses of mRNAs regulated by lncRNAs in melanoma using the DAVID database and IPA software.
The results highlighted that a number of mRNAs related to the candidate lncRNAs are involved in the occurrence and development of cancer, including cell division and focal adhesion. KEGG analysis revealed that the lncRNAs are involved in metabolic pathways and positively regulated the ERK1/ERK2 cascade and PI3K/AKT signaling pathway. The up-regulated mRNAs of SNHG16 and LINC01502 were mostly enriched in the metabolic pathways (Supplemental Data 4).
In addition, we conducted a pathway analysis on IPA to confirm whether these genes were involved in BRAF inhibitor resistance–related pathways, which have already been characterized and validated. The details of the BRAF-inhibitor resistance risk pathway enrichment analysis from the IPA database are shown in Supplemental Data 5. The results showed that copious genes were involved in BRAF inhibitor resistance–related pathways, including the ERK/MAPK signaling pathway, PI3K/AKT pathway, mTOR pathway, PTEN signaling, and FAK signaling.
Expression levels of specific genes in clinical samples. We also analyzed the expression levels of lncRNAs and mRNAs in melanoma tumor tissue samples (n=461) and non-tumor tissue samples (n=558). Compared to the levels in normal tissues, SNHG16, LINC0150, GPI, FAH, and MITF had higher expression levels in melanoma, whereas NDUFV2-AS1 and GNAL showed lower expression levels (Figure 4A-G).
Verification of expression levels and correlation of candidate genes in melanoma patients. (A-G) Gene expression levels of SNHG16 (A), NDUFV2-AS1 (B), LINC01502 (C), GPI (D), FAH (E), GNAL (F), and MITF (G) in TCGA_SKCM melanoma tumor and non-tumor tissue analyzed using GEPIA. (H-K) Correlation of candidate genes in melanoma patients from TCGA database.
To validate the regulation of mRNAs from the network by the candidate lncRNAs, we investigated the correlation between these mRNAs and their corresponding lncRNAs in a melanoma patient cohort from TCGA database. The correlation between lncRNAs and mRNAs indicated that all candidate lncRNAs, but NDUSFV2-AS1, were significantly up-regulated (Figure 4H-K), concomitant with the expression level of mRNAs participating in the same regulatory axis from the network in melanoma patients.
Kaplan–Meier survival analysis of candidate lncRNAs and mRNAs. Kaplan–Meier and log-rank tests were used to determine the relationship between specific gene expression levels and overall survival (OS) among the BRAF-mutated melanoma patient groups. Gene expression levels in each group were divided into high- and low-expression groups based on the median expression value.
In the BRAF-mutated melanoma group, high expression levels of MITF and FAH were significantly associated with worse prognosis, whereas GPI and lncRNAs exhibited no impact on prognosis (Figure 5).
Kaplan–Meier survival curves in melanoma patients with BRAF mutation, divided expression levels of candidate lncRNAs, and mRNAs. (A) SNHG16; (B) NDUFV2-AS1; (C) LINC01502; (D) GPI; (E) FAH; (F) GNAL; (G) MITF.
Discussion
In this study, three candidate lncRNAs, SNHG16, NDUFV2-AS1, and LINC01502, showed a positive correlation with dabrafenib IC50 in BRAF-mutated melanoma cell lines in our correlation analysis.
GO and KEGG pathway analysis revealed that mRNAs related to the candidate lncRNAs were involved in the occurrence and development of cancer, many of which are involved in BRAF inhibitor resistance–related pathways reported in past studies. Notably, the up-regulated genes of SNHG16 and LINC01502 were mostly enriched in metabolic pathways and were also found in cancer development pathways.
And in silico analysis according to the ceRNA hypothesis suggested that each of the three candidate lncRNAs may control the expression level of four mRNAs, GPI, FAH, GNAL and MITF.
According to previous studies, MITF, GPI, and FAH are involved in the lncRNA–miRNA–mRNA regulatory network, and all of them play an important role in the metabolic process. Considering that BRAF-targeted drugs can modulate the metabolic state of malignant melanoma cells (27), these genes in the network may affect the efficacy of drugs by regulating cancer metabolic processes.
SNHG16 has been widely reported to be enriched in colorectal cancer cell lines (28), bladder cancer (29) and other cancer cells (30) however, its relationship with BRAF has not yet been revealed. From the network, both GPI and FAH significantly positively correlated with SNHG16 in the melanoma patient group as target mRNAs of the lncRNA. GPI, a glucose-6-phosphate isomerase, plays an important role in cancer metabolism by acting as a glycolytic enzyme. Under hypoxic conditions, GPI expression can be increased, promoting glycolysis to provide ATP to cells. Glucose-6-phosphate isomerase deficiency results in mTOR activation (31) which leads to the emergence of BRAF inhibitor resistance (32). FAH encodes the last enzyme involved in tyrosine catabolism. This gene is activated by the ERK/MAPK pathway, promoting the TCA cycle, oxidative phosphorylation (OXPHOS), and production of ATP (33). It also activates the ERK pathway and induces mitotic abnormalities and genomic instability (34).
Mitochondrial microphthalmia-associated transcription factor (MITF), a target mRNA gene of LINC01502, is an important transcription factor that is involved in the regulation of mitochondrial function (35). It has been clearly shown to be related to the occurrence and development of melanoma, and MITF is a factor that contributes to BRAF inhibitor resistance (36). MITF is up-regulated following BRAF/MEK inhibitor treatment, and MITF overexpression confers resistance to both BRAF and MEK inhibitors in patients (37). In addition, a previous study indicated that MITF is phosphorylated in the mitochondria due to ERK1/2 activation. MITF then dissociates from pyruvate dehydrogenase (PDH), thereby increasing PDH functionality to increase mitochondrial OXPHOS activities (38). MITF directly regulates the expression of PGC1α. Activation of the BRAF/MAPK pathway can suppress MITF and PGC1α expression, thereby decreasing oxidative metabolism in melanoma. On the other hand, oxidative phosphorylation is promoted by MITF overexpression due to the use of BRAF inhibitor (39). Therefore, MITF is one of the most important factors in the development of BRAF inhibitor resistance.
Actually, the OS analysis showed that higher expression of MITF and FAH were significantly associated with a worse prognosis in the BRAF-mutated melanoma group in our study. But LINC01502 and SNHG16 expression levels were not significantly correlated with prognosis.
In fact, low LINC01502 expression was observed in many melanoma patients (Figure 4K); therefore, there was no statistically significant difference in prognosis when the cutoff value was set as the median. However, the expression of LINC01502 had a significantly positive correlation with MITF in melanoma patients at high expression levels, suggesting that a high expression of LINC01502 may induce a high expression of MITF in some melanoma patients. We hypothesized that LINC01502 can regulate MITF expression to confer resistance to BRAF inhibitors in melanoma by competing with shared miRNAs. Hence, we expect LINC01502 to be a potential biomarker that influences BRAF inhibitor resistance by regulating MITF expression. Moreover, additional therapies targeting lncRNAs such as oligonucleotide therapeutics may prevent the development of BRAF inhibitor resistance during treatment in melanoma patients.
Study limitations. First, in the genome-scale screening using CRISPR-Cas9, we did not successfully transduce all sgRNAs into the A375 cell line; hence, some lncRNAs that are worthy of further exploration may have been screened out at the outset. Second, although the regulatory network of the candidate lncRNAs has been constructed, the precise mechanism of some lncRNAs remains under investigation. The IPA results also show that the expression of the lncRNAs isolated in this study is correlated with the expression of genes involved in ERK/MAPK signaling. Previous reports have also shown that activation of the ERK/MAPK signaling pathway in melanoma is associated with malignancy and metastatic potential. EGF can influence BRAF inhibitor resistance through ERK/MAPK signaling pathway that makes proliferation of conjunctival malignant melanoma (40). Moreover, the overexpression of RCS2, which is negative regulator of MAPK and AKT signaling, can suppress melanoma growth (41). Additionally, the drug KS28 that inhibit RAF dimer formation correlated with overcoming BRAF inhibitor resistance in melanoma by down-regulating the MEK/ERK pathway (42). However, in this study, we could not find any possibility that the isolated lncRNAs and the molecules that make up the ERK/MAPK signal function as ceRNAs. It is possible that candidate lncRNAs play a role other than that of ceRNAs. Further validation experiments are required to determine the underlying mechanisms. Third, in the OS analysis, we were unable to determine the effect of gene expression on prognosis in the cohort of patients who received BRAF inhibitor treatment due to the lack of treatment information. Therefore, the significance of the expression of the identified lncRNA remains as a biomarker for predicting prognosis, rather than as a biomarker for predicting the effect of BRAF inhibitors.
In conclusion, we successfully identified three lncRNA genes (SNHG16, NDUFV2-AS1, and LINC01502) associated to BRAF inhibitor resistance using a genome-scale CRISPR screening technology. We also constructed a lncRNA-miRNA-mRNA network based on ceRNA hypothesis to predict the regulatory interaction mechanism. Bioinformatic analysis validation were carried out in both BRAF-mutated melanoma patients and cell lines. Most notably, 3 up regulated mRNAs (MITF, GPI and FAH) of candidate lncRNAs from our network are tightly associated with BRAF inhibitor resistance. In addition, MITF and FAH were significantly associated with worse prognosis in BRAF-mutated patients. The regulatory network will provide the novel idea to figure out the complete gene expression mechanism. The data presented here can provide a new strategy for further development of targeted drugs and help overcome drug resistance in melanoma patients with BRAF V600E mutation.
Acknowledgements
This work was supported by the JSPS KAKENHI Grant Number JP19K07740. The authors thank the Laboratory of Molecular and Biochemical Research, Biomedical Research Core Facilities, Juntendo University Graduate School of Medicine, Tokyo, Japan for technical assistance.
Footnotes
Authors’ Contributions
XW, MHa and RT conducted the experiments, analyzed and presented the data. MHo, MO, TF, SY, SK analyzed and validated the data. XW wrote original draft. SK was responsible for conceptualization, data curation, supervision, funding acquisition, and writing of the original draft. All Authors read and approved the submitted version of the manuscript.
Supplementary Material
Supplementary material can be obtained at:<https://www.dropbox. com/scl/fo/wtnuaci8islbhtaqfqth8/h?rlkey=4s1qtnm9cmr6j2g62jho2f 0z3&dl=0>
Conflicts of Interest
The Authors have no conflicts of interest directly relevant to the content of this article.
- Received March 21, 2024.
- Revision received April 6, 2024.
- Accepted April 17, 2024.
- Copyright © 2024 The Author(s). Published by the International Institute of Anticancer Research.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC-ND) 4.0 international license (https://creativecommons.org/licenses/by-nc-nd/4.0).
