An Integrative Gene Expression and Mathematical Flux Balance Analysis Identifies Targetable Redox Vulnerabilities in Melanoma Cells

Targeted therapy has been a major breakthrough for patients with melanoma harboring BRAF^V600 mutations, present in 50% of malignant melanomas (1). However, the clinical responses are highly variable and short lived, and relapse is almost universal (1). Overcoming reduced sensitivity and acquired resistance to targeted therapy is a major goal of current melanoma research. Several mechanisms of reduced sensitivity have been proposed (2), which have led to the development of several combination regimens in melanoma either with other targeted therapies (3) or in conjunction with immunotherapies (4). While these therapies improve responses, treatment outcomes still vary and benefits remain transient and unpredictable (4).

Some resistance can be attributed to genetic mutations (5), but accumulating evidence indicates that nongenetic processes play a critical role in response of cancer cells to drug treatment (6–8). Metabolic reprogramming, recognized as a hallmark of cancer (9), has recently emerged as a potential nongenetic process that contributes to the emergence of drug-tolerant cancer cells. Warburg first reported a link between tumor and metabolism in his influential observation that cancer cells convert most intracellular glucose to lactate, that is, aerobic glycolysis (10). Recent studies have shown that cancer cells can utilize glycolysis, mitochondrial respiration, or both, depending on their environment (11). Furthermore, cancer cells can adapt metabolically in response to external perturbations (12). This metabolic flexibility provides cancer cells with energy, and necessary intermediates for biosynthetic processes required for survival and to maintain redox balance under changing environments (12, 13). Furthermore, metabolic pathways are complex and interconnected, warranting a systems level approach to examine their relative importance, and how they change in cancer cells. Flux balance analysis (FBA) is the most commonly used mathematical modeling approach for genome-scale metabolic network reconstructions to estimate the role of metabolic reactions in a network (14). Such a quantitative approach can be utilized to predict global metabolic states of cancer cells under various conditions.

Melanoma cells upon BRAF inhibition have been shown to induce an enhanced oxidative phosphorylation (7, 15). Byproducts of augmented mitochondrial activity are reactive oxygen species (ROS), and this metabolic switch promotes oxidative stress in cells (16). However, the cellular response to ROS is rather complex: low levels facilitate intracellular signaling, while high levels may cause cell death (17). Therefore, cells require a robust antioxidant defense system to respond to an accumulation of ROS. By invoking antioxidant systems, cancer cells can clear excess ROS levels within a range that is not detrimental to them—a high redox homeostatic state (18). NADPH and glutathione (GSH) are two major antioxidants that maintain redox homeostasis in cells (19, 20). NADPH is produced via the pentose phosphate pathway, and other sources (21), and acts as a shared substrate for both GSH regeneration and ROS production, thus maintaining an optimum redox balance within cancer cells. However, optimum redox balance in cancer is context dependent, with tissue type and stroma greatly affecting how cancer cells survive, metastasize, and respond to therapies (20). Even within cancer subtypes, this redox state could be heterogeneous (22), motivating further research to examine how this balance could be altered for therapeutic benefits.

The link between oxidative stress and drug response in melanoma is beginning to be explored. BRAF inhibitor (BRAFi)-resistant melanomas were shown to upregulate NRF2-mediated antioxidant response to maintain cell survival (16). Nonetheless, it still remains to be examined how redox potential of melanoma cells affects their drug sensitivity, and how it is maintained under BRAF inhibition. These are important considerations, raising the possibility that redox balance can be modulated to enhance the effects of existing therapies.

Here, we show that antioxidant capacity of melanoma cells is linked to their drug sensitivity. Using an integrative approach through bioinformatics and FBA, we show that melanoma cells with reduced sensitivity to BRAFi exhibit an enhanced capacity for antioxidation and redox buffer, specifically through NADPH-oxidizing enzymes. By directly quantifying the redox-related metabolites, we confirm that drug-insensitive melanoma cells can maintain higher levels of antioxidant metabolites during treatment. Furthermore, we report that pharmacologic disruption of redox axis involving GSH enhances the effects of targeted therapies. We also extended our analysis to other melanoma subtypes in The Cancer Genome Atlas (TCGA), and found that elevated redox capacity could be a general feature of melanoma. Our results thus provide a proof of principle that redox vulnerabilities could be exploited for therapeutic benefits in melanoma.

All BRAF-mutated, NRAS-mutated, and NF1-mutated melanoma cells were engineered to express either histone 2B-mRFP using pHIV-H2B-mRFP plasmid, H2B-GFP, and/or fluorescent ubiquitination cell cycle indicator using Geminin1–110 monomeric Azami Green plasmid as described previously (23). Single cell–derived clonal populations of BRAF-mutated melanoma cells (SKMEL5) were selected by limiting dilution as described previously (8). Single cell–derived BRAF-mutated SKMEL5 subclones were derived as described previously (8). BRAF-mutated melanoma cells (SKMEL5, WM88), including the SKMEL5 subclones, NRAS-mutated melanoma cells (SKMEL2), and NF1-mutated melanoma cells (MeWO) were grown and cultured in DMEM and Ham F-12 media (DMEM:F12, 1:1, catalog no. 11330–032). Primary epidermal melanocytes, normal human adult (HEMa; PCS-200–013) were obtained directly from ATCC and cultured in Dermal Cell Basal Medium (ATCC-PCS-200–030) supplemented with Adult Melanocyte Growth Kit (ATCC-PCS-200–042). Media were obtained from Gibco and supplemented with 10% FBS. All cells were cultured in humidified incubators that were CO₂ and temperature (37°C) controlled. Cells were passaged 1–2 times per week and were maintained as exponentially growing cultures for a maximum of less than 20 passages. All cells were tested for Mycoplasma, and tested negative. Chemical reagents used are listed in Supplementary Information.

RNA-seq of melanoma cell lines (SKMEL5 subclones) was performed as reported previously (24). Transcript count data from Gene Expression Omnibus (accession number: GSE122041, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE122041) was used for downstream analysis performed in R (https://www.r-project.org). Differential expression analyses were performed using Bioconductor (https://www.bioconductor.org) packages DESeq2 and edgeR. Transcripts counts for BRAF-mutated melanoma cell lines were obtained from Cancer Cell Line Encyclopedia (CCLE) database. Refer to Supplementary Information for more details.

FBA is a computational approach for predicting genome-scale steady-state metabolic fluxes in samples of interest (14). Reaction fluxes throughout the metabolic network are predicted by solving the optimization problem:

where if m and r are the number of metabolites and reactions in the metabolic network, respectively, then v is a r × 1 vector of reaction fluxes to be solved for, S is the m × r stoichiometric matrix of the metabolic network, v_min,j and v_max,j are constraints imposed on each individual reaction flux, and f is an objective function that is maximized to maximize a particular metabolic phenotype of interest. Recon3D version 3.01 was used as the core metabolic network (25).

To maximize the oxidation of NADPH to NADP⁺ in the metabolic network, the following objective function was used in all cellular compartments:

This artificial sink on NADP⁺ in turn maximizes the flux through existing metabolic reactions that oxidize NADPH to NADP⁺. Individual reaction flux values were predicted by averaging the minimum and maximum predicted fluxes from flux variance analysis (FVA), which computes the minimum and maximum obtainable flux for each reaction under the model constraints and while maximizing the objective function above.

The dose–response curves are generated using a 2-fold dilution of single drugs at concentration indicated (highest) to zero (solvent the drug was dissolved in, either DMSO or H₂O), all concentrations contain equal percentage of the solvent used. For combination with PLX4720, 8 μmol/L PLX4720 was used as diluent, and serially diluted for the second drug from highest concentration to lowest. Direct measurements of cell counts were made using Cellavista software and ImageJ macros as described previously (8, 23, 26). Primary melanocytes were seeded into 96-well plates (∼3,000 cells/well) and drug treatments applied the following day, including DMSO or PBS control (all concentrations contained equal percentage of DMSO or PBS solvent). A 96-well plate was fixed at time 0 and was used as control for the treatment plates. Drug combination assay was performed as above and cells were treated for 24, 48, and 72 hours. At each time point, cells were fixed with 3% paraformaldehyde (PFA) for 20 minutes at room temperature, washed with PBS three times, and stored in a dark room until all plates were processed. When all plates were fixed, cells were stained with NucBlue Hoechst 33342 dye (Invitrogen, Ref# R37605, Lot# 2139329) at a concentration of 1 drop/mL of PBS as per the manufacturer's instructions. Cells were then imaged on the Evos M700 (Invitrogen) microscope using the DAPI and trans-luminescence channels for nuclear and whole cell imaging, respectively. The center 25% area of the well was imaged for all the plates. Hoechst-stained nuclei were segmented with a custom macro in FIJI (ImageJ), and cell nuclei counts for all treatments were then transformed on a log₂ scale and normalized to log₂ cell count at time 0.

For gene silencing experiments, ON-TARGETplus Human NOX5 siRNA SMARTPool (catalog no. L-010195–00) was used. ON-TARGETplus Non-Targeting pool (catalog no. D-001810–10–05) was used as scrambled control. NOX5 siRNA was used at the final concentration of 25 nmol/L. Both were obtained from Dharmacon and stored at −20°C until use. Transfection was carried out according to manufacturer's instructions using DharmaFECT1 Transfection Reagent (catalog no. T-2001–02). For combination with PLX4720, either scrambled mix or siRNA mix was used as diluent and the dose–response curves are generated using 2-fold dilution of PLX4720 at concentration indicated (highest, 32 μmol/L) to zero (either scrambled or siRNA alone), all concentrations contained equal percentage of the solvent used, and equal concentrations of either scrambled mix or siRNA mix. Melanoma cells were treated with siRNA or scrambled at day 0, and replaced with fresh reagents at day 3. Direct measurements of cell counts were made using Cellavista software and ImageJ macros as described previously (8, 23, 26). qRT-PCR quantification of NOX5 expression was done as described previously (24). For all drug responses including combination assay, primary melanocytes, and siRNA-drug combination, the drug-induced proliferation (DIP) rates are calculated using the slope of the log₂-normalized population curves after 24 hours.

CellRox DeepRed Reagent was used according to the manufacturer's instructions. Refer to Supplementary Information for more details.

For drug combinations indicated, the dose–response curves were generated as described above. For the indicated concentrations, two known antioxidants [reduced GSH (5 mmol/L) and ferrostatin-1 (FER, 5 μmol/L)] were added to the combinations at day 0, and replenished at 72 hours along with growth medium, inhibitors combination for the duration of the assay up to 6 days.

Global untargeted–targeted metabolomic analysis on BRAF-mutated melanoma cells (WM88, SKMEL5-SC10) treated with either DMSO or 8 μmol/L PLX4720 for 24 hours were done as described in Supplementary Information.

All data, calculated DIP rates, DEGs between subclones, SCGs for CCLE melanoma panel, expression data, codes used to generate the figures are available in the GitHub repo (https://github.com/paudelbb/Paudel_et_al_2019_Melanoma). Metabolomics data are available at the NIH Common Fund's National Metabolomics Data Repository website, the Metabolomics Workbench, https://www.metabolomicsworkbench.org where it has been assigned study ID (ST001339). RNA-seq of melanoma cell lines (SKMEL5 subclones) is available at the Gene Expression Omnibus (accession number: GSE122041).

We recently reported that BRAF-mutated melanoma cell lines, including isogenic single cell–derived subclones, exhibit varying drug sensitivities to a small-molecule BRAF kinase inhibitor (BRAFi; refs. 8, 26). Using DIP rates (Supplementary Table S1; ref. 23), as a measure of drug effects, we examined the molecular correlates of BRAFi sensitivity, and reported that top 200 differentially expressed genes (DEG) in drug-insensitive cells were enriched in processes and functions related to redox metabolism (24). Here, we extended our analysis to examine all significant DEGs among isogenic subclones (Supplementary Table S2) using their RNA-seq profiles. Among 2,165 DEGs, 1,361 (⅔) were upregulated, while 804 (⅓) were downregulated (Fig. 1A and B). Next, we projected all DEGs onto 10 BRAF-mutated melanoma cells from CCLE datasets (27), and calculated their correlation to DIP rates. We selected significantly correlated genes using a criterion as outlined in Fig. 1A, and identified both positively and negatively correlated genes (Fig. 1C; Supplementary Table S3). Gene Ontology (GO) analysis on positively correlated genes in both subclones (Fig. 1B) and CCLE melanoma panel (Fig. 1C) showed enrichment of molecular functions related to redox balance, coenzyme metabolic process, and oxidoreductase activity (Fig. 1D; Supplementary Table S4). Specifically, we observed enrichment of gene signatures [CYBA, NADPH oxidase 5 (NOX5), HTATIP2, SLC7A11, highlighted in magenta in Fig. 1B and C] related to redox balance and modulation, oxidoreductase activity, and reactions that utilize or consume antioxidant NADPH as a substrate (Fig. 1E). To quantify the contribution of redox enzymes toward NADPH oxidation, we next predicted steady-state fluxes through several major NADPH-oxidizing reactions using a previously published framework for integrating transcriptomic, kinetic, and thermodynamic data into FBA models of individual cancer cell lines (14, 21). The drug-insensitive SC10 model showed significantly greater total conversion of NADPH to NADP⁺ compared with the drug-sensitive SC01 model. While most reactions did not display significant differences in predicted flux distributions, the SC10 model had significantly higher fluxes through glutathione disulfide reductase (GSR), dihydrofolate reductase (DHFR), and NADPH oxidase (NOX1/3/5/CYBB; Fig. 1F; Supplementary Table S5A). Consistent with the results in isogenic subclones, the subset of CCLE melanoma cell lines that is insensitive to MAPK-pathway inhibition also had higher predicted fluxes through GSR, DUOX1/2 reactions (Fig. 1G; Supplementary Table S5B–S5C). Taken together, these results suggest that melanoma cells with reduced drug sensitivity, as measured by DIP rates, exhibit enhanced antioxidation and redox balance. Specifically, we found that NADPH oxidation reactions coupled to GSH exhibited higher gene expression and predicted fluxes in drug-insensitive melanoma cells.

To directly assess the metabolite abundance, we performed a combined global untargeted–targeted LC/MS-MS metabolomics analysis in BRAF-mutated melanoma cells under BRAF-inhibition. In particular, we selected one drug-sensitive (WM88) and one drug-insensitive (SC10) melanoma cell line, and compared their metabolite profiles both in the absence and presence of BRAFi (Supplementary Table S6 and S7). In total, 2,687 features or metabolites (unique retention time and m/z pairs) were detected, quantified, and analyzed in these studies. Unsupervised principal component analysis across all detected metabolites clustered sensitive melanoma cells away from insensitive cells at baseline along principal component 1 (PC1), while the drug-treated samples are separated along principal component 2 (PC2; Fig. 2A). In general, the cells exhibit distinct metabolite profiles, with BRAF inhibition inducing a marked metabolite profile shift (Fig. 2A and B). In these studies, 18 validated (Level 1) redox metabolites were quantified (Supplementary Table S8). These metabolites were validated through accurate mass measurements (<5 ppm error), similarities in isotope distribution (>90%), matched retention times, and matched fragmentation spectra based on in-house databases (28). In sum, we observed 2,687 metabolites and using the combined global metabolomic profiles were able to mine the data for high confidence metabolite identification for molecules of interest (redox, GSH synthesis, and energy metabolites).

To probe the treatment-dependent changes in metabolite abundances, we also performed pairwise comparisons between BRAFi-treated and control samples within each cell line. In particular, we focused on the redox-related metabolites based on our earlier results on antioxidation and redox balance. Differentially expressed metabolites were determined on the basis of a statistical cutoff of FDR-adjusted P-values < 0.001, and log-fold change greater or equal to 1. PLX4720 treatment significantly depleted the levels of energy metabolites such as GTP, UTP, ATP, and GDP in both cell lines (Fig. 2C; Supplementary Fig. S1A). Among the redox metabolites, there was an increase in monophosphate nucleotides such as UMP and AMP in SC10 (Fig. 2C). Interestingly, the GSH level was significantly reduced in WM88 but not in SC10 upon BRAFi treatment (Fig. 2C). GSH levels, however, were comparable between two cells at baseline (Fig. 2D; Supplementary Fig. S1B). Because GSH can exist in both reduced (GSH) and oxidized form (GSSG), we also quantified the levels of GSSG in both WM88 and SC10 cells. GSSG levels were lower in SC10 compared with WM88 at baseline, but similar in PLX4720 treatment (Fig. 2D). The ratio of reduced and oxidized GSH, often used as an indicator of oxidative stress in cells, was also calculated. These data show that SC10 cells had a significantly higher GSH:GSSG ratio at baseline compared with drug-sensitive WM88 cells, providing them a robust antioxidant capacity (Fig. 2D). Upon treatment, SC10 showed a slight decrease in GSH:GSSG ratio (μ_DMSO = 1.20, μ_PLX = 0.732, decrease of 1.6-fold). In contrast, WM88 showed an approximately 10-fold decrease in its ratio upon BRAF inhibition (μ_DMSO = 0.67, μ_PLX = 0.067; Fig. 2D).

In addition, we also quantified the abundance of cysteine and cystine—precursors of de novo GSH synthesis. Transport of cystine into cells to form cysteine depends on the activity of cystine/glutamate antiporter, SLC7A11, one of the genes positively correlated with DIP rates. Both WM88 and SC10 have relatively similar low abundance values for the cystine molecule. Because cystine is at/near the limit of detection and below the limit of quantification of this particular LC/MS-MS analytical method, we did not detect cystine; therefore, we quantified the ratio of cystine and cysteine in drug-treated conditions. The ratio of cystine and cysteine under treatment was significantly higher in WM88 than in SC10, indicating a rapid accumulation of extracellular cystine (Supplementary Fig. S1C). In addition, BRAFi treatment significantly increased the levels of cysteine in SC10 and not in WM88 (Fig. 2C). In addition, serine and glycine can be used for a substrate for GSH synthesis (29). Levels of serine were substantially higher in SC10 compared with WM88 at baseline, and upon BRAFi treatment were maintained in SC10, and significantly decreased in WM88 (Fig. 2C; Supplementary Fig. S1A and S1D). Furthermore, we examined the levels of cAMP at baseline and drug-treated conditions in these cells, and found that in both sensitive (WM88) and drug-insensitive (SC10), BRAFi treatment significantly decreased the levels of cAMP (Supplementary Fig. S1E). Interestingly, the levels of cAMP are much higher in SC10 compared with WM88 at baseline (Supplementary Fig. S1E). Among other energy metabolites, glycine increased in abundance in both cells upon BRAFi treatment, whereas glutamate decreased in both (Supplementary Fig. S1F and S1G). Collectively, these results indicate that drug-insensitive cells have a robust antioxidation and redox system that can be maintained under BRAF inhibition, perhaps through both an efficient regeneration and de novo synthesis of GSH.

On the basis of our results involving antioxidants NADPH (Fig. 1) and GSH (Fig. 2), we carefully constructed a schematic of reactions and identified the known inhibitors of the redox components (Fig. 3A). Because our data suggest a positive correlation between enhanced antioxidation balance in melanoma cells and reduced drug sensitivity, we examined whether we could alter sensitivity of melanoma cells to BRAFi by using redox inhibitors. To evaluate the role of different redox components, we tested BRAFi in combination with known redox inhibitors such as DPI, RSL3, buthionine sulfoximine (BSO), Erastin, and FK866. We subjected two cell lines, one sensitive (WM88) and one insensitive (SC10) to increasing concentrations of either DPI, RSL3, or BSO alone or in combination with 8 μmol/L PLX4720 (Fig. 3B). While both cells exhibited concentration-dependent antiproliferative effects to single redox inhibitors (DPI, RSL3, BSO), enhanced effects of combination were observed only in the insensitive, and not in cells already sensitive to PLX4720 (Fig. 3B). Similar effects were seen in A2058, considered as largely insensitive BRAF-mutated melanoma cells (Supplementary Fig. S2A). Depending on the drugs, the combination either increased efficacy (maximum effect), or potency (more effect with less drug) or both in drug-insensitive SC10. Because of the known off-target effects of the pharmacologic agents, we used siRNA-mediated knockdown of the identified targets to validate the effects of the combinations. We selected two drug-insensitive cell lines, A2058 and SC10, based on their NOX5 expression at baseline (Supplementary Fig. S2B). Consistent with the pharmacologic inhibition, siRNA knockdown of NOX5 enhanced the effects of BRAFi in NOX5^hi BRAF-mutated melanoma cells (Fig. 3C; Supplementary Fig. S2C). Interestingly, the combination increased both potency and efficacy in drug-insensitive cells (SC10, A2058) and was similar for most concentrations in drug-sensitive cells (WM88) (Fig. 3C; Supplementary Fig. S2C). However, not all redox inhibitors enhanced the effects of PLX4720. Erastin (SLC7A11 inhibitor) alone did not exhibit any effects on SC10, and its combination with PLX4720 was similar to the effects of PLX4720 alone (Fig. 3D, black dotted line). Similarly, the combination of FK866 (which inhibits the enzyme NAMPT and thus depletes NAD⁺) with PLX4720 was not advantageous compared with single agents alone. At lower concentrations, the effect of combination was similar to PLX4720 alone, and at higher concentrations, similar to FK866 alone, indicating that the combination did not enhance the effects of single agents (Fig. 3D). This is an interesting result, as it suggests that inhibition of only certain redox components may enhance the effects of BRAFi in BRAF-mutated melanoma cells.

Our drug combinations suggest that coadministration of certain redox inhibitors with PLX4720 enhance the benefits of targeted therapies in melanoma (Fig. 3). To investigate how coadministrations affect redox balance in melanoma cells, we stained BRAF-mutated melanoma cells (SKMEL5) treated with either single agents or combinations with CellROX Deep Red Reagent. Specifically, we tested three drug combinations that showed increased drug effects with BRAFi: DPI, RSL3, and BSO. Unlike in BSO and RSL3, the total cellular ROS levels were increased upon PLX4720 treatment (P < 0.001; Fig. 4A; Supplementary Fig. S3A), consistent with the results in earlier reports (19). Although counterintuitive, DPI treatment as a single agent induced a significant increase in ROS intensity compared with DMSO control (Fig. 4A; Supplementary Fig. S3A). This result is consistent with studies that suggest DPI is a nonspecific NOX inhibitor, is shown to inhibit xanthine oxidase, eNOS, and causes an enhanced oxidative stress in cells (30). Similarly, all three drug combinations (BSO+PLX, RSL3+PLX, and DPI+PLX) resulted in a significant increase in total cellular ROS compared with DMSO control, PLX4720 alone, or the respective single agents (Fig. 4A; Supplementary Fig. S3A). To examine whether ROS induction in drug combination differs by cell lines, we compared total ROS intensity in four different melanoma cells: two drug-insensitive cells (SC10, A2058) and two drug-sensitive cells (WM88, SC01). Our analyses suggest that the ROS induction due to combination varies by cell lines. In the concentrations tested, we saw a general trend toward an increase in ROS intensity in combinations in all cell lines. While only PLX4720+DPI induced significant ROS compared with DMSO in three cells (A2058, SC10, SC01), other treatment groups exhibited significant increase only in SC01 (Supplementary Fig. S3B). Because of substantial cell death in WM88, we did not have enough cells to examine total ROS. Our results reveal that coadministration of redox inhibitors with PLX4720 disrupts redox balance in cells, and causes an increase in oxidative stress that varies among cell lines based on their inherent antioxidant capacity. Therefore, we sought to test whether the lethal effects of these drug combinations could be mitigated by exogenous antioxidants. We excluded the combination of DPI and PLX4720, because of the broad off-target effects of DPI. For the other two combinations, we cotreated melanoma cells with two known antioxidants: reduced GSH and lipophilic antioxidant FER. The effects of PLX4720 alone could not be rescued by either GSH or FER (Supplementary Fig. S4), indicating that the antiproliferative effects of PLX4720 is not due to oxidative stress alone. Interestingly, supplementation of the culture medium with GSH substantially rescued both RSL3 and BSO-induced cell death in combination (Fig. 4B), while FER rescue was specific to RSL3 only, confirming the GSH levels as a mechanistic link for reduced drug sensitivity in BRAF-mutated melanoma cells.

While melanomas with reduced sensitivity to BRAF inhibition show enhanced antioxidative capacity potentiating redox-targeting therapies, it remains unclear whether other melanoma subtypes display similar or varying redox phenotypes. To investigate whether alterations in redox metabolism are common attributes across melanomas, FBA models of 103 melanoma tumors from TCGA and 860 normal skin tissues from the Genotype-Tissue Expression (GTEx) project were developed using their respective gene expression profiles (21, 31, 32). Melanoma tumor models displayed greater steady-state fluxes through major NADPH-oxidizing reactions including GSR and DHFR compared with normal tissue models (Fig. 5A). Interestingly, normal tissue models predicted significantly greater flux through DUOX1/2, an NADPH-oxidizing enzyme involved in cellular ROS generation; thus, increased activity of antioxidant-generating GSR and decreased activity of ROS-generating DUOX1/2 likely contribute to the enhanced antioxidative capacity in melanoma tumors compared with normal tissues. These findings are similar to the comparison in predicted metabolic fluxes between drug-insensitive SC10 cells and drug-sensitive SC01 cells (Fig. 1F), suggesting a possible progression in redox capacity from normal skin tissues (low capacity) to drug-sensitive melanomas (medium capacity) to drug-insensitive melanomas (high capacity).

To measure redox fluxes under MAPKi treatment, we integrated MAPK-pathway inhibition in our FBA models. On the basis of recent studies (33), we assumed that BRAFi inhibits glycolysis in implementing our models. Parmenter and colleagues demonstrated that BRAFi inhibits glycolysis by suppressing expression of glucose transporters (GLUT1/3) as well as hexokinase (HK2). To simulate this effect in FBA models, we examined the effect of both eliminating cellular glucose uptake and hexokinase activity on the ability of melanoma TCGA samples to oxidize NADPH to NADP⁺. We observed a large stratification when comparing predicted BRAFi sensitivity (maximal NADPH oxidation before vs. after BRAFi implementation) between melanoma tumor models (Supplementary Fig. S5A). We then performed FVA to identify individual reaction fluxes that correlate with predicted BRAFi sensitivity. Fluxes through DHFR, NOX1/3/5/CYBB, and NOX4 were significantly greater in predicted BRAFi-resistant tumors (Supplementary Fig. S5B). In addition, using our FBA models of TCGA tumors, a simulated metabolome-wide gene knockout screen was performed, where for each of the 3,268 genes in the Recon3D metabolic reconstruction, the effect of gene knockout on total NADPH oxidation in each melanoma model was determined (Supplementary Fig. 5B). Many of the gene knockouts with greatest effects on NADPH oxidation across TCGA melanoma models included major genes involved in the reduction of NADP⁺ to NADPH (e.g., GLUD1, GLUD2, IDH2, MTHFD1), as well as genes utilizing reduced NADPH to increase cellular antioxidant stores or decrease levels of ROS (e.g., GSR, CBR1; ref. 21). To determine whether simulated gene knockdown results are biased toward genes proximal to NADPH reduction in the metabolic network, next we calculated the distance between each metabolic reaction to the closest NADPH-reducing reaction. Gene distance has a mean value of 1.31 with a maximum value of 3. Although all of the top 50 genes from the simulated knockdown screen have distance values of 0 or 1, we observed small distance values for all genes, suggesting a high degree of connectivity to NADPH metabolism for all genes across the human metabolic network (Supplementary Fig. S5C). Together, these analyses suggest that our results are not biased by network distance and it is relevant that genes encoding enzymes of Fig. 5A are featured in the top 50 genes shown in Fig. 5B.

Next, we sought to test whether our identified redox-targeting therapies work in other subtypes of melanoma. We applied the drug combinations (BSO+PLX and RSL3+PLX) in NRAS-mutated (SKMEL2), and NF1-mutated (MeWo) melanoma cells. The combinations increased efficacy of BRAFi in both SKMEL2 and MeWo cells, consistent with results above (Fig. 5C). Although PLX4720 has some off-target effects at moderate to high concentrations in non-BRAF–mutated melanoma cells (34), we observed these effects are minimal in terms of growth rates (DIP rates) in PLX4720 alone in both SKMEL2 and MeWo cells and only enhanced when it is combined with redox inhibitors (Fig. 5C). Furthermore, to investigate the effects of these agents on normal cells, we treated primary HEMa (PCS-200–013) with increasing concentrations of BSO, RSL3, and combined them with 8 μmol/L PLX4720. The effects of drugs, both single, and double agents are statistically insignificant compared with DMSO control, suggesting that the primary cells are largely unaffected by these agents (Supplementary Fig. S6). Collectively, these findings suggest that altered redox metabolism and enhanced antioxidant production are common features among melanomas, and that inhibition of redox enzymes may be a viable treatment strategy for melanomas in general.

Overall, we report here that higher antioxidant ability of cells correlates to lower drug sensitivity to targeted therapies in melanoma. Our analyses reveal an association between an elevated redox capacity, involving two reducing cofactors NADPH and GSH, and reduced drug sensitivity in melanoma. Furthermore, we show that drug-insensitive melanoma cells not only possess a robust redox system, but also maintain key antioxidant metabolites during treatment. Indeed, we show that the disruption of this redox buffer improves the outcomes of targeted therapies. In summary, our findings suggest that modulation of cancer antioxidant defense could be exploited to augment the benefits of existing therapies in melanoma (16). This strategy may extend beyond melanoma, as antioxidant pathways have been implicated in tumor progression and drug resistance (22, 35, 36).

How ROS affects cancer cells has been a contentious issue for years. ROS likely elicits a broad spectrum of cellular responses (13, 17). Although still debated, its effects on cancer cells can best be described as hormesis (17); at low levels, it can mediate sustained network signaling, leading to enhanced proliferation (17), while at high levels, it can induce cellular damage and cell death (37). The level of ROS that a cell can handle, thus, likely depends on how functional its antioxidant machinery is. In general, tumors can sustain a much higher level of ROS compared with normal tissues by invoking their antioxidant systems, thereby creating a high redox homeostatic state (18). Melanoma tumors could possibly exist in such a high redox homeostatic state, as we see an increased activity of antioxidant-generating GSR, and decreased activity of ROS-generating DUOX1/2 in TCGA melanoma tumors compared with normal skin. However, there exists a significant variability even within tumors, captured well in melanoma cell lines, suggesting redox capacity perhaps could be on a continuum: low capacity (normal skin), moderate capacity (drug-sensitive melanomas), and high capacity (drug-insensitive melanomas).

On the basis of modeling by FBA, the high redox capacity in drug-insensitive melanomas is predicted to stem from elevated fluxes through metabolic reactions that generate important antioxidant molecules including GSH. This computational approach toward predicting genome-scale metabolic phenotypes has many advantages. By integrating the most comprehensive model of human metabolism to date with transcriptomic data from cell lines or tumors of interest, quantitative predictions on the metabolic status of these samples can be obtained, which would not be available from gene expression analysis alone (21, 25). In addition, although FBA models can only make steady-state flux predictions, the interconnections between 13,000+ metabolic reactions in the human metabolic network often yield predictions with greatly improved accuracy and utility compared with smaller dynamic metabolic models that are missing these genome-scale interactions (14). Because of these advantages, FBA has led to many successful predictions in cancer metabolism, providing powerful insights into both cancer pathophysiology and treatment strategies (38–40). In addition, we used FBA to predict metabolic fluxes for reactions in central carbon metabolism [glycolysis and tricarboxylic acid (TCA) cycle] in TCGA melanoma tumors, and compared them with GTEx normal skin tissues (Supplementary Fig. S7A). Our analysis suggests that many reactions had a statistically significant difference between melanoma tumors and normal skin, with most glycolytic and some TCA cycle fluxes significantly increased in melanoma tumors (Supplementary Fig. S7A). Interestingly, similar analyses between BRAF-mutated TCGA melanoma tumors compared with TCGA non-BRAF–mutated melanomas show no difference in predicted metabolic fluxes through central carbon metabolism (Supplementary Fig. S7B). By coupling FBA, gene expression, and metabolomics screen, we show that the drug-insensitive melanomas have an increased dependence on antioxidant systems. This increased dependence could present a therapeutic opportunity in melanomas. Indeed, our own results and recent reports inform that redox modulation could be a viable strategy to improve the treatment outcomes (41).

GSH is the most abundant antioxidant in the cell. It can be synthesized either de novo or regenerated using NADPH as the substrate (42). De novo GSH biosynthesis requires cysteine, which is imported into the cells by cystine/glutamate antiporter, SLC7A11 (35, 36). Our data suggest that drug-insensitive melanoma cells have both upregulated SLC7A11 expression, and higher predicted fluxes through GSH regenerating GSR, suggesting an efficient maintenance of GSH. Intriguingly, drug-insensitive melanoma cells can also maintain higher levels of GSH, and maximally utilize cysteine during treatment. In addition, the level of serine is higher in drug-insensitive cells, and increased upon treatment. Because serine is a precursor of cysteine, it may help cells maintain necessary levels of GSH (29). Therefore, we next questioned whether the modulation of GSH level affects drug sensitivity in melanoma. Inhibition of SLC7A11, however, did not improve the effects of BRAFi in melanoma cells, indicating that melanoma cells could maintain GSH levels through NADPH regeneration. In contrast, GSH depletion by BSO enhanced the effects of combination. Thus, GSH modulation appears as an exciting therapeutic target in melanoma. Efficient regeneration of GSH requires sufficient NADPH, another important cofactor in cells. Our work focused on NADPH oxidation (consumption) and examined the contributions of several enzymes that convert NADPH to NADP⁺. Sixteen of the top 50 enzymes that had the most effect on the predicted fluxes through NADPH are known to be involved in NADPH production pathways. Examining the importance of NADPH production in melanoma drug sensitivity is a possible area of future investigation. We speculate that drug-insensitive melanoma cells may also possess a robust NADPH production pathway. It is worth noting that our current work refers to the early response of melanoma to BRAFi, and examines the relationship between inherent redox potential in a cell and its initial drug sensitivity. Our drug-insensitive cells are not drug resistant, instead they respond to MAPKi albeit at varying levels, such that the drug responses are in a continuum. When chronically made resistant to BRAFi, melanoma cells will go through certain adaptations, including the emergence of drug-resistant phenotypes. When melanoma cells adapt to long-term perturbations and oxidative stress, they exhibit a strong activation of NRF2-mediated oxidative response and regeneration of reduced GSH (16). Melanoma cells undergo metabolic modifications including the changes in their mitochondria, redox alterations, and metabolic plasticity to survive prolonged exposure to MAPKi (43) and ultimately develop resistance.

Consistent with earlier reports (44, 45), we observed that BRAFi alone also induced ROS in melanoma. However, BRAFi-induced effects on melanoma cells could not be rescued by known antioxidants, suggesting that effects of BRAFi on drug-insensitive melanoma cells are not solely due to ROS. Another possibility is that MAPKi-insensitive melanoma cells perhaps utilize their robust antioxidant system to tolerate the induced ROS. It is implicated that MAPKi inhibits glycolysis in melanoma cells (15, 33, 43) and leads to an increase in ROS production (44). We postulate that inhibition of glycolysis might make melanoma cells overreliant on antioxidants. We believe that elevated antioxidant usage may be a strategy for melanoma cells to thrive initially under MAPKi. However, this very strategy makes them vulnerable to redox inhibition when redox pathways are cotargeted with MAPKi, consistent with recent studies that show metabolic vulnerability in SLC7A11^hi cancer cells (36). Interestingly, this vulnerability could be a general feature of melanoma tumors regardless of the mutational status, as indicated by FBA analysis on TCGA samples, and enhanced effects of PLX4720 when combined with redox inhibitors in SKMEL2 (NRAS-mutated) and MeWo (NF1-mutated) melanoma cells. Interestingly, the effects of these drugs and their combinations are statistically insignificant in primary melanocytes. Therefore, an increase in oxidative stress due to the combinations identified here, possibly disrupts that redox balance. One such agent for combination, RSL3, is known to induce ferroptosis in cancer cells (46). Our results are, thus, consistent with recent reports that suggest ferroptosis-inducing drugs may enhance the current treatment options for melanoma (47). Ferroptosis has recently been exploited in different cancer types as an alternative way to reduce the fraction of drug-tolerant persister cells (48). Recently, redox metabolite cystine and its antiporter SLC7A11 have been linked to ferroptosis in pancreatic cancer (46). Yet, one important consideration is that ROS is not one, but a family of chemical species. Therefore, it remains to be determined whether different species of ROS have different effects on melanoma cells. Interestingly, inhibition of NOX5 by genetic knockdown (siRNA) enhanced the effects of PLX4720 in NOX5^hi cells, suggesting the NOX5-generated ROS might exhibit proproliferation signaling. This is consistent with recent studies that indicate NOX5-mediated ROS-induced downstream signaling may contribute to the growth of melanoma cells and other human tumors (49). It is tempting to speculate the role of NADPH oxidases in maintaining cellular bioenergetics as described recently (50). Although additional experiments are necessary, redox signaling activated in part due to altered phosphatase function has been implicated in NOX5-mediated cell proliferation (49). Further studies are required to ascertain the species of ROS that have detrimental effects on melanoma cells. Taken together, these results combined with recent reports support the hormetic nature of ROS-mediated signaling. Molecular basis for hormesis would be an exciting area of future investigation. One could imagine this could be cell line– or tumor specific, as the gene regulatory network has been linked with the ability of cells to undergo metabolic plasticity (13). Future work seems necessary to also investigate the compensatory antioxidant pathways, and how they help maintain redox homeostasis in melanoma cells.

In summary, our integrative approach demonstrates that enhanced redox balance capacity provides an early survival advantage in melanomas treated with MAPKi. This finding extends beyond BRAF-mutated melanomas (Fig. 5) and, we speculate, other cancer types in which an overactive MAPK pathways is the oncogenic driver. Our analyses indicate that this finding is actionable and identify new vulnerabilities in melanoma cells—coadministration of redox inhibitors would significantly increase the effect of MAPKi in insensitive tumors. In virtually all cases, we observed both enhanced potency and efficacy (24) in such drug combinations compared with single agents. Intriguingly, no increased effects were observed in sensitive cells. These results suggest strategies for personalized interventions. Thus, biomarkers for redox potential in melanomas could be used to stratify patients, to identify the ones that would benefit from a triple drug combination (MAPKi plus redox inhibitors). Naturally, additional studies in large datasets are necessary to conclusively firm up correlations between BRAFi responsiveness and gene expression or fluxes involving redox pathways. We envision such studies as eminently feasible, integrating data from transcriptomics, epigenomics, and epigenetics, and possibly involving machine learning algorithms to substantiate the goodness of predictive biomarker sets.

B.B. Paudel reports grants from NCI (predoctoral and postdoctoral salary from grants listed in acknowledgements) during the conduct of the study; in addition, B.B. Paudel has a patent relating to the subject matter of this manuscript pending. P.E. Stauffer reports grants from NIH (T32 Training grant) during the conduct of the study. J.A. McLean reports other compensation from Waters Corporation (our laboratory is a Waters Center of Innovation) and Agilent Technologies Inc. (our laboratory is an Agilent Thought Leaders Laboratory) and personal fees and other compensation from MOBILion Systems Inc. (service on the scientific advisory board of MOBILion Systems Inc. and our laboratory is a beta-testing site for their instrumentation) outside the submitted work. V. Quaranta reports grants from NCI (listed in acknowledgments) during the conduct of the study and nonfinancial support from Parthenon Therapeutics, Inc. (as equity holder and academic cofounder) outside the submitted work; in addition, V. Quaranta has a patent for U.S. Prov. Appl. No. 1 pending (as an inventor relating to the subject matter of this manuscript) and a patent for U.S. Prov. Appl. No. 2 pending (as an inventor relating to the subject matter of this manuscript). No potential conflicts of interest were declared by the other authors.

B.B. Paudel: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing-original draft, project administration, writing-review and editing. J.E. Lewis: Software, formal analysis, investigation, visualization, methodology, writing-original draft, writing-review and editing. K.N. Hardeman: Investigation, writing-review and editing. C.E. Hayford: Investigation, writing-review and editing. C.J. Robbins: Investigation, writing-review and editing. P.E. Stauffer: Investigation, writing-review and editing. S.G. Codreanu: Formal analysis, investigation, methodology, writing-review and editing. S.D. Sherrod: Formal analysis, investigation, methodology, writing-review and editing. J.A. McLean: Formal analysis, supervision, methodology, writing-review and editing. M.L. Kemp: Formal analysis, supervision, funding acquisition, methodology, project administration, writing-review and editing. V. Quaranta: Conceptualization, resources, formal analysis, supervision, funding acquisition, methodology, project administration, writing-review and editing.

We are grateful to Meenhard Herlyn (Wistar Institute), Kim Dahlman, and Ann Richmond (Vanderbilt University) for kindly providing BRAF-, NRAS-, and NF1-mutated melanoma cell lines; to Joshua P. Fessel for insightful discussions on NADPH oxidases; to Darren R. Tyson, Christian T. Meyer, David J. Wooten, and Leonard A. Harris for useful discussions; to Jing Hao and Kevin Contreras for support in reagent procurement, and experimental preparation. We would also like to thank Prof. Kevin A. Janes (University of Virginia) for critical review of this manuscript.

This work was supported by the US NIH grants U54 CA217450, U01 CA215845, R01 CA186193, and U01 CA174706 (to V. Quaranta); U01 CA215848 (to M.L. Kemp); Vanderbilt Institute for Clinical and Translational Research (VICTR) grants 16721 and 16721.1 (to B.B. Paudel). J.E. Lewis is supported by NIH/NCI F30 CA224968. C.E. Hayford is supported by NIH/NCI F31 CA221147. P.E. Stauffer is supported by 2T32DK007563-31. K.N. Hardeman is now at Center for Human Nutrition, The University of Texas Southwestern Medical Center, Dallas, TX. C.J. Robbins is now at Yale University School of Medicine MD-PhD Program, New Haven, CT.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.