Transcriptional, chromatin, and metabolic landscapes of LDHA inhibitor–resistant pancreatic ductal adenocarcinoma

Metabolic reprogramming, due in part to the overexpression of metabolic enzymes, is a key hallmark of cancer cells. Lactate dehydrogenase (LDHA), a metabolic enzyme that catalyzes the interconversion of lactate and pyruvate, is overexpressed in a wide variety of cancer types, including pancreatic ductal adenocarcinoma (PDAC). Furthermore, the genetic or pharmacological inhibition of LDHA suppresses cancer growth, demonstrating a cancer-promoting role for this enzyme. Therefore, several pharmacological LDHA inhibitors are being developed and tested as potential anti-cancer therapeutic agents. Because cancer cells are known to rapidly adapt and become resistant to anti-cancer therapies, in this study, we modeled the adaptation of cancer cells to LDHA inhibition. Using PDAC as a model system, we studied the molecular aspects of cells resistant to the competitive LDHA inhibitor sodium oxamate. We performed unbiased RNA-sequencing (RNA-seq), assay for transposase-accessible chromatin with sequencing (ATAC-seq), and metabolomics analyses of parental and oxamate-resistant PDAC cells treated with and without oxamate to identify the transcriptional, chromatin, and metabolic landscapes of these cells. We found that oxamate-resistant PDAC cells were significantly different from parental cells at the levels of mRNA expression, chromatin accessibility, and metabolites. Additionally, an integrative analysis combining the RNA-seq and ATAC-seq datasets identified a subset of differentially expressed mRNAs that directly correlated with changes in chromatin accessibility. Finally, functional analysis of differentially expressed metabolic genes in parental and oxamate-resistant PDAC cells treated with and without oxamate, together with an integrative analysis of RNA-seq and metabolomics data, revealed changes in metabolic enzymes that might explain the changes in metabolite levels observed in these cells. Collectively, these studies identify the transcriptional, chromatin, and metabolic landscapes of LDHA inhibitor resistance in PDAC cells. Future functional studies related to these changes remain necessary to reveal the direct roles played by these changes in the development of LDHA inhibitor resistance and uncover approaches for more effective use of LDHA inhibitors in cancer therapy.


Introduction
Metabolic reprogramming is a key hallmark of cancer cells and represents a therapeutic vulnerability in various cancer types (1)(2)(3). Notably, a number of enzymes and transporters that promote glycolysis and glucose metabolism are overexpressed in cancer cells relative to normal cells (4,5). One such enzyme, lactic acid dehydrogenase A (LDHA), is known to be overexpressed in several cancer types and associated with poor prognosis (6)(7)(8)(9)(10).
LDHA encodes the A subunit of lactate dehydrogenase, a metabolic enzyme that catalyzes the interconversion of lactate and pyruvate (11,12). The LDHA dependency of cancer cells has been demonstrated in several xenograft and genetically engineered mouse models, in which the genetic ablation of LDHA significantly reduces tumor growth (13)(14)(15). These studies prompted the development and testing of LDHA inhibitors (LDHAi) for cancer treatment (16)(17)(18)(19). One such inhibitor, sodium oxamate (hereafter referred to as oxamate), is a pyruvate analog and thus a competitive LDHAi that halts lactate production through LDHA inhibition (19). Therefore, similar to the genetic inhibition of LDHA, oxamate treatment also results in tumor suppression due to the loss of LDHA activity (19,20).
Cancer cells have been shown to depend upon various nutrients as either energy source or building blocks and these are shown to be necessary for their survival (2). In particular, depending upon whether cancer cells utilize primarily lipid or glucose they are categorized as lipogenic or glycolytic (21)(22)(23). A previous study used metabolic profiling to stratify pancreatic ductal adenocarcinoma (PDAC) into subtypes with distinct sensitivities to metabolic inhibitors (23). This prior study identified two major PDAC classes: those that are sensitive to lipogenesis inhibitors, such as stearoyl CoA desaturase (SCD) inhibitors, and those that are sensitive to glycolytic inhibitors, such as oxamate (23). LDHA is overexpressed in PDAC and predicts poor prognosis (24,25), and LDHA inhibition has been shown to suppress PDAC growth and progression (26,27).
In this study, we modeled the mechanism of resistance to the LDHAi oxamate, aiming to understand the mechanisms underlying cancer cell adaptation to glycolysis and LDHA inhibition, which may lead to improved glycolysis targeting and LDHAi-based cancer therapies. We used a glycolytic PDAC cell line MIAPaCa2 as a model system and employed several integrated and unbiased approaches, including RNAsequencing (RNA-seq)-based mRNA profiling, assay for transposase-accessible chromatin with sequencing (ATACseq)-based chromatin accessibility profiling, and large-scale untargeted global metabolomics analysis to identify the transcriptional, chromatin, and metabolic features of oxamatesensitive and oxamate-resistant PDAC cells and characterize their responses to oxamate treatment.

Results
Transcriptome-wide RNA-sequencing analysis identifies distinct mRNA expression profiles between oxamate response and oxamate resistance in PDAC cells To identify alterations that occur in pancreatic cancer cells able to survive LDHA inhibition, we treated the oxamatesensitive glycolytic PDAC cell line MIAPaCa2 (parental cells) with 10 mM oxamate for 4 weeks to generate an oxamateresistant PDAC cell line ( Figure 1A). MIAPaCa2 is a glycolytic cell line as observed by a previous study used metabolic profiling to stratify pancreatic ductal adenocarcinoma (PDAC) into subtypes with distinct sensitivities to metabolic inhibitors (23). Oxamate resistance in MIAPCa2 cells that survived oxamate treatment was confirmed using a clonogenic assay ( Figure 1B). We then performed a series of large-scale analyses, including RNA-seq, ATAC-seq, and untargeted global metabolomics, to uncover the differences between parental and oxamateresistant MIAPaca2 cells and their responses to oxamate treatment ( Figure 1C).
We first performed RNA-seq to monitor changes in mRNA expression in parental and oxamate-resistant MIAPaCa2 cells treated with either vehicle or oxamate. RNA-seq analysis showed that the treatment of parental MIAPaCa2 cells with oxamate led to significant changes in mRNA expression, evidenced by the identification of 1025 differentially expressed mRNAs following oxamate treatment compared with vehicle treatment (516  upregulated and 509 downregulated mRNAs; Figures 1D, E  and Table S1). However, oxamate-resistant MIAPaCa2 cells were largely transcriptionally non-responsive to oxamate treatment, with only 129 genes identified as differentially expressed following oxamate treatment compared with vehicle treatment (113 upregulated and 16 downregulated mRNAs;  Figures 1D, E and Table S2). Consistent with the finding that parental and oxamate-resistant MIAPaCa2 cells display significantly different mRNA expression profiles, we observed 2605 differentially expressed mRNAs (1897 upregulated and 708 downregulated mRNAs) in oxamate-resistant MIAPaCa2 cells  compared with parental MIAPaCa2 cells (Figures 1D, E and  Table S3). The differences were even more pronounced when both cell lines were treated with oxamate, resulting in 4040 differentially expressed mRNAs (2706 upregulated and 1334 downregulated mRNAs) in oxamate-resistant MIAPaCa2 cells compared with parental MIAPaCa2 cells (Figures 1D, E and  Table S4). Reactome Pathway Analysis of differentially regulated mRNAs across various comparisons revealed the enrichment of many distinctly regulated biological pathways ( Figures 1F-I).
Collectively, these results demonstrate that parental and oxamate-resistant PDAC cells show significantly different transcriptional profiles and significantly different responses to oxamate, resulting in the enrichment of distinct biological pathways.
ATAC-seq reveals significant differences in chromatin accessibility patterns between oxamate response and resistance in PDAC cells The upregulation or downregulation of mRNA expression can correlate with increased or reduced chromatin accessibility, respectively (28). Therefore, we measured changes in chromatin accessibility in parental and oxamate-resistant MIAPaCa2 cells treated with and without oxamate using ATAC-seq. The ATACseq results showed that parental MIAPaCa2 cells treated with oxamate underwent a 4-fold increase in chromatin accessibility changes (9906 differentially expressed peaks) compared with oxamate-resistant MIAPaCa2 cells treated with oxamate (2362 differentially expressed peaks) (Figures 2A, B). These results were consistent with the larger number of differentially expressed mRNAs observed in oxamate-treated parental MIAPaCa2 cells than in oxamate-treated oxamate-resistant MIAPaCa2 cells. Similarly, parental MIAPaCa2 cells treated with oxamate were significantly different from oxamateresistant MIAPaCa2 cells treated with oxamate and showed 32,999 changes in chromatin accessibility. A greater than 2fold increase in mRNAs were differentially expressed in the comparison between parental MIAPaCa2 cells and oxamateresistant MIAPaCa2 cells than in the comparison between vehicle-treated and oxamate-treated parental cells. However, a similar number of changes in chromatin accessibility was observed between vehicle-(8762 differential peaks) and oxamate-treated parental MIAPaCa2 cells (9906 differential peaks) ( Figure 2A). As expected, the Integrative Genomics Viewer (IGV) tracks were consistent with the expected changes in chromatin accessibility; two such examples, for the genomic regions encoding MARS1 and FAR2, under the indicated sample conditions are shown in Figure 2C. Under all conditions, a large majority of chromatin accessibility changes were associated with promoter regions, including approximately 35% in oxamate-treated parental MIAPaCa2 cells, approximately 38% in vehicle-treated parental MIAPaCa2 cells, approximately 40% in vehicle-treated oxamate-resistant MIAPaCa2 cells, and approximately 42% in oxamate-treated oxamate-resistant MIAPaCa2 cells ( Figure 2B). Collectively, these results demonstrate that similar to the differences in mRNA expression patterns, parental and oxamate-resistant M I A P a C a 2 c e l l s d e m o n s t r a t e d i s t i n c t c h r o m a t i n accessibility profiles.
Integration of RNA-seq and ATAC-seq data reveals direct correlations between differential mRNA expression and changes in chromatin accessibility To identify differentially expressed mRNAs associated with changes in chromatin accessibility, such as upregulated mRNAs in areas of increased chromatin accessibility and vice versa, we integrated the RNA-seq and ATAC-seq data across various comparisons. Of 1025 differentially expressed mRNAs identified in the comparison between vehicle-and oxamatetreated parental MIAPaCa2 cells, only 79 mRNAs (approximately 8%) were associated with similar directional changes for both mRNA expression and chromatin accessibility ( Figure 3A). Only 1 of 129 (less than 1%) differentially expressed mRNAs in the comparison between vehicle-and oxamate-treated oxamate-resistant cells displayed similar directional changes for both mRNA expression and chromatin accessibility ( Figure 3A). Furthermore, 310 of 2605 differentially expressed mRNAs in the comparison between parental and oxamate-resistant cells (approximately 12%) showed similar directional changes for both mRNA expression and chromatin accessibility ( Figure 3A). Of 4040 differentially expressed mRNAs identified in the comparison between oxamate-treated parental and oxamate-resistant MIAPaCa2 cells, 193 (approximately 5%) showed similar directional changes for both mRNA expression and chromatin accessibility. Collectively, the integrated analysis revealed a subset of mRNAs with expression changes that occurred concurrently and in the same direction as changes in chromatin accessibility.
Next, we performed Reactome Pathway Analysis for only those differentially expressed mRNAs associated with similar directional changes in chromatin accessibility. The integrated analysis of changes in mRNA expression and chromatin accessibility revealed no specific biological pathways enriched in oxamate-treated oxamate-resistant MIAPaCa2 cells compared with vehicle-treated oxamate-resistant MIAPaCa2 cells, indicating that these cells were largely non-responsive to oxamate treatment. However, various distinct biological pathways were enriched in the other comparisons, as shown in Figure 3B.
We then aimed to identify transcription factors potentially involved in the regulation of mRNA expression of genes associated with similar directional changes in chromatin accessibility. To this end, we used the Enrichr platform together with the Encyclopedia of DNA elements (ENCODE) and chromatin immunoprecipitation (ChIP) Enrichment Analysis (ChEA) to identify consensus transcription factors using the ChIP-X option (29). This particular setting of Enrichr allows for the identification of consensus target genes regulated by transcription factors found in the ENCODE and ChEA databases. For transcription factors with data from multiple experiments, set intersection was applied to obtain consensus. For each comparison, multiple transcription factors were identified as potentially direct transcriptional regulators for differentially expressed mRNAs ( Figure 3C and Tables S5-S7). In particular, in the comparison between parental and oxamate-resistant MIAPaCa2 cells, androgen receptor (AR), SRY-box transcription factor 2 (SOX2), transcription factor 3 (TCF3), SUZ12, and RELA were identified as the top-most significantly enriched transcription factors. In the comparison between vehicle-and oxamate-treated parental MIAPaCa2 cells, estrogen receptor 1 (ESR1), CCAAT enhancer-binding protein beta and delta (CEBPB and CEBPD, respectively), CCCTC-binding factor (CTCF), and RAD21 were identified as the top-most significantly enriched transcription factors. In the comparison between oxamate-treated parental and oxamate-resistant MIAPaCa2 cells, upstream binding transcription factor (UBTF), E74-like ETS transcription factor 1 (ELF1), breast cancer 1 (BRCA1), TATA-box-binding protein-associated factor 1 (TAF1), and PBX homeobox 3 (PBX3) were identified as the top-most enriched transcription factors.

Untargeted global metabolomics analysis identifies metabolic patterns in parental and oxamate-resistant PDAC cells
We next decided to examine metabolic changes in parental and oxamate-resistant MIAPaCa2 cells treated with and without oxamate. We first performed a Seahorse analysis to measure glycolytic functions in parental and oxamate-resistant MIAPaCa2 cells. We found that oxamate treatment resulted in reduced extracellular acidification in oxamate-resistant MIAPaCa2 cells compared with parental MIAPaCa2 cells, suggesting that oxamate-resistant cells showed inhibition of glycolysis ( Figure 4A). This finding is consistent with the RNA-seq and ATAC-seq data results ( Figure 1D Figure 3A), which showed that oxamate-resistant cells were largely non-responsive to oxamate treatment at the mRNA expression and chromatin accessibility levels, respectively.
We then performed untargeted metabolomics analyses in parental and oxamate-resistant MIAPaCa2 cells treated with and without oxamate using capillary electrophoresis time-of-flight mass spectrometry (CE-TOFMS), which detected a total of 228 different metabolites ( Figure 4B and Table S8). As expected, the metabolomics analysis identified significant differences in lactic acid levels, which were significantly reduced in oxamate-treated parental MIAPaCa2 cells compared with untreated parental cells. Oxamate-resistant MIAPaCa2 cells also showed significantly reduced lactic acid levels compared with parental cells ( Figure 4C). Furthermore, significant differences in the levels of other glycolytic pathway metabolites ( Figures 4D, E) and those associated with other metabolic pathways were also observed ( Figures 5-7). These results indicate that, similar to the differences observed in mRNA and chromatin accessibility profiles, metabolic profiles reveal major differences between parental and oxamate-resistant MIAPaCa2 cells and their responses to oxamate treatment.

Alterations in the mRNA expression levels of metabolic genes in parental and oxamate-resistant PDAC cells and their enrichment in various metabolic pathways
We next used our RNA-seq data to specifically examine the expression of metabolic genes and transporters using a combination of metabolic genes defined by the Kyoto Encyclopedia of Genes and Genomes (KEGG) and a list described in a previously published study that includes transporters able to impact metabolic pathways (Table S9) (30). The comparison between parental and oxamate-resistant MIAPaCa2 cells revealed the differential expression of 339 metabolic genes or transporters (Table S10 and Figure 8A). Reactome Pathway Analysis of these 339 metabolic genes or transporters revealed that these differentially expressed genes were significantly enriched in pathways associated with the metabolism of lipids, carbohydrates, nucleotides, and inositol phosphate and the reversible hydration of carbon dioxides, among others ( Figure 8B and Table S11). Similarly, the FIGURE 5 Metabolites generated during central carbon metabolism in parental and oxamate-resistant cells treated with or without sodium oxamate. Changes in the indicated metabolites generated during central carbon metabolism under the indicated conditions. comparison between vehicle-treated and oxamate-treated parental MIAPaCa2 cells revealed 115 metabolic genes or transporters that were differentially expressed between these two conditions ( Figure 8C and Table S12). Reactome Pathway Analysis of these 115 metabolic genes or transporters revealed that these differentially expressed genes were significantly enriched in pathways associated with the metabolism of carbohydrates, vitamins, and cofactors, among others ( Figure 8D and Table S13).
Next, we attempted to identify potential transcription factors that regulate differentially expressed metabolic genes. We performed Enrichr analysis to identify differentially expressed metabolic genes in the comparison between parental and oxamate-resistant MIAPaCa2 cells and in the comparison between vehicle-and oxamate-treated parental MIAPaCa2 cells. In the comparison between parental and oxamateresistant MIAPaCa2 cells, NFE2-like BZIP transcription factor 2 (NFE2L2), chromobox 3 (CBX3), TCF3, AR, and GATAbinding protein 1 (GATA1) were identified as the top-most significantly enriched transcription factors ( Figure 8E and Table  S14). In the comparison between parental and oxamate-treated MIAPaCa2 cells, CEBPB, zinc finger MIZ-type containing 1 A B C FIGURE 6 Metabolites generated during lipid, amino acid, and nucleotide metabolism in parental and oxamate-resistant cells treated with or without sodium oxamate. (A-C) Changes in the indicated metabolites generated during lipid, amino acid, and nucleotide metabolism under the indicated conditions.
In the comparison between vehicle-and oxamate-treated oxamate-resistant MIAPaCa2 cells, only 13 genes were identified as differentially expressed ( Figure 9A and Table S16); although these genes were not enriched in any major metabolic pathways ( Figure 9B), some of these genes were enrichment in unconventional metabolic pathways, such as cellular hexose transport and vitamin C (ascorbate) metabolism (Table S17). Finally, in the comparison between oxamate-treated parental and oxamate-treated oxamate-resistant MIAPaCa2 cells, 540 differentially expressed metabolic genes or transporters were identified ( Figure 9C and Table S18). Reactome Pathway Analysis of these 540 metabolic genes or transporters revealed the significant enrichment of pathways associated with the metabolism of carbohydrates, vitamins, and cofactors, among others ( Figure 9D and Table S19).
Next, we identified transcription factors potentially involved in the regulation of differentially expressed metabolic genes. To do so, we performed an Enrichr analysis of the differentially expressed metabolic genes in the comparisons between vehicleand oxamate-treated oxamate-resistant MIAPaCa2 cells and between oxamate-treated parental and oxamate-treated oxamate-resistant MIAPaCa2 cells. In the comparison between vehicle-and oxamate-treated oxamate-resistant MIAPaCa2 cells, GATA1 and repressor element 1-silencing transcription factor (REST) were identified as the top-most significantly enriched transcription factors ( Figure 9E and Table S20). In the comparison between oxamate-treated parental and oxamateresistant MIAPaCa2 cells, E2F transcription factor 6 (E2F6), Myc-associated factor X (MAX), upstream transcription factor (USF)1, USF2, and nuclear transcription factor Y subunit beta (NFYB) were identified as the top-most significantly enriched transcription factors ( Figure 9F and Table S21).

Integration of metabolic and RNA-seq data to identity co-regulated pathways at both the transcriptional and metabolic levels
To identify metabolic pathways and mRNA expression levels that are co-regulated, we performed an integrative analysis combining RNA-seq and metabolomics data. We focused on all metabolic genes identified as altered across the various comparisons (parental versus oxamate-resistant; vehicle-versus oxamate-treated parental; vehicle-versus oxamate-treated oxamate-resistant; and oxamate-treated parental versus oxamate-resistant). To conduct the integrated analysis, we used the joint pathway analysis tool in MetaboAnalyst 5.0 (31). The results of this integrated analysis are shown in Figures 10A-D and summarized in Tables S22-S25. For example, aminoacyl tRNA biosynthesis, arginine biosynthesis, and proline biosynthesis were identified as the three top-most significantly altered metabolic pathways in the integrated analysis for parental versus oxamate-resistant MIAPaCa2 cells. We also analyzed differences among comparisons to identify specific changes associated with specific pairwise comparisons. The Venn diagram for this analysis is shown in Figure 10E. For example, 75 metabolic pathways were identified as significantly altered in the integrated analysis for parental versus oxamate-resistant MIAPaCa2 cells, whereas 81 metabolic pathways were identified as significantly altered in oxamate-treated parental versus oxamate-resistant MIAPaCa2 cells. Of these, 75 were common between the two comparisons, 0 were specific to the comparison between parental and oxamate-resistant MIAPaCa2 cells, and 6 were specifically altered in the comparison between oxamate-treated parental and oxamate-resistant MIAPaCa2 cells. These results demonstrate that the integration of metabolomics and RNA-seq data can reveal information that is not possible to identify by looking at either RNA-seq or metabolomics data alone.

Discussion
Targeting metabolic enzymes has emerged as a new approach for cancer therapy (3,32). However, based on our knowledge of the response of cancer cells to other therapeutic agents, including the development of eventual therapeutic resistance, cancer cells are expected to eventually adapt to metabolic enzyme inhibition. Therefore, understanding the adaptation process at the molecular level is necessary to optimize the use of metabolic enzyme-targeting agents and improve the treatment or prevent the development of drug-resistant cancers.
The overexpression of the metabolic enzyme LDHA is necessary for tumor growth and progression (6)(7)(8)(9)(10)(13)(14)(15). LDHA typically exists as a tetramer (LDH-5) and functions primarily to interconvert pyruvate to lactate and transform NADH to NAD+. However, LDHA is a multifunctional protein and has been localized not only to the cytoplasm but also in mitochondria and the nucleus (33,34). Outside of the nucleus, LDHA plays a role in the regulation of glycolysis, whereas LDHA inside the nucleus functions as a singlestranded DNA-binding protein (SSB) (34). Additionally, nuclear LDHA has been shown to sense reactive oxygen species, leading to the production of alpha-hydroxybutyrate in human papillomavirus (HPV)-induced cervical cancer (35).
Here, using the LDHAi oxamate as a metabolic enzyme inhibitory drug and PDAC as a model system, we studied the impact of LDHA inhibition on the transcriptional, chromatin, and metabolomics landscapes. We also performed integrated analyses to identify common transcriptome and chromatin accessibility features and common transcriptome and metabolic features in parental and oxamate-resistant PDAC cells treated with and without oxamate. These studies revealed several differences in the mRNA expression profiles, chromatin accessibility profiles, and metabolic profiles across parental and oxamate-resistant PDAC cells treated with or without oxamate.
RNA-seq analysis showed that parental cells were highly responsive to oxamate treatment, whereas oxamate-resistant cells showed only marginal responses to oxamate treatment. Consistent with this finding, vehicle-and oxamate-treated parental cells displayed large differences in their mRNA expression profiles. These results are consistent with the rationale that once oxamate resistance develops, oxamateresistant cells become largely inert in response to oxamate treatment. In conclusion, these results demonstrate that oxamate resistant cells may acquire such a state as a mechanism of adaptation to survive in the presence of oxamate and to become oxamate resistant. Notably, also pathway analysis identified distinct set of genes and pathways gets enriched in each of treatment conditions, in which oxamate treated parental cells very clearly separated from highly nonresponsive oxamate resistant cells.
Furthermore, ATAC-seq analysis and its integration with the RNA-seq analysis revealed that only a small subset of mRNAs (1% to 12%) associated with differential expression were associated with changes in chromatin accessibility within promoter regions. These results are consistent with those of other studies, which also reported an approximately 10% overlap between RNA-seq and ATAC-seq data (36). Furthermore, the usefulness of integrated analyses was highlighted by the ability of integrated RNA-seq and ATAC-seq data to identify changes in chromatin accessibility associated with consequent changes in mRNA expression, which allowed for the identification of potential transcriptional regulators associated with changes in mRNA expression. Additionally, using the integration data with that of transcriptional regulator identification data allows for identifying distinct transcription factors that might contribute to distinct transcriptional response to oxamate treatment. Future studies will allow us and others to evaluate the role of these individual transcription factors in driving the oxamate resistant state and their role in glucose metabolism regulation in part via glycolysis regulation.
Finally, the integration of metabolomics data with RNA-seq data for metabolic genes and transporters allowed changes in mRNA levels to be linked with changes in metabolites. These studies thereby provide an overall reason for changes in metabolites under distinct oxamate treatment or parental and oxamate-resistant states, and demonstrate that they in part are driven due to changes at transcriptional of enzymes that in turn regulate various metabolic pathways.
In conclusion, these studies highlight the importance of performing integrative analyses of RNA-seq, ATAC-seq, and metabolomics data to identify changes that would not be possible to identify using any single omics-based approach. Further studies using other functional genomics approaches are likely to reveal the driving roles played by these transcriptional and metabolic changes in conferring resistance to LDHAis. Similar approaches can be employed to identify differences in the response to other metabolic inhibitors, which may eventually lead to the optimization of metabolic enzyme-targeted therapies to achieve the maximum clinical benefits in cancer patients.
Future studies using a larger number of PDAC cell lines, in vivo models, and other cancer types of glycolytic nature will allow us to further identifying unifying features of oxamate resistant cells. Similarly, use of other LDHA inhibitors and comparing the genetic inhibition of LDHA with that or its pharmacological inhibition will reveal the impact of nonenzymatic and non-canonical functions of LDHA in promoting growth of cancer cells.

Cell culture
MIAPaCa2 cells were obtained from the American Type Culture Collection (Manassas, VA, USA) and were grown in Dulbecco's modified Eagle medium (Life Technologies, ThermoFisher Scientific, Waltham, MA, USA), supplemented with 10% fetal bovine serum (Life Technologies, ThermoFisher Scientific) and 1% penicillin/streptomycin (Life Technologies), at a CO2 concentration of 5%.

Generation of sodium oxamate-resistant cell lines
Sodium oxamate-resistant cells were generated by continuously culturing MIAPaCa2 cells in the presence of 10 mM sodium oxamate for 4 weeks. The medium was replaced every 3 days and supplemented with fresh sodium oxamate until these cells began to form distinct colonies. The sodium oxamateresistant phenotype in isolated polyclonal populations was further validated by clonogenic assay.

Clonogenic assays
The clonogenic ability of MIAPaCA2 parental and sodium oxamate-resistant cells was measured under untreated and sodium oxamate-treated conditions. For these assays, 5 × 10 3 cells/well were seeded in 6-well culture plates. For drug treatment experiments, cells were treated with 10 mM of sodium oxamate. The medium containing sodium oxamate was changed every 3 days. After 10-14 days of treatment, surviving colonies were stained with a solution containing 40% methanol, 10% acetic acid, and 0.005% Coomassie Brilliant Blue R-250 (Sigma-Aldrich, USA), and the plates were imaged using an Epson Perfection V850 Pro Photo Scanner (USA).
Seahorse analysis for the measurement of extracellular acidification rates Extracellular acidification rates were measured using the XFe24 analyzer (Seahorse Bioscience). Briefly, parental and oxamateresistant MIAPaCa2 cells (untreated and treated with sodium oxamate) were seeded at a density of 40,000 cells/well in triplicate. After 24 h, the culture medium was replaced with assay medium, and extracellular acidification rates were measured for 2 h.

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, Carlsbad, CA, USA), 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(dT) 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 (Agilent Technologies, Palo Alto, CA, USA) and quantified using a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA) 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 the entire read sequences. BAM files were generated during this step. Unique gene hit counts were calculated by using feature counts from the Subread package v.1.5.2. Only unique reads that fell within exon regions were counted.
The reads were first mapped to the latest UCSC transcript set using Bowtie2 version 2.1.0 (37), and the gene expression level was estimated using RSEM v1.2.15 (38). Differentially expressed genes were identified using the DESeq2 program (39). Genes showing altered expression with p < 0.05 and fold changes > 1.5 were considered to be differentially expressed. Goseq (40) was used to perform the gene ontology (GO) enrichment analysis, and Kobas was used to perform the pathway analysis (41). The reads were first mapped to the latest UCSC genome set using Bowtie2 version 2.1.0 (37). Mitochondrial reads, duplicate reads, and non-unique reads were removed before peak calling. MACS2 was used for peak calling using BAMPE mode (42). Differentially expressed peaks were identified using the DEseq2 program (39). Peaks showing altered expression with p < 0.05 and fold change > 1.5 were considered differentially expressed.

ATAC-sequencing and data analysis
Downstream genes of the differential peaks were used for GO and pathway enrichment analysis. Goseq (40) was used to perform the GO enrichment analysis, and Kobas was used to perform the pathway analysis (41).
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 Reactome Pathway Analysis.

Enrichment analysis for transcription factor binding sites
Transcription factor enrichment analysis was performed by Enrichr, a web-based enrichment analysis using ENCODE and the ChEA consensus TFs from the ChIP-X option for identifying potential regulatory transcription factors involved in the differential expression of metabolic genes between parental and oxamate-resistant MIAPaCa2 cells treated with and without sodium oxamate. The gene set of differentially expressed genes was used as input. To identify transcription factors using data from multiple experiments, set intersection was applied to obtain consensus. For each comparison, multiple transcription factors were identified as potentially direct transcriptional regulators for differentially expressed mRNAs.

Reactome pathway analysis
The differentially expressed genes identified in the RNA-seq or ATAC-seq analyses of metabolites identified by metabolomics analysis were evaluated using Reactome Pathway Analysis using the Reactome Pathway Analysis tool version 80 (https://www. reactome.org). Briefly, common differentially expressed genes among parental and oxamate-resistant MIAPaCa2 cells treated with and without sodium oxamate were uploaded into the Reactome Pathway Analysis tool. In this analysis, all nonhuman identifiers were converted to their human equivalents. The 25 most significantly altered pathways specific to Homo sapiens were identified using all resources and sorted by p-values.

Metabolomics analysis
Parental and oxamate-resistant MIAPaCa2 cells treated with and without oxamate were analyzed for alterations in metabolic pathways using the CE-TOFMS-based scan profiling method developed by Human Metabolome Technologies (Boston, MA, USA) using an Agilent CE-TOFMS system (Agilent Technologies Inc.). Parental and oxamate-resistant cells were grown in the absence and presence of 10 mM sodium oxamate in triplicate, and 1 × 106 cells were analyzed for each condition. Samples were prepared according to the recommendations of Human Metabolome Technologies. The samples were mixed with 800 µL methanol, followed by the addition of 450 µL Milli-Q water containing internal standards (10 µM), and mixed thoroughly. The extract (1,000 µL) was centrifuged (2,300 × g, 4°C, 5 min). The supernatant (350 µL) was filtrated through a 5-kDa cutoff filter (ULTRAFREE-MC-PLHCC, Human Metabolome Technologies, Yamagata, Japan) to remove macromolecules. The filtrate was concentrated by centrifugation and resuspended in 50 µL ultrapure water immediately before the measurement. Metabolome analysis was performed in samples of cultured cells using CE-TOFMS for cationic and anionic metabolites. For data analysis, peaks detected during spectrometric analysis were extracted using the MasterHands version 2.17.1.11 automated integration software (developed at Keio University, Tokyo, Japan) to determine the mass/charge ratio (m/z), migration time, and peak area. Peak area was converted to relative peak area using the following equation: Relative peak area = Metabolite Peak Area/ Internal Standard Peak Area × Sample Amount. The peak detection limit was determined based on a signal-to-noise ratio of 3. Putative metabolites were assigned based on the m/z value and migration time using Human Metabolomic Technologies' standard library and known-unknown peak library. The tolerance was ±0.5 min for migration time and ±10 ppm for m/z. Hierarchical cluster analysis and principal component analysis were performed using statistical analysis software (developed at Human Metabolome Technologies). In total, 228 metabolites were detected (121 metabolites in cation mode and 107 metabolites in anion mode) and annotated using Human Metabolomic Technologies' standard library and knownunknown peak library. All metabolite concentrations were calculated by normalizing the peak area of each metabolite to the area of the internal standard and by comparing with standard curves obtained from a 100 µM single-point calibration. The peak profiles of putative metabolites were represented on metabolic pathway maps using the Visualization and Analysis of Networks containing Experimental Data (VANTED) software (http://vanted. ipk-gatersleben.de/). The pathway map was prepared based on metabolic pathways known to exist in human cells according to information in the KEGG database (http://www.genome.jp/kegg/). The list of fold changes for altered metabolites identified in each comparison is shown in Table S8.
Integrated analysis for RNA-seq and metabolomics data Differential genes and metabolites with fold changes were combined and used as the input file in the MetaboAnalyst 5.0 joint pathway analysis tool. MetaboAnalyst supports raw MS spectra processing, comprehensive data normalization, statistical analysis, functional analysis, meta-analysis, and integrative analysis with other omics data. The most significantly enriched pathways were selected and plotted.

Data availability statement
The RNA-seq data and ATAC-seq data was submitted to Gene Expresison Omnibus. (GEO). The accession number for the RNA-seq data is GSE164976 and the accession number for ATAC-seq data is GSE164974. Metabolomics data and other relevant spread sheets can be found in the article/ Supplementary material.

Author contributions
NW and RG conceived the project. PM, VR, RG, and NW designed the experiments. PM and VR performed the research. PM, VR, RG, and NW analyzed and interpreted the data. PM, RG, and NW wrote the paper. All authors reviewed, made suggestions, approved the paper, and provided comments.