Divergent organ-specific isogenic metastatic cell lines identified using multi-omics exhibit differential drug sensitivity

Background Monitoring and treating metastatic progression remains a formidable task due, in part, to an inability to monitor specific differential molecular adaptations that allow the cancer to thrive within different tissue types. Hence, to develop optimal treatment strategies for metastatic disease, an important consideration is the divergence of the metastatic cancer growing in visceral organs from the primary tumor. We had previously reported the establishment of isogenic human metastatic breast cancer cell lines that are representative of the common metastatic sites observed in breast cancer patients. Methods Here we have used proteomic, RNAseq, and metabolomic analyses of these isogenic cell lines to systematically identify differences and commonalities in pathway networks and examine the effect on the sensitivity to breast cancer therapeutic agents. Results Proteomic analyses indicated that dissemination of cells from the primary tumor sites to visceral organs resulted in cell lines that adapted to growth at each new site by, in part, acquiring protein pathways characteristic of the organ of growth. RNAseq and metabolomics analyses further confirmed the divergences, which resulted in differential efficacies to commonly used FDA approved chemotherapeutic drugs. This model system has provided data that indicates that organ-specific growth of malignant lesions is a selective adaptation and growth process. Conclusions The insights provided by these analyses indicate that the rationale of targeted treatment of metastatic disease may benefit from a consideration that the biology of metastases has diverged from the primary tumor biology and using primary tumor traits as the basis for treatment may not be ideal to design treatment strategies.


Introduction
Breast cancer is the most common malignant neoplasm among women in the United States. Recently, the American Cancer Society reported a 5-year survival rate near 86-99% for regional and local breast cancer, respectively [1]. On the other hand, the 5-year survival for metastatic breast cancer that involves distant organs is only 27% [1]. This latter low survival rate is likely due in part to a difference in clonal divergence of the metastatic tumor growth as well as to the use of primary tumor characteristics as the rationale for treatment strategies of metastatic disease (e.g., [2]). Thus, there is a growing consensus that matched primary breast tumor and metastatic lesion biopsy samples often exhibit divergent expression of markers, for example, hormone receptors (HR: estrogen (ER) & progesterone (PR)) as well as HER2, which influences outcomes [3][4][5][6][7][8]. In addition, genomic sequencing studies are providing strong corroborating evidence that metastatic progression represents evolutionary processes that results in distinct biological entities at metastatic sites that have greatly diverged from the primary tumor [9][10][11]. Consequently, the current practice of primary tumor-based selection of chemotherapy is limited with respect to patient specific precision therapeutic targeting of a patient's metastatic lesions as well as general therapeutic resistance.
Notably, monitoring and treating metastatic progression remains a formidable task due to many gaps in our knowledge including an inability to monitor specific differential molecular adaptations that allow the cancer to survive and thrive within different tissue types. This is a consequence of the fact that visceral organs differ vastly from one another with unique attributes of metabolism, developmental programs, microenvironments, and function resulting in defined physiological identities. Hence, in order to develop optimal treatment strategies for metastatic disease, an important consideration is the divergence of the metastatic cancer growing in visceral organs from the primary tumor [2,4,5,11]. Accordingly, in order to gain insight into new treatment regimens aimed at controlling and ablating metastatic progression, there is an urgent need to evaluate the distinct molecular differences that exist between isogenic tumor cells growing at different metastatic sites and their sensitivity to different chemotherapeutic agents. To address this issue, as previously reported [12,13], we established isogenic (analogous to patient metastatic cells) human metastatic breast cancer cell lines and have now included two additional cell lines starting from a well-established aggressive breast cancer cell line (MDA-MB-231). In all cases, these metastases spontaneously arose through dissemination from the primary mammary fat pad tumor site in a mouse model system. The resulting six metastatic cell lines are representative of the common metastatic sites of lymph node, lung, bone, liver, and brain observed in breast cancer patients [14]. This model system has provided excellent data that supports our principal hypothesis that organ-specific growth of malignant lesions is a selection process that results in cancer cells that have adapted distinct biochemical and molecular attributes that allow them to thrive at each site outside of the primary tumor site [12,13].
Dissecting metastatic cancers based on objective molecular markers and metabolic pathways remains an outstanding challenge. Standard histochemical techniques are limited to identification of relatively few markers in, most frequently, primary tumor sections that may no longer be present at metastases [4][5][6][7][8]. Consequently, a goal of the presented study was to obtain proteomic and RNAseq data sets for our isogenic metastatic cell lines and analyze the resulting proteomes and transcriptomes with pathway analysis tools. Analyses of both data sets revealed that the metastatic cell lines had diverged from the primary tumor, consistent with our previous studies [12,13]. We found, relative to the primary tumor, overlapping general common metastatic associated pathways, but importantly unique organ-specific pathways have also been uncovered. We hypothesize that the latter reflect processes of adaptation, i.e., the gaining of site (organ) specific attributes that reflects local microenvironmental influences resulting in selected gene expression and protein pathway preferences for each organ. Similar, to previous studies of comparisons of proteomic and RNAseq data [15][16][17][18], we discovered discordance between proteomes and transcriptomes as well as some common characteristic pathways. Consistent with the RNAseq data, RT-PCR, in general, only provided confirmation of very few of the expressed proteins discovered during proteomic screening. In addition, pathway analyses of metabolomes provided confirmations corresponding to associations to the proteomic-based and transcriptomic-based pathways. Finally, in vitro drug efficacy assays showed significant differential responses of ten cell lines, i.e., two parental cell lines, two primary tumor cell lines, and six metastatic cell lines, to the drugs that were tested. This latter data has provided evidence that chemotherapeutic regimes based on primary tumor markers may result in ineffective control of metastatic tumors due to the changes that have occurred within the tumor cells at those metastatic sites that affect a drug's killing ability at the site.

Phenotypic characterization of isogenic cell lines
As we previously reported for the isogenic cell lines generated from the MDA-MB-435 cells: parental-435 and its primary-tumor (1˚-tumor), brain, liver, lung, and spine isogenic cell lines (S1 Fig) [12], the parental MDA-MB-231 (parental-231) and its isogenic cell lines (1˚-tumor, lung, and lymph node) exhibited somewhat different morphologies when grown on plastic (All bright-field images: Fig 1 and S2 Fig). As has been reported [19], parental-231 cultures had a mixture of several cell morphologies that differed in shape and size and that grew in a chaotic overlapping manner (Far-left bright-field images: Fig 1 and S2 Fig). 1˚-tumor-231 cells grew in a similar manner but appeared to have fewer cell morphologies with the majority of the cells in these cultures being relatively large, broad, and elongated with some extended spindle characteristics and few or no small cells (Middle-left bright-field images: Fig 1 and S2 Fig).
In contrast, lung-231 cell cultures had even fewer small cells and the majority were medium sized spindly cells with some broadly larger elongated cells (Middle-right bright-field images: Fig 1 and S2 Fig). On the other hand, the majority of lymph node-231 cells appeared to be relatively more epithelial-like and formed "cobble-stone"-like monolayers with some overlapping board-spindly cells (Right bright-field images: Fig 1 and S2 Fig). Also, reflected in the brightfield images of the four cultures, which were plated at the same time and in the same numbers, were differences in growth-rates between the lung-231 and lymph node-231 cell lines and between these two cell lines relative to the parental and 1˚-tumor cell line (Lower panels of Fig  1). Thus, the lung-231 cells had the lowest growth-rate (36 hr: Lower panels of Fig 1), which is visible as a sparser covering of the plastic in the bright-field images, along with an apparently higher death-rate/senescence as indicated by the stable or plateauing/declining growth-rate by day 4 and 5. In contrast, the lymph node-231 cells had the fastest growth rate (22 hr: Lower panels of Fig 1) and these cultures (in the bright-field images) appeared to be nearly confluent relative to the other three cultures. We previously reported similar differential growth-rates in growth-rate comparisons for cultures of 435 isogenic cell lines (S3 Fig) [12]. Overall, the characteristics of the in vitro growth of all the isogenic cell lines has indicated that fundamental underlining inherent molecular changes had occurred during in vivo growth and, at least during these early passages on plastic, these differential characteristics endured as reflected in differing morphologies and growth-rates.

General "omics"
Other than for the principle component analysis (PCA) and hierarchical clustering analyses (indicated below), we focused on the expression level fold changes (FCs: range -1.25 to +1.25) of the individual proteins or transcripts (genes) or aqueous metabolites from each isogenic metastatic cell line relative to their cognate counterparts from the 1˚tumors. Our hypothesis is that metastatic progression, starting from the 1˚tumor, is an evolutionary process that evolves in situ at each specific tissue site under the influence of inherent microenvironmental signaling and site-specific selection pressures. As such, we have been interested in understanding the Phase-contrast images of the parental-231, primary tumor (1t umor)-231, lung-231, and lymph node-231 cell lines are shown in the top two rows. The top row images were photographed using a X10 objective coupled with a X4 phase-contrast ring while the second-row images were duplicate images obtained with a X10 objective coupled with a X10 phase-contrast ring. The optical configuration of the top row gave 3D images. The black scale bars = 50 μM. Growth curves of each cell line are presented below the images of the corresponding cell line with the mean growth-rate given in the bottom right-hand corner of the curves. https://doi.org/10.1371/journal.pone.0242384.g001

Proteomic-based PCA and proteomic-transcriptomic-based hierarchical clustering
Similar to our previous metabolomic and Raman spectroscopic based PCA and hierarchical clustering analyses, proteomic-based PCA indicated that 231 isogenic cell lines (Fig 2A) as well as 435 isogenic cell lines ( Fig 2B) were separated into distinct tissue defined clusters, which was an indication that each isogenic cell line differs from its isogenic counterparts at the proteome level. The heat map shown in Fig 2C is complementary evidence that all isogenic cell lines have distinct proteomes. Thus, similar to the relative distances indicated in Fig 2B and 2C indicates that the proteome of parental-435 cells was most closely related to that of the 1˚tumor-435 with the latter being more closely related to the liver-435 cell line, while the proteomes of the spine-435 and lung-435 cell lines formed a subclade and that brain-435 cells had a proteome that was least related to the other isogenic family of cell lines. However, Fig 2C shows that the parental-231 and 1˚tumor-231cells formed a subclade, which was not apparent from Fig  2A, but, as in Fig 2A, the lymph node-231 (LN-231) and lung-231 cell lines were closely related relative to their distance from the parental-231 and 1˚tumor-231 cell lines' subclade. From Fig  2C it was found that the proteomes of the 435 cell lines and 231 cell lines formed two distinct general separate clades that was likely due to their parental cells being from different individuals. This is also likely why the two lung proteomes were not closely related, which provided some evidence that proteomes from distinct individuals remain generally discrete even after growth within very similar microenvironments.
The transcriptomic-based clustering analysis depicted in Fig 2D is complementary to Fig  2C and again indicates that the 435 cell lines and 231 cell lines formed two distinct general separate clades. In addition, the two 231 isogenic metastatic cell lines (lung-231 and lymph node (LN)-231) had, similar to their proteomes, transcriptomes that were closely related while in this case, unlike the proteomic-based analysis (Fig 2C), the parental-231 and 1˚tumor-231 cell line transcriptomes were more distantly related. A similar distinction between proteomicbased and transcriptomic-based clustering was exhibited by the parental-435 and 1˚tumor-435 cell lines with the latter forming a subclade with brain-435 and these three grouped separately from the lung-, liver-and spine-435 cell lines with the latter two of these forming a subclade. Hence, although proteomic-based and transcriptomic-based PCA/hierarchical analyses gave consistent complementary result with respect to providing evidence that each cell line display distinct transcriptomes and proteomes there was not an exact match in clustering patterns Transcriptomic-based hierarchal clustering of all cell lines. All analyses indicated that each cell line had distinct proteomic/transcriptomic signatures, which resulted in the cell lines' clustering into separate groups/clades. As shown beneath the heat maps, proteins (C) or transcripts (D) distributed across rows have been presented between to the two data sets. This was consistent with the known phenomena that proteomic and transcriptomic data sets do not generally exhibit large expression overlaps between transcripts with their protein products [17,18,20], which has been described as being due to a variety of regulatory distinctions associated with mRNAs and proteins [15,16,18,[20][21][22][23][24][25].
transport' (S2 Table), 'metabolism' as well as 'hemostasis' (S3 Table), 'cellular responses to stress' along with 'metabolism of carbohydrates' (S4 Table), and 'metabolism of RNA', 'ribosome' and 'innate immune system' (S5 Table). Nevertheless, many pathways listed (S1-S6 Tables) were relatively specific. Moreover, the pathway analyses allowed for an inspection of up-and down-regulated pathways that were common to 2 or more metastatic sites and thus, provided a means to assess those pathways that promote general metastatic processes, regardless of being up-or down-regulated. The abridged datasets already include several examples of shared pathways, such as up-regulated 'spliceosome' of brain-435 and lymph node-231 and very closely related 'mRNA splicing-major pathway' along with 'mRNA splicing' of lung-231, up-regulated 'lysosome' brain-435 and liver-435, down-regulated 'cell cycle' of brain-435 and liver-435, down-regulated 'pyrimidine metabolism' brain-435 and liver-435, up-regulated 'TCA cycle' of brain-435 and the very closely related 'TCA cycle & respiratory electron transport' of lung-435, up-regulated 'metabolism of RNA', 'cell cycle', and 'cell cycle, mitotic' of lung-231 and lymph node-231, and down-regulated 'EGFR1', 'neutrophil degranulation', 'metabolism', 'vesicle-mediated transport', 'membrane trafficking' as well as 'posttranslational protein phosphorylation' of lung-231 and lymph node-231. Also, there were examples of pathways that were up-regulated at one site while being down-regulated at other sites. Examples, (S1-S6 Tables) are: up-regulation of 'insulin-like growth factor (IGF) transport & uptake by IGF binding proteins (IGFBPs)' and 'neutrophil degranulation' in liver-435 and the down-regulation of both pathways in lung-435 and the former is in lymph node-231 and the latter is in spine-435, lung-231, and lymph node-231.
At the same time, pathway analysis allowed for the discovery of up-and down-regulated pathways that are unique to each metastatic site, which provides insights into how the cells may have evolved by adapting attributes of their site of growth that would have been induced by signaling cascades that were inherent to each tissue. Similar to the analysis of all pathways given above, the top 10 unique up-regulated and top 10 unique down-regulated proteome based pathways for each isogenic metastatic cell line, ranked on the lowest to highest q-values (q � 0.05), are provided in Tables 1-6 with the complete pathway lists given in S19 Spreadsheet. The general top 10 summaries shown in Tables 1-6 indicate that: in brain-435 cells the 'mitochondrial TCA' pathway was up-regulated, which may have occurred either as a response to energy needs, i.e., as an energy associated pathway via metabolism of glucose to pyruvate (see the listed superpathway) and then to acetyl-CoA for use in the TCA cycle or was up-regulated in response to an anaplerosis need, the 'mitochondrial fatty acid β-oxidation' (energy generating) pathway was also up-regulated while a down-regulation of mitochondrial biogenesis and DNA repair pathways was found for brain-435; in the liver-435 cell line relatively liverspecific pathways were up-regulated including 'fibrin clot formation', general 'hemostasis', heparan sulfate-glcNAc-glcA (HS-GAG) degradation along with scavenger pathways while 'mitotic checkpoint', glucose uptake, apoptosis, and 'RNA metabolism' pathways were downregulated; in lung-435 cells components of the innate immune system (lung, similar to skin, must safe-guard against airborne pathogens), i.e., interferon signaling pathways, were up-regulated along with oxidative phosphorylation pathways, which were likely induced by the high O 2 tension of the lung, as well as disease pathways, such as Parkinson's, Alzheimer's, etc., which have been associated with mitochondrially generated reactive oxygen species while down-regulations included the translation (EIF-4e) but proapoptotic (p70s6) pathway, the hypoxia driven VEGF pathway, which was likely due to the high O 2 tension of the lung as well as estrogen and androgen signaling; the spine-435 cell line had an up-regulation of histone modification, amino acid and oligopeptide solute transport, 'terpenoid backbone biosynthesis', 'FOXA1 transcription factor network', with a down-regulation of several toll-like receptor (innate immune response) pathways; the pathways up-regulated in lung-231 cells were several involved with translation including mitochondrial translation and DNA repair, with downregulated pathways involved with extracellular matrix degradation, apoptosis, IL-7 (hematopoietic growth factor), epidermal growth factor (EGF) signaling, and stress induced heat shock protein; in the lymph node-231 cell line RUNX3 transcriptional regulation and HDAC Class I signaling were up-regulated along with up-regulated 'hematopoietic stem cell regulation by GABP-α/β complex', which might reflect influences on the cells during lymph node site growth as these pathways are involved with regulation of hematopoietic lineages, RHO GTPases, 'EPHA-mediated growth cone collapse', 'TGF-β receptor', and 'glucocorticoid receptor regulatory network' were also up-regulated with 'mitochondria β-oxidation of short chain saturated fatty acids', 'urea cycle and metabolism of arg, pro, glu, asp, & asn', and one-carbon metabolism were down-regulated.
In order to gain an appreciation that the lists of proteome pathways are not isolated independent biochemical reaction entities but rather are interconnected systems, we leveraged the ConsensusPathDB analysis program to construct overlapping pathway interconnection maps [26,27]. This analysis revealed that due to the multifunctional attribute of one or more of the proteins of a pathway two or more pathways have shared proteins that connect pathways into larger networks of biochemical reaction systems with overlapping functions. Examples of such maps are shown in In this example, one of the interconnected maps (Left-hand side, Fig 4) was dominated by heat shock factor (HSF1) associated pathways with a relatively high amount of overlapping protein components and is loosely connected with an insulin/growth hormone signaling group of pathways that have fewer overlapping protein partners. On the right-hand side of the upper portion of Fig 4 there emerged a clustered network of strongly associated TCA pathways  Transcriptome-based pathway discovery. Analogous to the proteome-based analyses, lists of transcripts from each isogenic metastatic cell line with FCs � -1.25 and � +1.25, relative to their 1˚tumors (S7-S12 Spreadsheets) were analyzed with ConsenusPathDB. Abridged data sets, i.e., the top 10 up-and down-regulated pathways for each isogenic metastatic cell line, (q � 0.05) are provided in S7-S12 Tables with the complete pathway lists given in S20 Spreadsheet. Although transcriptome-based pathway discovery analyses uncovered pathways common with those of the proteome-based pathways (see below), overall, the majority of transcriptome-based pathways differed from the above proteome-based pathways.
Similar to the proteomic based pathway lists, we submitted the transcriptomic based lists of up-and down-regulated pathways to hierarchical clustering analysis and the resulting heat maps are shown in Fig 3B. Unlike the proteomic-based heat maps, both the up-and down-regulated pathway data sets of these analyses separated the lung-231 and lymph node-231 cell lines. In the case of the up-regulated heat map the lung-231 cell line was linked to a clade made up of closely related spine-435 and lymph node-231 cell lines while at the same time being closer to the lung-435 cell line that formed a clade with the brain-435 cell line with this clade being more distantly linked to the liver-435 cell line. In the case of the down-regulated heat map, lung-231 was linked to the clade of closely related spine-435 and liver-435 cell lines while the lymph node-231 cell line was linked in a separate clade to the liver-435 cell line and more distantly, both were linked to the brain-435 cell line. As with the ConsensusPathDB analyses of the proteomic data sets, several of the integrated databases of ConsensusPathDB uncovered the identical or overlapping similar pathways from the submitted gene lists. For example, in metastatic brain-435 (S7 Table) the 'ECM-receptor interaction' pathway (Source: KEGG) was up-regulated and confirmed as the 'extracellular matrix organization' pathway (Source: Reactome) and (S7 Table) the down-regulated 'DNA replication' pathway (Source: Reactome) was replicated (Source: Wikipathways). Several of these types of examples can be seen in the S7-S12 Tables. General vague pathway designations were also again observed, such as 'extracellular matrix organization', 'cell cycle', 'cytokine signaling in immune system', 'metabolism of RNA', 'gene expression (transcription)', 'RNA polymerase II transcription', generic transcription pathway', 'chromatin modifying enzymes', 'axon guidance', 'muscle contraction', and 'neutrophil degranulation'.
The top 10 up-and down-regulated transcriptome-based unique pathways for each isogenic metastatic cell line are listed in the S13-S18 Tables. Pathway duplications or lists of similar pathways from different integrated Sources along with vague pathways designations were again obvious. Pathways that provided possible examples of metastatic cell lines acquiring tissue specific assimilations included: for brain-435 up-regulated 'presynaptic depolarization & calcium channel opening', 'NCAM1 interactions', and phenylethylamine degradation I' pathways; for liver-435 up-regulated 'interferon-α/β signaling' and 'chondroitin sulfate/dermatan sulfate metabolism' pathways; for spine-435 (given its neuronal component) up-regulated 'striated muscle contraction', 'val, leu, and Ile degradation', 'ketogenesis', and 'ion channel transport' pathways; and, for lymph node-231 up-regulated 'non-genomic actions of 1,25-dihydroxyvitamine D3'. Similar to the case of the proteome pathways, we explored the interconnected pathway networks formed by those transcriptome-based pathways that were unique to each isogenic cell line. Examples of such networks are shown in Fig 6 (brain-435) and Fig 7 (lung-231) as well as in S8-S11 Figs. However, in these cases the number of pathways that were significantly up-or down-regulated, i.e., q � 0.05, was not always greater-than or equal to 20. Thus, for those cases with less-than 20 pathways: brain-435 cells had 9 up-and 7 down-regulated pathways (  (Fig 7), and lymph node-231 cells had 8 down-regulated (and 5 of these were trending down with q = 0.0.55; S11 Fig). Not only were there often fewer significant (q � 0.05) transcriptomic pathways in the networks but, in all cases, the connections between the pathways within the networks were weaker, that is, the connections consisted of fewer shared transcript products. Moreover, several more 'orphan', i.e., non-connected pathways emerged during these analyses.
Common proteomic and transcriptomic pathways. As a means of confirming the proteome-based and transcriptome-based pathway discovery analyses, we compared the protein lists to their corresponding transcript lists (FCs � -1.25 and � 1.25) to match the proteins to their identical transcripts. This showed that the resulting matched protein-transcript lists were generally only a small percentage of either of their source lists, i.e., RNAseq-based transcripts  10.0 or 6.3%, 17.8 or 5.3%, 12.7 or 6.8%, 38.9 or 16.9%, and 34.5 or 17.4% for brain-435, liver-435, lung-435, spine-435, lung-231, and lymph node-231 respectively. These resulting lists were submitted to ConsensusPathDB to find those pathways that were common to both "omics" analyses. Representatives of these common pathways are given as the top 10 up-and down-regulated pathways (q � 0.05) in S19-S24 Tables with the complete common pathway lists given in SPT1 Spreadsheet. In the case of the lung-435 cell line we found only 5 up-regulated pathways that trended (S21 Table-gray shading: q = 0.055) as common. In addition, during all pathway discovery analyses, we used the default settings of Con-sensusPathDB, which used � 2 proteins to define a 'pathway'. Consequently, several of the listed common pathways in S19-S24 Tables were composed of only a few (2-4) proteins/genes and yet remained statistically significant, i.e., q � 0.05, and therefore, have been included but their relevant biological significance will require further investigations. Common as well as unique proteomic-transcriptomic pathways. The top 10 up-and down-regulated pathways (q � 0.05) that were common as well as unique proteome-transcriptome pathways of each isogenic metastatic cell line are presented in Tables 7-12. Examples of possible tissue specific associated pathways were those of liver-435: the up-regulated general 'immune system' and especially the 'innate immune system', 'neutrophil degranulation', 'MHC class II antigen presentation', and general 'hemostasis', as well as with 'metabolism of

PLOS ONE
Divergent pathways of isogenic cell lines alter drug efficacies

PLOS ONE
Divergent pathways of isogenic cell lines alter drug efficacies

PLOS ONE
Divergent pathways of isogenic cell lines alter drug efficacies fat-soluble vitamins'. The immune system and its sub-pathways that are involved with exogenous pathogen threats are prevalent components of the liver, which receives the largest portion of its blood supply from the portal vein and hence the gut where pathogen levels are relatively high [28]. In addition, the liver processes and stores fat-soluble vitamins [29]. For spine-435 cells such a pathway was the up-regulated 'cholesterol biosynthesis, regulation, and transport' pathway, which is important for bone growth hemostasis as well as neuronal function [30,31]. Surprisingly, reported tumor/metastatic suppressor pathways: 'transcriptional regulation by RUNX3', 'death receptor signaling', and 'regulation of PTEN stability and activity' were upregulated in lymph node-231 cells and as such deserve further investigation as to the function of these pathways in normal lymph nodes/lymphocytes/hematopoiesis, and reticular cells.
Metabolomic-based pathway discovery. As a further means of confirming the pathways revealed by ConsensusPathDB when using the protein or transcript lists from the proteomic and RNAseq studies we submitted lists of aqueous phase metabolites that were acquired from an earlier metabolomics study of the 435 cell lines, which included hierarchical clustering and principle component analyses but no pathway analyses [12], to ConsensusPathDB. The resulting complete metabolomic pathway discovery and comparison pathways along with the metabolite CAS numbers are given in S21-S28 Spreadsheets. The top 10 up-and down-regulated pathways (q � 0.05) for each of the isogenic metastatic cell lines are presented in S25-S28 Tables. For the brain-435 and liver-435 cell lines the up-regulated pathway lists were generally based on a small number of metabolites (2-4 metabolites) as defining the pathways, which, were statistical significant, but as noted above (Common proteomic and transcriptomic pathways) their biological relevance needs further investigation. Also, these analyses revealed a possible limitation to the CensusPathDB platform as it was found; e.g., in the upregulated pathway list of lung-435 cells that the SMPDB Source gave 5 apparent independent pathways, but these were based on the same 6 metabolites (S25 Spreadsheet). Thus, as noted above during the proteomic-based pathway discovery analyses, online pathway designations Table 7. Unique pathways found to be common in both proteomic and transcriptomic analyses for the metastatic brain-435 cell line. can result in lists of similar related or identical pathways that have been given different labels, which appears to have been the case with this example of SMPDB designations. The latter was again apparent for all 4 lists of down-regulated pathways shown in S25-S28 Tables. Metabolomic-based unique pathways. In these cases (S29-S32 Tables), the majority of pathways were either considered as being weakly defined (2-4 metabolites/pathway) or the different pathway designations were based on identical metabolite lists; e.g., Source SMPDB in the down-regulated list for the liver-435 cell line (S30 Table and S24 Spreadsheet) and in the up-and down-regulated lists of lung-435 (S31 Table and S25 Spreadsheet).

Source Up Pathways # of Gs-Ps ǂ in Set # of Obs. Gs-Ps Obs. Gs-Ps (%) q-value
Common metabolomic and proteomic pathways. An analysis of these common pathways revealed only two (e.g., in brain-435, lung-435, and spine-435) or no (liver-435) up-regulated pathways but those that were found were unique to each cell line (S33-S36 Tables). In addition, the two up-regulated transport pathways (S36 Table) of spine-435 may indicate some specific influence of the neuronal component of spine during growth in the spine. On the other hand, most of the top 15 common down-regulated pathways (S33-S36 Tables) were shared between cell lines. For example, the 'cell cycle, mitotic' (Source: Reactome) pathway was down in brain-, liver-, and spine-435 cell lines while the 'pyrimidine metabolism" (Source: KEGG) pathway was down in brain-, liver-, and lung-435 cell lines. Those pathways that were down-regulated in brain-435 and liver-435 cell lines included: 'cell cycle', 'translation', as well as 'S phase' (Source: Reactome) along with the 'pyrimidine metabolism' (Source: Wikipathways) pathway, and 'purine metabolism' (Source: KEGG) pathway. The 'post-translational protein modification' (Source: Reactome) pathway was shared by brian-435 and lung-435 while the 'metabolism of nucleotides' (Source: Reactome) pathway was shared by brain-435 and spine-435. The vaguely defined 'metabolism' (Source: Reactome), 'pyrimidine nucleotides nucleosides metabolism (Source: INOH) along with 'DNA replication' and 'selenoamino acid metabolism' (Source: Reactome) pathways were unique to brain-435. This type of pattern of shared pathways can be seen when S33-S36 Tables were compared. Unique pathways for liver-435 (S34 Table) were: the 'superpathway of purine nucleotides salvage' (Source: Human-Cyc) pathway along with 'DNA replication', 'teleomere C-strand (lagging strand) synthesis', and 'TCA cycle and respiratory electron transport' (Source: Reactome) pathways. For the lung-435 cell line (S35 Table) the unique pathways were the: 'pentose phosphate pathway', 'glycogenosis, type IA, von Gierke Disease', and 'pyrimidine metabolism' (Source: SMPDB) pathways, and 'interconversion of nucleotide di-and triphosphate' and 'asn N-linked glycosylation' (Source: Reactome) pathways as well as the 'glucagon signaling pathway' (Source: KEGG). The unique down-regulated pathways for spine-435 were: 'gluconeogenesis', Fanconi-Bickel syndrome', and 'oncogenic action of succinate' (Source: SMPDB) pathways, along with 'citrate cycle' and 'aminosugars metabolism' (Source: INOH) pathways, and the 'pyruvate metabolism and TCA cycle' (Source: Reactome) pathway. Common metabolomic and transcriptomic pathways. These pathways are presented in S37-S40 Tables. Similar to S33-S36 Tables, no, or only one or two common up-regulated pathways were found. Of these the single up-regulated pathway: 'arg and pro metabolism' (Source: INOH) found for the brain-435 cell line was shared with the spine-435 cell line while the 'metabolism of amino acids and derivatives' (Source: Reactome) and the 'his, lys, phe, tyr, pro, and trp catabolism' (Source: Reactome) pathways were uniquely up-regulated in lung-435 and spine-435 cell lines respectively. Again, similar to what was found in S33-S36 Tables, several of the top 12-15 down-regulated pathways were shared between cell lines but more of the pathways were unique to each cell line. Thus, for the brain-435 cell line the unique downregulated pathways were: 'pyrimidine metabolism' (Source: Wikipathways) along with  Importantly, a comparison of S33-S36 Tables with corresponding S37-S40 Tables indicated many of the pathways in these sets of tables; e.g., S33 Table vs S37 Table, etc., were the same. Consequently, we found identical pathways, albeit with all but one being down-regulated pathways, that had been discovered by the use of three separate data sets, i.e., proteomic, transcriptomic, and metabolomic. Thus, a comparison of these tables showed a convergence of the three data sets onto common pathways. For the brain-435 cell line these were down-regulated: 'cell cycle', 'cell cycle, mitotic', 'pyrimidine metabolism' (Source: Wikipathways), 'S phase', 'DNA replication', and 'extension of telomeres'; for the liver-435 cell line these were down-regulated: 'pyrimidine metabolism' (Source: EHMN), 'nucleotide di-and triphosphates', 'S

qRT-PCR verification
It has been noted above that proteomic derived protein lists could only be sparingly matched to their corresponding transcripts from the RNAseq analyses. Our qRT-PCR results were similar in that only a small percentage of the protein and transcript levels found during the proteomics and RNAseq analyses could be confirmed with qRT-PCR. Representative genes analyzed by qRT-PCR are given in Table 7 along with their fold changes in their respective protein and transcript (gene) lists. Bar graphs of the qRT-PCR results are presented in Fig 8A and 8B for the isogenic metastatic 231 and 435 cell lines respectively. In all cases, the direction (up-or down-regulated) matched the up-and down-regulation of their proteins and transcripts (Table 13).

In vitro drug testing
From our earlier studies [12,13] and reports from other labs [5,6,8,10,11,32], it is becoming established that metastatic lesions have diverged from their primary tumors at several levels: Table 11. Unique pathways found to be common in both proteomic and transcriptomic analyses for the metastatic lung-231 cell line. genetically, proteomically, and metabolically. Thus, we have shown here that isogenic cell lines derived from different organs have distinct divergences between their proteomes, transcriptomes, and metabolomes that make each of them fundamentally different biological entities. Given this, an important question arises about how cancer cells that have adapted to growth in different organs respond to the same therapeutic regimes. To ascertain whether the same drug will kill each cell line with the same efficacy, we tested 4 drugs on each of the isogenic cell lines: one RK-33 (DDX3X inhibitor) [33][34][35] along with three FDA approved and breast treatment established (doxorubicin (DOX) [36], gemcitabine (GEM) [37], and paclitaxel (PAC) [38] drugs, using an in vitro cell culture assay. Table 14    https://doi.org/10.1371/journal.pone.0242384.g008

PLOS ONE
Divergent pathways of isogenic cell lines alter drug efficacies lung-231 and lymph node-231 (p = 0.005) and trending to significance in comparison to 1t umor-435 (p = 0.056); DOX in the case of brain-435 with an IC 50 value higher-than those of lung-231 (p = 0.009), lymph node-231, and 1˚tumor-435; DOX in the case of lung-435 with an IC 50 value greater-than 1˚tumor-435. In addition, we determined the linear fold change in any differences in efficacy between the 1˚tumors and their isogenic cell lines, which are shown in Table 15 where the values in bold-type and underlined are statistically significant (2-sided Student's t-test p < 0.001 or as indicated with p � 0.05 considered significant).

Discussion
Presently, despite advances in therapies, metastatic breast cancer remains incurable [1,39]. To address the reasons as to why this is the case, several independent laboratories have provided evidence, over the course of several decades that metastatic tumors have, to varying degrees, diverged from their primary tumors [3, 4, 6-8, 11, 32]. Consequently, it is now realized that therapies that are effective against regional breast cancer, which are based on a few molecular markers that have been used to define breast cancer subtypes, can have minimal efficacy against metastases that have diverged from the primary tumor [4,5,7,11]. Accordingly, with the advent of more recent molecular ("omics") characterizations, it is becoming accepted that careful clinical evaluations of the molecular (proteome and transcriptome) as well as metabolism of metastatic lesions would be helpful and likely necessary if efficacious treatment of metastatic disease is to be developed [9,11]. As described in previous reports [12,13] and extended here, our approach to this problem of gaining a better understanding of the molecular and metabolic changes that occur at metastatic lesions relative to their primary tumors, has been the development and characterization of model systems of isogenic cell lines that have been cultured directly from metastatic organ samples. The advantage of this model system with respect to those that rely on several cell lines

PLOS ONE
Divergent pathways of isogenic cell lines alter drug efficacies from different individuals is the isogenic nature of the cell lines. From this perspective, phenotypic, molecular, and metabolic divergences from an isogenic primary tumor, as is the case in clinical settings, can be assessed within the context of the isogenic background of these cell lines. As described above and discussed below, we have found that these cell lines represent unique biological entities that have diverged from their primary tumors in growth characteristics in culture, proteomics, transcriptomics, as well as metabolomics. A principle goal has been to study the sensitivity of these different cell lines to chemotherapeutics as well as in future in vivo studies.
Here we tested a DDX3X (DEAD box helicase) inhibitor (RK-33) [33][34][35] and three clinically established breast cancer chemotherapeutics: gemcitabine (GEM) [37], paclitaxel (PAC) [38], and doxorubicin (DOX) [36]. As shown in Table 14, we found several differences in efficacy across drug treatments of the cell lines. For example, relative to all other cell lines, the lung-231 cell line was the least sensitive to GEM and exhibited a critical~8-fold decrease in sensitivity as compared to its primary tumor cell line. On the other hand, lung-435 and spine-435 were the most sensitive to GEM and again differed (being more sensitive) relative to their primary tumor cell lines. In addition, lung-231 and lymph node-231 cell lines were more sensitive than their primary tumor cell line to DOX but the reverse was the case for the four metastatic 435 cell lines, which were less sensitive to DOX than their primary tumor cell line. Other distinctions were the lung-435 cell line being less sensitive to PAC than was the case for its primary tumor cell line. Although significant differences were observed in the sensitivity to RK-33 the changes were relatively less pronounced, and this latter characteristic of RK-33 may be a therapeutic advantage. That is, as there is a molecular dependency for DDX3X expression in cancer cells to maintain cellular and bioenergetic homeostasis [34,40,41], it is less likely to undergo marked changes during growth and establishment of metastatic tumors. This, in part, could explain why RK-33 doses required to kill the different isogeneic cell lines was the least variable, particularly evident in the 231 cell lines (Table 14), amongst the different chemotherapeutic agents used in this study.
To attempt an explanation for the observed differences in drug efficacies it needs to be noted that the cellular mechanisms involved with GEM, PAC, and DOX are complex and multiple pathways and several proteins must be taken into consideration [36][37][38]. This problem is exemplified from an evaluation of changes in some of the single proteins that might be involved with the decreased efficacy of GEM against lung-231 cells [37]. Thus, there was a 1.7-fold increase in a protein inhibited by GEM (ribonucleotide reductase 1: RRM1) in lung-231 cells, which might be evaluated as a requirement of more drug against this target, i.e., a decrease in sensitivity. However, there was a simultaneous 1.8-fold increase in a solute carrier (SLC29A1) that transports GEM into cells as well as decreases (-1.9-fold in each case) in inactivating enzymes (cytidine deaminase and cytosolic 5' nucleotidase). Consequently, in this case, a simple consideration of proteins involved with GEM metabolism in these cells provides an ambiguous conclusion as to the lack of sensitivity of lung-231 cells to GEM treatment. Contributing to this ambiguity is that in cell lines more sensitive to GEM, such as brain-435 and liver-435 (Table 14), SLC29A1 was found to be decreased by -1.8 and -1.35 respectively, which when coupled with decreases in activating enzymes: deoxycytidine kinase (-1.30-and -1.67-fold in brain-435 and liver-435 respectively), UMP/CMP kinase (-1.25-and -1.35-fold in brain-435 and liver-435 respectively) and nucleoside-diphosphate kinase (-1.50-fold in both cell lines) would tend to lead to the conclusion that GEM ought to relatively ineffective in these cell lines rather than relatively more effective. Similar decreases in activating enzymes were found in the lung-435 and spine-435 cell lines and yet these were the most sensitive to GEM. A similar evaluation of DOX's efficacy in these cell lines also produced similar findings. Thus, a consideration of 18 proteins [36] involved with the transport, export, ability to detox reactive oxidative species (ROS), or repair DNA indicated that in most cases there was minimal or no changes in the levels of these proteins across all cell lines. Exceptions were about a 1.7-fold increase in both a DOX exporting protein (ABCB1) and a DNA repair enzyme (MSH2) in lymph node-231 cells, which were the most sensitive to DOX treatment, i.e., there was not a diminished sensitivity to DOX as compared to any of 435 metastatic cell lines. Importantly, the latter cell lines showed mostly no changes in any of the evaluated proteins except for an increase (~1.35-fold) in two enzymes (SOD1 and CAT) involved with a response to ROS in brain-435 cells and thus, again, little evidence at the protein level for the observed differences in sensitivity to DOX. Finally, exploring changes in tubulins (the major target of PAC) [38] in all cell lines showed no differences that could provide a clue as to the observed differences in sensitivity to PAC. Thus, it appears that differences other than those at the single protein levels, such as at the pathway or pathway network levels, i.e., a combination of proteins and pathways that differ between cell lines brings about the variances in the observed drug sensitivities across these cell lines.
Along these lines, our drug assay results are in line with a recent report that demonstrated a link between differences in protein networks across 41 breast cancer cell lines and changes in the sensitivity of these cell lines to drug treatments [42]. However, such a result may have been expected as the cell lines used were derived from separate individuals as well as being from across all subtypes of breast cancer. Thus, it is well known that each subtype is susceptible to different therapeutics [43] and certainly different individuals often have different responses to chemotherapies, which is always recorded in clinical trial generated patient survival curves. Accordingly, as pointed out above, this was also apparent in a comparison of our isogenic cell lines (triple negative subtype) with a reversal in sensitivity to DOX when the 231 isogenic metastatic cell lines are compared to the 435 derived isogenic metastatic cell lines. More importantly, the report of the efficacy of drugs across 41 cell lines did not take into account, as we have, how metastatic spread to visceral organs likely alters the efficacy of clinical chemotherapies and the authors did not address the treatment of metastatic disease. However, findings of molecular discordances between primary tumors and their metastasis continues to provide evidence that selecting therapeutic regimes that have been based on a characterization of the primary tumor but are aimed at ablating visceral metastatic lesions will likely be ineffective [11]. It is becoming evident that treatment strategies for metastatic disease will likely be more effective if these are based on the fundamental genetic and molecular characteristics of the metastatic lesions. The latter conclusion has been put forth in a recent clinical breast cancer study (reported while our study was in progress) that described an evaluation of evolutionbased mutational changes at metastatic sites that occurred independently from any primary tumor clonal evolution and as such it was suggested that organ-specific microenvironments were driving such changes [11]. Given this evidence the authors suggested that clinical characterization of metastatic lesions ought to be carried out prior to treatment determinations of metastatic cancer [11]. Thus, our hypothesis has been strongly supported by this clinical study. However, in the present study, we have collected transcriptomic, proteomic, and metabolomic data sets and focused our analyses on pathways and their networks, which is a distinct difference from gene/protein mutational analyses that have defined the cited clinical study.
Although mutational evolutionary analyses are providing important insights into the clonal (gene-based) divergences associated with tissue dependent metastatic adaptation(s) (evolution) [9,11,[44][45][46], as previously reported and expanded upon here, our isogenic model system has revealed isogenic cell line specific pathways that indeed have been influenced by the microenvironment of the cell line's organ of origin. Examples have been pointed out in the Results section (Tables 1-6). Thus, this model provides complementary evidence (relative to the clonal evolutionary evidence) that fundamental molecular and metabolic divergences of metastatic tumors from their primary tumor are an unavoidable consequence of growth within a tissue specific microenvironment that differs vastly from that of the breast epithelial microenvironment. Another finding from the pathway analyses approach has indicated that pathways don't necessarily fit a simple binary on/off (up/down) model but instead are likely in a state of homeostasis or steady state of regulation with pathway defining proteins being both up-and down-regulated relative to the primary tumor (S41-S46 Tables). The therapeutic implications of this finding is exemplified by the EGFR1 pathway, which has been considered as a clinical therapeutic target in breast cancer [47] but was found here to not necessarily be in a overexpressed "on" mode across 5 of 6 of the isogenic cell lines studied as this pathway also exhibited down-regulated "off" protein components across all 6 isogenic cell lines (S41-S46 Tables; The down-regulated EGFR1 pathway of lung-231 was not included as these tables show only those pathways found to be simultaneously up-and down-regulated.). In fact, EGFR1 protein was not found to be up-or down-regulated in any of these cell lines with the pathway being defined by several of the other 457 protein members of the pathway, such as up-regulated ASAP1 in liver-435, spine-435, and lymph node-231 cell lines as well as PRKCZ in brain-435, lung-435, spine-435, and lymph node-231 cell lines or the down-regulated APPL2 in lung-435, spine-435, lung-231, and lymph node-231 cell lines as well as ENO1 in all six cell lines. Moreover, there were roughly twice as many down-regulated proteins relative to the up-regulated proteins of this EGFR1 pathway (S41-S46 Tables), which indicates that the identification of a single up-regulated (overexpressed) target of such a complex pathway and using a therapeutic against it may have minimal impact on the pathway, i.e., on treatment. In addition, as exemplified in Figs 4-7, several different pathways can be interconnected into large integrated networks that are likely all regulating each other. Consequently, it is apparent that targeting a single component (protein) in what might be thought of as a single 'key' pathway may be ineffective due to the self-regulation of the pathway or the overlapping function of the interacting network(s). These results provide a partial explanation as to why our analysis of the single proteins involved with the effectiveness of GEM. PAC, or DOX did not show a clear association to the sensitivity of the cell lines to these drugs, i.e., their efficacies are likely based on complex pathway dynamics rather than any single protein.
Finally, it was important to analyze the possible clinical associations that our pathway approach achieved. This was hindered by the fact that available comparative human breast cancer databases (e.g., [48]) report survival as a function of an overexpression (relative to normal tissue) of markers/genes associated with the primary tumor while our studies have been focused on metastases. Nevertheless, we crossed referenced proteins from pathways that were found to be up-regulated (relative to the primary tumors) across multiple isogenic cell lines. Thus, Table 16 shows a randomly selected list of up-regulated pathways and hence proteins found to be common across 2-5 cell lines. It is noteworthy from the clinical data, i.e., elevated expression of these genes at the primary tumor, was associated with both a poor survival (e.g., Fig 9: FLNB [49] and H1F0 [50] genes) as well as an enhanced or better survival (e.g., CDC42 & HLA-A, lower portion of Table 16). Thus, this analysis indicates that an overexpression of proteins at metastatic sites, relative to primary tumor levels as well as, from the clinical data, relative to normal breast tissue levels rather than normal tissue of origin levels can complicate/ contradict the interpretation of disease free survival. That is, the latter implies no metastatic progression and as such the basis of the disease free survival data does not reflect the status of the markers at metastatic sites. Consequently, this lends support to the conclusion that further studies are required to better understand how analyses of pathways at metastatic sites can contribute to a better understanding of the pathology of the metastases as well as of therapeutic options that may enhance survival.

Conclusions
The insights provided by these analyses indicate that the rationale of targeted treatment of metastatic disease may benefit from a consideration that the biology of metastases has diverged from the primary tumor biology and using primary tumor traits as the basis for treatment may not be ideal to design treatment strategies. Thus, exploring an interconnected integrated pathway analysis approach as an alternative to the single gene/protein marker evaluations now in use may provide a better understanding of which pathways are participating in metastatic cancer survival at a specific site. In addition, compiling normal expression levels of markers/pathways specific to different tissues would greatly aid with the discovery of changes in these levels in the metastatic lesions and pave the way for explorations as to how these changes affect treatment outcomes as well as direct future studies aimed at controlling and ablating metastatic disease.

Generation of isogenic metastatic cell lines from specific organs
The  against cell lines at the ATCC bank. Generation and characterization of the parental MDA-MB-435-tdTomato (hence referred to as 435) cell line has been previously described [53]. MDA-MB-231 (hence referred to as 231) cells were not genetically modified and thus were the parental cell line for this line's isogeneic primary tumor and metastatic (lung & lymph node) cell lines. Primary tumors and subsequent tissue specific isogeneic cell lines were generated/cultured as previously described [12]. For the 231 cell lines, generation of growth curves and growth rate analysis was as previously described for the 435 cell lines [12]. All culturing was done in standard humidified incubators at 37˚C and 5% CO 2 . Media were: DMEM-10% FBS for parental cell lines and DMEM:F12 (50:50)-5% FBS for all primary tumor and metastatic cell lines.

Optical microscopy
Phase contrast microscopy was done on a Nikon ECLIPSE TS 100 microscope (Nikon Instruments, Inc.) equipped with a Photometrics CoolSnap ES digital camera (Roper Scientific). Images were collected with NIS-Elements F3.2 software and processed with ImageJ.

Protein preparation
Total protein solutions were prepared by directly lysing cells cultured on 100 mm dishes, which, in all cases, were at about 70-80% confluency. Lysis buffer (200 μl) was: 100 mM Tris pH 6.8, 12% glycerol, and 2% SDS, 1 mM EDTA, and 1:200 dilution of a protease cocktail (Sigma, I1386) (added immediately prior to use). Lysates were placed in 0.5 ml microcentrifuge tubes and sonicated (12-15 bursts) on ice and frozen at -80˚C until use.

Protein concentration estimates
Aliquots (100 μl) from the frozen stocks (thawed on ice) of total protein preparations were placed into 0.5 ml microcentrifuge tubes. Protein concentration estimates were carried out using room-temperature samples diluted (1:10-1:15) in a phosphate-free saline solution (NaHCO 3 (45 mM), NaCl (95 mM), KCl (4.5 mM), CaCl 2 (0.24 mM), MgCl 2 (0.08 mM) pH 7.35). The diluted protein solutions were assayed using a BioRad RC DC kit according the manufacturer's protocol and BSA for the standard curve. This kit was chosen as the protein precipitate formed during step-1 of the protocol is free of compounds that interfere with the step-2 color reagent, such as EDTA, amino acids, lipids, and nucleic acids.

Proteomics
Protein pellets (100 μg each) were submitted to the Mass Spectroscopy and Proteomics Facility at Johns Hopkins University Medical School for routine differential proteomics analyses. The Director: Dr. Robert Cole, oversaw all analyses. State-of-the-art TMTs (tandem mass tags) were applied to digested samples for direct comparison of all 10 samples in a single tandem MS experiment. The mass spectroscopy spectra output was analyzed with Proteome-Discover for peptide identification and, as such, mapped to specific protein identifiers and quantified. Data was further processed to identify fold changes in protein expression levels from isogenic metastatic cell lines relative to their primary tumor cell lines. Briefly, for each sample, the multiple spectra values for each peptide were summed to single values per unique peptide and then the many different peptide values normalized across all the samples to minimize possible technical variation. These quantile normalized log2 values were compared to determine differential peptide expression levels. In addition, all peptides were mapped to their cognate genes, which facilitated annotation and possible downstream functional analyses.

RNAseq
RNA was prepared from frozen cell pellets (-80˚C stored). The frozen stocks were from lots of cell lines that were at the same passage as the stocks that were used for protein preparations or one to two passages later. RNA concentrations and quality control spectrophotometric determinations were done on a NanoDrop microvolume spectrophotometer (ThermoFisher Scientific) and only samples with 260/280 & 260/230 ratios of 2.0-2.1 and 1.8-2.2 respectively were used. RNASeq was performed by a commercial entity (BGI Americas, San Jose, CA). Briefly, total RNA was checked for quality (RIN > 9) and libraries were constructed. Libraries were 50bp single-end sequenced on a BGISEQ-500 instrument to a standard depth of 20 million reads per sample. Sequencing data was filtered and supplied as differential gene expression data sets.

Metabolomics
Metabolomes were generated as previously described [12]: Briefly, metabolite data from all samples were acquired using Agilent 6540 Quadrupole-Time-of-Flight (Q-TOF) mass spectrometer with Agilent 1290 HPLC at the Metabolomics Facility at Johns Hopkins Medical Institution. Data was analyzed using Agilent Mass Hunter and Agilent Mass Profiler Professional (MPP) version 13.1.1 and Agilent Qualitative and Quantitative Analysis Software packages (version 6.00) to determine the metabolic profile of each sample. Aqueous phase metabolites were used in pathway analyses.

Quantitative real-time polymerase chain reaction (qRT-PCR)
RNA was isolated from cell lines and transcribed into cDNA using manufacturer's protocols (Qiagen, Germantown, MD and Bio-Rad, Hercules, CA). Diluted cDNA was used as template for qRT-PCR to amplify target genes in replicates of two on a thermal cycler with primer sequences given in S29 Spreadsheet. Relative change in target gene expression was calculated using the 36B4 gene as housekeeper [54].

Principal component analyses and hierarchical clustering
To reduce the complexity of the large amount of data that was generated from the proteomics (S1-S6 Spreadsheets) and RNA sequencing (S7-S12 Spreadsheets), and to infer relationships between the data sets, we performed principal component analysis (PCA) and hierarchical clustering at both the protein and transcript (gene) level. Hierarchical clustering was performed using Morpheus (https://software.broadinstitute.org/morpheus).

In vitro drug assays
The established FDA approved chemotherapuetic drugs used were: paclitaxel (TSZ CHEM, Cat# RS036, Lot# 061916), doxorubicin (Cayman Chemical, Item# 15007, Lot# NA), and Gemcitabine (Sigma, Cat# G6423-10mg, Lot# 026M4704V). In all cases, cells were plated at 2000 cells per well onto 96 well plates and 24 hrs later treated with each drug over a serial dilution range of drugs: paclitaxel (PAC, 0.01-50 nM), doxorubicin (DOX, 0.01-5 μM), gemcitabine (GEM, 0.001-10 μM), and an in-house DDX3X inhibitor drug RK-33 (1-25 μM). Each concentration of drug was added to cells (wells) in quadruplicate along with no drug added control wells. Two to three biologic replicates were done. Standard colorimetric MTS assays (addition of 10% MTS reagent in medium with a 2 hr incubation) were done 72 hrs after drug treatment. Plots of the spectrophotometric outputs (absorbance vs log[drug]) were used to determine the IC 50 values of each drug.

Statistical methods
As described above, all Proteomics and Metabolomics source datasets were generated at core facilities at the Johns Hopkins University while an outside company generated the RNAseq datasets. As such, we received datasets with completed statistical analyses and all p-and or qvalues presented in all Tables and Spreadsheets were obtained from the source datasets. In the case of Pathway analyses, we utilized the online database, ConsensusPathDB, which has published the statistical methods used [26,27], as described above. For the In vitro drug assay-IC 50 values dataset, we applied F tests to determine unequal or equal variances and then the appropriate two-sided Student's t-tests (p � 0.05) were utilized to evaluate significant differences.