Over the past decade, metabolic reprogramming has been recognized as a hallmark in oncogenic transformation in PDAC and other cancers [17–19]. The well-known aerobic glycolysis or Warburg effect (upregulated glucose uptake and lactate production) has been shown to confer proliferative and survival advantages in multiple cancer types by supplying sufficient bioenergetic precursors and NADPH [20]. PDAC cells also display elevated glutaminolysis to maintain redox balance and scavenge extracellular fatty acids/amino acids to survive in a hostile microenvironment with limited fuel supply [18,21–23]. Heterogeneity exists between metabolic profiles of different cancer types and between different cancer cell populations within the same tumor, owing to the context-specific oncogenic signaling events and micro-environmental factors present [24,25]. The process of EMT involves major changes to the gene expression network and cellular phenotype. It is, therefore, likely to be accompanied by metabolic alterations to accommodate the shift in cell's priority from proliferation to invasion of neighboring tissues and to adapt to changes in the environment. Indeed, a wide range of metabolic alterations have been observed with induction of EMT status in breast, lung, ovarian, cervical and prostate cancers, although the nature of the metabolic reprogramming varies widely across the studies [26–40]. We have previously shown augmentations of glucose consumption and lactate output in Panc1 cells undergoing tumor necrosis factor-α (TNFα)- and transformation growth factor-β (TGFβ)-induced EMT, with differential molecular changes observed in the two models [41].

Alongside changes to cell metabolism in our previous study [41], we also observed the induction of Snail expression concurrently with change in EMT status in Panc1 cell upon treatment with TGFβ or TGFβ combined TNFα. In view of the importance of Snail-dependent EMT in underlying PDAC-related lethality, we sought to induce EMT in Panc1 and the non-tumorigenic human pancreatic ductal epithelial (HPDE) cells [42] via Snail overexpression to investigate the metabolic consequences. Specifically, we chose to compare the effects of Snail overexpression in a pancreatic cell line already on the EMT spectrum (Panc1) to that of a purely epithelial pancreatic cell line (HPDE) to study the specific consequences of Snail induction at different points across the EMT differentiation spectrum. Here we report that EMT in both cell lines is associated with elevated glucose uptake and lactate excretion, as well as downregulation of proteins in the mitochondrial electron transport chain (ETC).

#### **2. Experimental Section**

#### *2.1. Antibodies and Reagents*

Antibodies used are listed as follows: Snail (3879), Vimentin (5741), LDH-A (2012, Cell Signaling Technology, Danvers, MA, USA); E-cadherin (sc-21791), N-cadherin (sc-7939), Beta-actin (sc-47778, Santa Cruz, TX, USA); Hexokinase II (ab37593), Total OXPHOS Human WB Antibody Cocktail (ab110411, Abcam, Cambridge, UK). All other reagents were from Sigma-Aldrich (Sydney, NSW, Australia) unless stated otherwise.

#### *2.2. Cell Culture*

Panc1 cells from ATCC were cultured in Dulbecco's modified Eagle's medium (DMEM) containing 4.5 g/L glucose supplemented with 10% fetal bovine serum and penicillin-streptomycin. The human pancreatic ductal epithelial (HPDE) cell line [42] were a kind gift from Dr. Phoebe Phillips at Lowy Cancer Institute, UNSW Australia. HPDE cells were cultured in keratinocyte serum-free (KSF) media (ThermoFisher Scientific, Waltham, MA, USA) containing 1.6 g/L glucose supplemented with epidermal growth factor (5 ng/mL), bovine pituitary extract (50 ug/mL) and penicillin-streptomycin.

#### *2.3. Overexpression of the SNAI1 Gene*

Plasmid containing human SNAI1 (encoding Snail) cDNA (Addgene plasmid 23347) was a kind gift from Bob Weinberg [43]. Stable SNAI1-expressing Panc1 cell line was generated using retroviral-mediated pBabe-puro-Snail infection. Briefly, 8 × 105 Hek293-FT cells were seeded on 10 cm culture dishes and allowed to attach overnight before being transiently transfected with gag/pol and VSV-G packaging plasmids, along with the pBabe-puro-Snail plasmid or empty pBabe-puro vector. The culture media for Hek293-FT was refreshed 12 hours later and media containing viral particles were harvested 24 hours and 48 hours after media refreshment. The 24 hours and 48 hours media was combined and applied to target cells at 40% confluence with polybrene (8 μg/μL) for 24 hours. A pooled cell population was used for experiments following puromycin selection.

#### *2.4. Cell Morphology*

Pictures of cells were taken using a phase contrast microscope (Nikon Eclipse TS100, Nikon, Tokyo, Japan) attached to a camera (Nikon digital sight) under 40× or 100× magnification.

#### *2.5. Western Blotting*

Cells (1 × 106) growing on 6-well plates were lysed in RIPA buffer containing protease inhibitors as described previously [44]. After denaturation, samples (20 ug protein) were resolved by 10% polyacrylamide gel electrophoresis (187 V, 1 h) and transferred to PVDF membranes (65 V, 65 min). For immunoblotting, membranes were blocked in 5% skim milk in tris-buffered saline containing Tween 20 (TBST), incubated with primary antibodies overnight at 4◦C, washed with TBST and incubated with secondary antibodies in 5% skim milk in TBST for an hour. After washing, membranes were developed with enhanced chemiluminescence reagents (Western Lighting Plus-ECL, Perkin Elmer, Waltham, MA, USA) and visualized under Las4000 imager (GE Healthcare, Chicage, IL, USA). Densitometry was performed using ImageJ software by obtaining the optical density of each band.

#### *2.6. Quantitative PCR*

Total RNA was extracted using the Roche High Pure RNA Extraction kit according to manufacturer's instructions (Roche, Basel, Switzerland). RNA (1 μg) was then reverse-transcribed using the Roche Transcriptor first strand cDNA synthesis kit. The resulting cDNA was mixed with primers (primer sequences sourced from Sigma KiCqstart or the Primer Bank [45]) and SYBR-Green (Roche, Basel, Switzerland) in 96-well plates. Quantitative PCR was performed using Roche 480 Light Cycler to obtain Ct values for each gene of interest and a housekeeper beta-actin. Analysis was conducted using the ΔΔCT method.

#### *2.7. Wound Healing Assay*

Cells were grown in 6-well plates and scratch wounds were created by scraping confluent cell monolayers with a sterile pipette tip on 3 sites on each well. The cells were then incubated under normal conditions with refreshed media under the Nikon Tie inverted time-lapse microscope for 24 hours. Migration at 24 hours was quantified by measuring the area closed between two moving borders of the cells from each scratch. Values from the 3 wounds on each of the triplicate wells were averaged and 3 independent experiments were carried out.

#### *2.8. Measurement of Cell Proliferation*

Cells were seeded in 6-well plates at a density of 1.5x105 per well. After 4 days, cells were rinsed with PBS and immersed in 0.5% crystal violet (w/v)/50% methanol (v/v) solution and left to fix for 20 min. After fixation, cells were gently rinsed to remove all the crystal violet solution and allowed to dry overnight. The next day, fixed cells were solubilized with 1% SDS at 37◦C and 50 μL of the solution was taken for absorbance measurements at 570 nm using a plate reader as an indication of cell number of the well.

#### *2.9. Glucose Uptake Assay*

Glucose uptake was assessed using the glucose analogue 2-deoxyglucose (2-DG). Cytochalasin B (25 μM) was applied to control wells for 15 min before the assay and during the assay to give a measure of background glucose uptake.

After washing with PBS, cells in 6-well plates were incubated in Ringer solution (140 mM NaCl, 20 mM HEPES, 5 mM KCl, 2.5 mM MgSO4, 1.2 mM CaCl2, pH7.4) containing 10 μM 2-DG and 0.5 μCi/mL radio-labelled 3H-2-DG for exactly 8 min. Following incubation, cells were washed with cold PBS and lysed in 1 M NaOH. The amount of 3H radioactivity in lysates was counted using a beta-counter (Tri-Carb liquid scintillation counter, Perkin Elmer, Waltham, MA, USA) from which background was subtracted. Protein concentrations of lysates were measured using BCA assay (Pierce BCA protein assay kit, Thermo Fisher Scientific, Waltham, MA, USA) for normalization.

#### *2.10. Lactate Assay*

Lactate concentrations in cell culture media (72 hour after plating) were determined in a reaction mixture containing hydrazine hydrate (0.4 M, pH 9.0), EDTA (10 mM, pH 9.0) and NAD+ (0.5 mM). Samples and standards were added into 96-well plates followed by lactate dehydrogenase (10 units/well) and the amount of lactate was assessed by measuring the amount of NADH formed at 340 nm after 2 hours of incubation at 37◦C (a timepoint when lactate conversion is complete).

#### *2.11. Bioenergetic Profiling of Oxygen Consumption and Extracellular Acidification in Snail Overexpressing and Control Cells*

Cells were seeded at 2 × 104 per well on a XF96 seahorse cell plate (Agilent Technologies, Santa Clara, CA, USA) in their respective growth media. The next day, cells were washed with Seahorse assay media containing 25 mM glucose, 2 mM L-glutamine, 1 mM sodium pyruvate, pH 7.4 equilibrated in the same media at 37 ◦C for 30 min in a CO2-free incubator. Bioenergetic profiling was performed by monitoring oxygen consumption and extracellular acidification rates at basal levels, followed by sequential injections of 1 μM oligomycin (an ATP synthase inhibitor), 0.5 μM Carbonyl cyanide-4-(trifluoromethoxy)-phenylhydrazone (FCCP) (a mitochondrial uncoupler) and 1 μM rotenone (a Complex I inhibitor) using the Seahorse XF96 Analyzer. The time-course of energetic profiles, as well as basal oxygen consumption, basal extracellular acidification rate and maximal oxygen consumption were calculated from the primary data.

#### *2.12. Glucose Oxidation*

Glucose oxidation was measured in cells seeded in 6-well plates (1 × 106 cells per well). Briefly, cells were washed with PBS and incubated in DMEM containing 1 g/L D-glucose and 2 μCi/ml 14C-glucose for 1 hour at 37 ◦C. After incubation, the culture media was added to 1 M perchloric acid and the CO2 released was absorbed in 1 M NaOH solution over 2 hours. The CO2 produced was quantified by counting 14C content in the NaOH solution using a beta-counter (Tri-Carb liquid scintillation counter, Perkin Elmer, Waltham, MA, USA).

#### *2.13. Measurement of Half-Maximal Inhibitory Concentration*

The half-maximal inhibitory concentration (IC50) was measured by crystal violet assay. Cells were seeded at 1 × 104 cells per well into a 96-well plate in triplicate and allowed to attach overnight. The next day cells were treated either with serial dilutions of gemcitabine or paclitaxel. The IC50 was measured after 48 hours by crystal violet assay as described above (*n* = 3 biological replicates), with cell viability being expressed relative to vehicle control (phosphate buffered saline for gemcitabine, 0.1% ethanol for paclitaxel). The IC50 was then calculated by non-linear regression by fitting the log-transformed drug concentration against relative cell viability.

For comparison under different glucose conditions, cells were allowed to adhere overnight in high glucose DMEM (i.e., 4.5 g/L glucose) before being treated with serial dilutions of gemcitabine spiked with an IC75 dose of paclitaxel in media containing either high or no glucose.

#### *2.14.13C metabolic Tracer Experiment and Metabolomics*

Triplicates of Panc1 and HPDE cells were cultured in 6-well plates in their respective glucose-free DMEM and KSF media as described earlier. Approximately 4.5 g/L and 2.9 g/L of uniformly labelled 13C6-glucose was added to DMEM and KSF media respectively and cells were cultured for 5 hours. To measure the accumulation and 13C enrichment of extracellular pyruvate and lactate, 50 μL culture media was harvested hourly. The collected media were centrifuged (300 *x g*, 4◦C) for 5 min and the supernatant stored at −30 ◦C until analysis by gas chromatography mass spectrometry (GCMS) using an extraction and derivatization described previously [46]. To measure 13C enrichment of intracellular metabolites, cells were quenched at the end of the 5-hour culture, and metabolites were then extracted and derivatized for GCMS analysis [41]. GCMS of derivatized metabolites was conducted using g a HP-5ms capillary column (0.25 mm i.d. × 30 m × 0.25 μm; Agilent J&W, Agilent Technologies, Santa Clara, CA, USA) installed in an Agilent HP 6890-5973 gas chromatography/mass selective detector.

#### *2.14. 13C Flux Analysis*

Flux modelling was performed to explain the activity of catabolic pathways used by Panc1 and HPDE to metabolize glucose. Metabolic fluxes, which is a measure of metabolite flows, can be estimated by quantitatively fitting a metabolic model to the metabolite data [46]. A simple 13C metabolic flux analysis model was used, comprising of glycolysis, pentose-phosphate pathway and TCA cycle [47]. The metabolite data used for the fit included the isotopomer abundances of extracellular pyruvate and lactate, and the enrichment fractions of intracellular pyruvate, lactate, malate, 2-oxoglutarate, citrate, succinate, alanine and aspartate measured by GCMS [48]. Fluxes were then estimated by least-square optimization such that simulated results gave the best fit to the experimental data. Due to the lack of absolute abundance data for intracellular metabolites, metabolite data were simulated under the assumptions of both metabolic (i.e., constant fluxes) and isotopic steady-state (i.e., maximum intracellular 13C enrichment) [46]. Flux changes due to Snail overexpression was quantified by Monte-Carlo analysis [46]. In this bootstrapping approach, the dataset was repeatedly corrupted with Gaussian noise 200 times, and fluxes were re-estimated each time. The resulting flux distributions were

then used to quantify flux changes. Full 13C flux analysis results are provided within the Supplementary Material, Table S1.

#### *2.15. Statistical Analysis*

Unless indicated otherwise, results comparing vector infected and Snail over-expressing Panc1 or HPDE cells were analyzed by student t-test and expressed as means ± standard error of the mean (SEM). Statistical significance was set at *p* < 0.05.

#### **3. Results**

#### *3.1. Comparison of Basal Levels of EMT Markers in Panc1 and HPDE Cells Establishes EMT Status in Panc1 Cells*

Prior to generation of Snail overexpressing Panc1 and HPDE cell lines, we first sought to determine their basal levels of EMT status. To achieve this, we performed immunoblotting on both Panc1 and HDPE cells cultured under normal conditions to look at basal markers of EMT status, including E-cadherin, N-cadherin, and vimentin (Figure 1). These preliminary immunoblotting experiments confirmed that Panc1 cells are natively somewhere along the EMT spectrum, displaying both markers of epithelial cell type (E-cadherin) as well as markers of mesenchymal status. Conversely, HPDE cells only displayed markers of epithelial status, indicating little to no induction of EMT.

**Figure 1.** Immunoblotting of basal levels of EMT markers E-cadherin (E-cad), N-cadherin (N-cad), and vimentin in Panc1 and HPDE cells. β-actin was used as loading control.

#### *3.2. Snail Overexpression Induced EMT in Panc1 and HPDE Cells*

To study the metabolic changes associated pancreatic cells either already on the EMT spectrum or pancreatic cells with little EMT induction, we overexpressed the principal EMT-inducing transcription factor Snail in the PDAC cell line Panc1 and in non-tumorigenic HPDE cells respectively. Cells were infected with either the empty retroviral pBabe-puro vector (vector) or vector containing human SNAI1 (Snail). Two weeks after puromycin selection, surviving cells of the Snail clones in both cell lines displayed distinct morphology compared to the vector control in that they were more spindle like and dispersed, suggesting the dissociation of tight junctions (Figure 2A or Figure 2E). In Panc1, the increase in Snail (15-fold, *p* < 0.01) was coupled with marked reductions of E-cadherin levels (*p* < 0.001) in Snail-overexpressed cells, while levels of mesenchymal markers (N-cadherin and vimentin) presented little change (Figure 2B). In HPDE cells, N-cadherin and vimentin, as well as Snail, were only present at negligible levels in vector control but were remarkably induced upon Snail overexpression (80-fold increase, Figure 2F). The overexpression of Snail in HPDE also resulted in significant decreases in E-cadherin levels (Figure 2F).

To assess the functional effect of Snail overexpression in terms of migratory capacity, vector control and Snail-overexpressed cells were subjected to wound healing assays. Migration as indicated by the area of wound closure over 24 hours surprisingly did not differ between vector and Snail Panc1 cells (Figure 2C), while Snail resulted in increased percentage of wound closure over 24 hours in HPDE cells (*p* < 0.05, Figure 2G). The proliferation of cells over a 4-day period, measured by crystal violet assays, was slightly but significantly (*p* < 0.01) slowed down by Snail overexpression in Panc1 but not in HPDE cells (Figure 2D or Figure 2H).

**Figure 2.** Snail overexpression induced EMT in Panc1 (**A**–**D**) and HPDE (**E**–**H**) cells. Vector control (V) and Snail-overexpressing (S) cells were generated in Panc1 via retroviral-mediated infections. (**A**,**E**) Representative cell images were taken under bright field microscopy. (**B**,**F**) Cell lysates were resolved by SDS-PAGE and immunoblotted with anti-E-cad, anti-N-cad, anti-vimentin, and anti-Snail antibodies with β-actin used as a loading control. (**C**,**G**) Cell migration as measured by wound healing assay. (**D**,**H**) Cell proliferation as measured by crystal violet assay. Results are shown as mean ± SEM with *n* = 3. \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001 for difference between vector control and Snail-overexpressing cells.

#### *3.3. Snail Overexpression Resulted in Increased Glucose Uptake and Lactate Secretion in Panc1 Cells*

The upregulation of aerobic glycolysis, or increased glycolysis and lactate production in the presence of sufficient oxygen, has been frequently observed during tumorigenesis and in some cancer cells undergoing EMT [26,27,32,34–36]. It was also one of the most pronounced changes seen with TGFβ-induced EMT in Panc1 cells, during which Snail was induced [41]. Here the overexpression of Snail in Panc1 cells also resulted in the upregulation of glucose uptake by nearly 2-fold (*p* < 0.01, Figure 3A). There were no associated changes in SLC2A1 (encoding Glut1), SLC2A3 (encoding Glut3) or HK2, the enzyme converting glucose to glucose-6-phosphate (Figure 3D–F). Snail overexpression in Panc1 cells also caused a 2-fold increase in secreted lactate over a 3-day period (*p* < 0.01, Figure 3B). Using the Seahorse XF96 Analyzer system, lactate production was measured in a more acute setting as the basal extracellular acidification rate sampled over 20 min (Figure 3C). The basal extracellular acidification rate (ECAR; i.e., first 20 mins of assay) was slightly higher in Snail overexpressing Panc1 cells (*p* < 0.05, Figure 3C,D), while no difference was reported in maximal glycolysis after inhibition of mitochondrial ATP production by oligomycin. This lack of difference maximal glycolysis meant that Snail-overexpressing Panc1 cells had lower glycolytic capacity (i.e., the difference between maximal and basal glycolysis; *p* < 0.05; Figure 3C,D). Despite the changes in lactate production, no alterations were seen in LDH-A levels or LDH-B, MCT1 expressions whereas MCT4 expression was enhanced slightly (*p* < 0.05) (Figure 3F–I).

**Figure 3.** Effects of Snail overexpression on glucose uptake and lactate production in Panc1 cells. (**A**) Rate of glucose uptake measured using by 3H-2-Deoxy-Glucose uptake over an 8 min period. (**B**) Lactate assay performed on cell culture media after 72 hours of culture. (**C**) Extracellular acidification rate (ECAR) measured using the Seahorse XF96 Analyzer at basal levels for 30 min followed by sequential injections of 1 μM oligomycin, 0.5 μM Carbonyl cyanide-4-(trifluoromethoxy)-phenylhydrazone (FCCP) and 1 μM rotenone at 30 min intervals. (**D**) basal and maximal glycolytic activity and glycolytic capacity calculated from Seahorse data, (**E**) Fold-change in total RNA for SLC2A1 (encoding Glut1) and SLC2A3 (encoding Glut3) measured by qPCR using β-actin as housekeeper, (**F**) Immunoblotting results for hexokinase II (HK2) and lactate dehydrogenase A (LDH-A) with β-actin as loading control. (**G**–**I**) Fold change in total RNA detected for lactate dehydrogenase-B (LDH-B) and monocarboxylate transporter 1 (MCT1) and 4 (MCT4). Results are shown as mean ± SEM with *n* = 3. \* *p* < 0.05, \*\* *p* < 0.01 for difference between vector control and Snail-overexpressing cells.

#### *3.4. Snail Overexpression Resulted in Increased Glucose Uptake and Lactate Production in HPDE Cells*

HPDE cells overexpressing Snail displayed elevated glucose uptake to a level comparable to Panc1 cells (2-fold, *p* < 0.001, Figure 4A). The augmentation of glucose uptake was accompanied by reduced SLC2A1 (Glut1) expression (*p* < 0.001, Figure 4D) and a 4-fold increase in SLC2A3 (Glut3) expression (*p* < 0.01, Figure 4E). The level of HK2 did not differ between vector and Snail clones (Figure 4F). Although levels of lactate production in vector control HPDE cells were higher than that observed for vector control Panc1 cells, Snail-overexpressing HPDE cells also showed an increase in lactate accumulation in the culture media (Figure 4B). This increase in lactate output was also apparent over the 30-minute period in which basal ECAR was measured using the Seahorse XF96 Analyzer (Figure 4C,D). Similar to Snail overexpressing Panc1 cells, there was no difference observed in maximal glycolysis rate leading to an overall decrease in glycolytic capacity in Snail overexpressing HPDE cells (*p* < 0.05). Among the lactate production (LDH-A and B) and secretion (MCT1 and 4) mRNA measured, only LDH-B transcript displayed a nearly 2-fold increase (*p* < 0.001) (Figure 4F–I).

**Figure 4.** Effects of Snail overexpression on glucose uptake and lactate production in HPDE cells. (**A**) Rate of glucose uptake measured by 3H-2-Deoxy-Glucose tracer over an 8 min period. (**B**) Lactate assay performed on cell culture media collected after 72 hours of culturing. (**C**) Extracellular acidification rate (ECAR) measured using the Seahorse XF96 Analyzer at basal levels for 30 min followed by sequential injections of 1 μM oligomycin, 0.5 μM Carbonyl cyanide-4-(trifluoromethoxy)-phenylhydrazone (FCCP) and 1 μM rotenone at 30 min intervals. (**D**) basal and maximal glycolytic activity and glycolytic capacity calculated from Seahorse data. (**E**) Fold-change in total RNA for SLC2A1 (encoding Glut1) and SLC2A3 (encoding Glut3) measured by qPCR using β-actin as housekeeper. (**F**) Immunoblotting results for hexokinase II (HK2) and lactate dehydrogenase A (LDH-A) with β-actin as loading control. (**G**–**I**) Fold change in total RNA detected for lactate dehydrogenase-B (LDH-B) and monocarboxylate transporter 1 (MCT1) and 4 (MCT4). Results are shown as mean ± SEM with *n* = 3. \*\* *p* < 0.01, \*\*\* *p* < 0.001 for difference between vector control and Snail-overexpressing cells.

#### *3.5. Snail Overexpression Impacted on Oxidative Metabolism in Both Panc1 and HPDE Cells*

Following the observations of enhanced aerobic glycolysis, we next investigated the overall and glucose-specific oxidative metabolism in Panc1 and HPDE cells. Using the Seahorse XF96 Analyzer, OCR at basal levels were not different in either the Panc1 or HPDE Snail overexpressing cells compared with their respective vector controls. Maximal OCR (after addition of the mitochondrial uncoupler FCCP that elicits maximal respiration) was significantly reduced upon Snail-overexpression-induced EMT in both Panc1 (*p* < 0.001) and HPDE (*p* < 0.05) cells (Figure 5B or Figure 5E). Notably, the maximally stimulated OCR in Snail-overexpressed Panc1 cells was nearly halved in comparison to vector control and was not higher than its basal level (Figure 5A). Alongside decreases in maximal respiration in Snail overexpressing Panc1 cells (*p* < 0.001), significant decreases were also observed in ATP production (*p* < 0.05), spare respiratory capacity (*p* < 0.001), and non-mitochondrial respiration (*p* < 0.01, Figure 5B). These observations were in line with significant (*p* < 0.05) decreases in the content of mitochondrial ETC subunits V, III, II and I in Snail-overexpressing Panc1 cells (Figure 5C). Glucose-specific oxidation measured using 14C-labelled glucose tracers was, however, not altered in Panc1 (Figure 5D).

**Figure 5.** Effects of Snail overexpression on oxidative metabolism in Panc1 and HPDE cells. (**A**,**E**) Oxygen consumption rate (OCR) was measured using the Seahorse XF96 Analyzer at basal levels for 30 min followed by sequential injections of 1 μM oligomycin, 0.5 μM Carbonyl cyanide-4-(trifluoromethoxy)-phenylhydrazone (FCCP) and 1 μM rotenone at 30 min intervals. Basal and maximal (following FCCP injection) values of OCR from the time-course data were used for statistical analysis. (**B**,**F**) Basal respiration, ATP production, maximal respiration, spare respiratory capacity, proton leak and non-mitochondrial respiration calculated from Seahorse trace data for Panc1 and HPDE cells respectively. (**C**,**G**) Immunoblotting results for ETC complex I, II, III, V antibodies, with β-actin used as loading control. Densitometry on western blots was performed using image J. (**D**,**H**) Glucose oxidation was measured using the U-14C-Glucose tracer over one-hour period. Results are shown as mean ± SEM with *n* = 3. \* *p* < 0.05, \*\* *p* < 0.01 for difference between vector control and Snail-overexpressing cells.

In contrast to Panc1, the OCR at basal level and after addition of oligomycin, FCCP and rotenone remained unaffected by Snail overexpression in HPDE cells (Figure 5E). Snail overexpressing HPDE cells displayed a lower maximal respiration compared with vector control (*p* < 0.05), and a trend towards decreased spare respiratory capacity (Figure 5F). The decrease observed in maximal respiration observation was in line with reductions in the levels of ETC complexes in HPDE cells upon Snail-driven EMT (*p* = 0.088 for complex III; *p* < 0.01 for complex I), albeit to a lesser extent than Panc1 (Figure 5G). Glucose specific oxidative activity was unaltered in Snail overexpressed cells, as indicated by the 14CO2 produced from 14C-glucose substrate (Figure 5H).

#### *3.6. 13C Flux Analysis Validated Observed Changes in Aerobic Glycolysis and TCA Cycle Activity*

A modelling approach was used to provide a coherent interpretation of metabolic fluxes using metabolite data obtained for Panc1 and HPDE cells with or without Snail overexpression (Figure 6A). Fluxes were estimated by fitting the metabolic model to the measured accumulation rate of extracellular

lactate and pyruvate, and to the 13C enrichment pattern of intracellular metabolites (Figure 6B,C). Despite similarities in the enrichment fractions (Figure 6B), our analyses accounted for the fact that media used between the two cell lines (i.e., DMEM and KSF media) were different. In DMEM, glucose was 95% labelled and lactate was present, whereas in KSF media glucose was only 59% labelled and had no lactate (Figure 6A or Figure 6C). Similarly, only DMEM contained free alanine (results not shown), which explained the significant dilution of 13C-enriched intracellular alanine in Panc1 compared to HPDE (Figure 6B). The assimilation of unlabeled pyruvate was added ad-hoc to better fit metabolite data from HPDE cells.

**Figure 6.** Metabolomics and flux analysis of the effects of Snail overexpression on glucose metabolism in Panc1 and HPDE cells. (**A**) 13C metabolic flux analysis model used to interpret metabolite data. Input substrates shown as gray. Average glucose 13C-enrichment in DMEM and KSF media shown to highlight differences in starting label. (**B**) 13C-enrichment of intracellular metabolites, showing measured fractional enrichments (left bars) and simulated results from the best fit (right bars). Data are shown with mean ± 1 SD with *n* = 3. (**C**) Abundance and 13C-enrichment trajectories of extracellular lactate and pyruvate. "m3" represents fraction of lactate and pyruvate labelled at all three carbons. Results from the best fit are shown as red dotted lines. Error bars represent ± 1 SD. (**D**) Box-and-whiskers plots show estimated fluxes. PDH, PC and MEM fluxes normalized to glutamine uptake flux. MEC flux expressed as a fraction of total ME flux. All flux changes were significant, unless indicated otherwise. Enzyme: MEC (cytoplasmic malic enzyme), MEM (mitochondrial malic enzyme), PC (pyruvate carboxylase), GLS (glutaminase), glycolysis (pyruvate kinase), PDH (pyruvate dehydrogenase). Metabolites: PYR (pyruvate), LAC (lactate), ALA (alanine), AKG (2-oxoglutarate), MAL (malate), ASP (aspartate), CIT (citrate), SUC (succinate), R5P (ribose 5 phosphate), AcCoA (acetyl-CoA), GLN (glutamine). NS *p* > 0.05.

Flux results confirmed that aerobic glycolysis increased in Snail overexpression, for both Panc1 and HPDE cells (Figure 6D). Basal lactate secretion rates were lower in Panc1 than HPDE cells; this is consistent with the results from the enzymatic lactate assay and ECAR (Figure 3B–D and Figure 4B–D). Likewise, the reduction of basal OCR in Panc1 but not HPDE (Figure 5A,E) was reproduced by the estimated pyruvate dehydrogenase fluxes (Figure 6D). Only Panc1 cells with Snail overexpression showed a net consumption of extracellular pyruvate (Figure 6C), although all cultures were secreting pyruvate from glucose.

Flux modelling revealed a few features not immediately observed from metabolite data. The overexpression of Snail in Panc1 cells reduced TCA cycle activity to a greater extent than HPDE cells. This was achieved by lowering both pyruvate dehydrogenase flux and complete oxidation of glutamine, with the latter indicated by a reduced mitochondrial malic enzyme activity and a concomitant shift towards cytoplasmic malic enzyme flux (Figure 6D). This metabolic configuration was required to reproduce the reduced enrichment of intracellular pyruvate and lactate, which could not be solely accomplished by the reversible exchange between the respective intracellular and extracellular pools. Overall, 13C flux analysis showed that Snail overexpression increased aerobic glycolysis and altered carbon flow in the TCA cycle, more so in Panc1 cells than in HPDE.

#### *3.7. Snail Overexpression in Panc1 Cells Does Not Increase Resistance to Gemcitabine or Combination Gemcitabine-Paclitaxel Therapy*

Given the known role of EMT in chemoresistance in PDAC [10], we sought to determine if Snail-overexpression and subsequent EMT induction in Panc1 cells could increase resistance to chemotherapies commonly used in PDAC treatment. To this end, we determined the half-maximal inhibitory concentration (IC50) of gemcitabine, paclitaxel and gemcitabine combined with an IC75 dose of paclitaxel (Table 1). Overexpression of Snail in Panc1 cells did not result in an increased resistance to gemcitabine alone or in combination with an IC75 dose of paclitaxel, but Snail-overexpressing Panc1 cells were slightly more sensitive to paclitaxel monotherapy. We hypothesized that Snail overexpression may give PDAC cells enhanced chemoresistance under the low nutrient microenvironment conditions typically found in PDAC tumors as a result of the dense stromal/desmoplastic reaction [49], but no statistically significance differences were observed between the IC50 of Snail-overexpressing and vector controls treated with combination therapy under limiting glucose conditions.

**Table 1.** The half-maximal inhibitory concentration (IC50) of gemcitabine, paclitaxel and combination gemcitabine with an IC75 dose of paclitaxel in vector control and Snail overexpressing Panc1 cells. Combination treatment was given under both high and no glucose media conditions.

IC50 (M) with 95% confidence interval. Vector compared with Snail overexpression: \* *p* < 0.05.

#### **4. Discussion**

The occurrence of EMT in response to micro-environmental factors partially underlines the malignant phenotype and chemoresistance of PDAC. High levels of Snail, a potent EMT-inducing transcription factor, closely correlate with lymph node invasion and distant metastasis in human PDAC samples [10–12,50–52]. Attenuation of Snail expression in PDAC cell lines resulted in the reversal of EMT, together with decreased sphere and colony formation capacity [16]. It has been shown in several studies that PDAC cells overexpressing Snail underwent EMT and exhibited EMT-associated invasive behaviors both in vitro and in vivo [10,13–15]. In the present study, stable overexpression of Snail in the PDAC cell line Panc1 and non-tumorigenic HPDE cells resulted in pronounced EMT-like phenotypic change as evidenced by alterations in morphology, epithelial/mesenchymal markers and, in the case of HPDE, enhanced migratory capacity (Figure 2). There were also several adaptations in glucose and oxidative metabolism observed with snail overexpression in both cell lines. Despite these alterations to metabolism, Snail-overexpression in PDAC cells did not result in enhanced resistance to gemcitabine or combination gemcitabine/paclitaxel therapy when cultured under either high or limited glucose conditions (Table 1).

Over the past decade, metabolic reprogramming has been increasingly recognized as a hallmark of oncogenic transformation in PDAC and other cancers [17–19]. A small body of literature has also emerged in the last 5 years uncovering additional metabolic alterations related to EMT in breast, lung, ovarian, cervical and prostate cancers, but the actual changes vary considerably across different EMT models and cancer types [26–40]. The EMT events were mostly accompanied by elevated levels of more than one EMT-inducing transcription factors, with Snail being a principal player in the majority of cases [28–33,36]. In the context of PDAC, Snail was highly induced in Panc1 during TGFβ induced EMT, which was associated with upregulated aerobic glycolysis [41]. The induction of EMT by Snail overexpression in the current study was also accompanied by augmented glucose uptake, lactate production and increased levels of basal glycolytic activity, with changes in the expression of transporters and enzymes involved in these processes, namely increased MCT4 expression in Panc1 and higher GLUT3 as well as LDH-B expression in HPDE cells (Figures 3 and 4). In addition, marked downregulations of mitochondrial ETC subunits content and, particularly in the case of Panc1, impaired overall oxidative metabolism were evident.

A direct role of Snail in regulating the glycolytic process has been reported in several occasions. Dong et al. (2013) observed an inverse correlation between levels of Snail and the gluconeogenic enzyme FBP1 in breast cancers [27], hence favoring glucose flux through glycolysis rather than the reverse direction. In basal-like breast cancer, an aggressive subtype containing abundant EMT features, the Snail-G9a-Mnmt1 complex was shown to directly bind to the FBP1 promoter, leading to DNA methylation and transcriptional silencing of the gene [27]. The Snail-mediated suppression of FBP1 was thought to promote glucose uptake and lactate production via improving insulin sensitivity and decreasing PDH (the enzyme catalyzing the conversion of pyruvate to acetyl-CoA and therefore mitochondrial oxidation) activity [27]. In the same vein, Snail levels were high in the more aggressive and castration-resistant subtype of prostate cancer, in which Snail depletion reduced glucose consumption and lactate production [38]. Snail was found to regulate metabolism through miRNA-126-mediated RPS6KB1/HIF1α/PKM2 signaling [38]. In Madin Darby Canine Kidney (MDCK) cells, Snail overexpression resulted in increased activity of phosphofructokinase, a rate-limiting glycolytic enzyme promoting the opposite process to FBP1 [37]. Increased PDK1 expression and the consequent reduction of PDH activity also pointed to the diversion of glycolytic flux towards lactate synthesis [37]. The glycolytic switch was also observed in the breast cancer cell lines MCF-7 and MDA-MB-231 when Snail was induced by Wnt signaling or E-cadherin knock-down [39].

The augmentations of glucose uptake, lactate production, and basal glycolytic activity seen in both Panc1 and HPDE, despite differences in molecular changes, could contribute to EMT-related functional properties especially increased migratory and invasive potential. Enhanced aerobic glycolysis has been associated with invasive cancers and several glycolytic enzymes have been shown to stimulate migration via signaling effects [53–55]. The strongest argument in favor of the glycolytic dependency of migration came from observations that mesenchymal prostate and breast cancer cells exhibited higher aerobic glycolysis, cytoskeletal remodeling and faster migration than epithelial counterparts while no difference in mitochondrial ATP production was found [56]. Migration was attenuated only by inhibition of glycolysis but not mitochondrial respiration [56]. As cell migration is an energy-expensive process involving major remodeling of the cytoskeletal network [57,58], one benefit of the profound upregulation of the glycolytic pathway in PDAC EMT models is presumably to maintain a steady and rapid supply of ATP for cellular migration. The increased lactate secretion by tumor cells could result

in an acidic peri-tumor microenvironment which induces MMP-9 expression and the release of other proteolytic enzymes to degrade components of the ECM [59,60]. This is of particular importance in PDAC, which exhibits a prominent desmoplastic reaction involving extensive proliferation of stromal cells and ECM deposition, constituting a physical barrier for tumor cell extravasation [49]. The flow of H+ along its concentration gradient to adjacent normal tissues could also lead to toxic effects in normal cells such as stromal cells but not cancer cells that developed resistance to low pH environment during carcinogenesis [61]. In addition, lactate has been reported to directly enhance tumor cell motility and contribute to tumor immune escape by inhibiting monocyte migration and cytokine release [62]. However, overexpression of Snail within a metastatic subclone of PC-3 prostate cancer cells reduced both glucose and lactate consumption and increased oxidative metabolism, indicating that expression of EMT features may not always coincide with a higher glycolytic phenotype across all cancer types [63].

There was evidence of downregulation of ETC complex subunits in both cell lines with Snail overexpression. These changes were more pronounced in Panc1 cells compared to HPDE and translated into a functional deficit where both the maximal OCR and the flow of carbons into the TCA cycle decreased in Panc1 cells. The inhibition of oxidative phosphorylation by Snail was implicated in Dong et al. (2013) where Snail-mediated reduction in FBP1 resulted in the loss of mitochondrial transcription factor B1M (TFB1M), leading to defects of protein translation in ETC complex I components [27]. As complex I and III are the main sites of ROS production, the downregulation of complex I level caused by Snail was accompanied by ROS reduction [27]. Lee et al. (2012) showed the direct binding of Snail to promoters of three Cytochrome c oxidase (COX) subunits of ETC complex IV [39]. Complex IV activity and mitochondrial respiration were impaired as a result but no change in ATP status was seen, possibly owing to the compensatory increase in glycolytic ATP production [39]. Given that cell migration is an energy consuming process and the increased wound closure was only observed in HPDE cells overexpressing Snail that exhibit only a small change in maximal respiration (Figure 5F), one could speculate that the large decrease in maximal respiration and ATP production seen in Panc1-Snail cells limited their migratory capacity. This phenomenon was observed in other models of EMT induction where loss of oxidative metabolism impedes cell migration [64]. The dissociation of mitochondrial ETC content and basal/maximal OCR in HPDE cells with or without Snail might be attributed to decreased electron donation to the ETC complexes and increased energy consumption in other EMT-related processes. Alternatively, that fact that Panc1 cells are already on the EMT spectrum prior to Snail induction may limit any additional gain in migratory capacity (Figure 1).

There have been indications in the literature that the process of EMT promotes chemoresistance in various carcinomas including PDAC [65–68]. While the induction of several transcription factors and increased stemness were suggested as possible mechanisms, it is not clear if EMT-associated metabolic programming or increased glycolysis plays a role [66,67]. In the present study, Snail-induced EMT changes did not alter Panc1 sensitivity to gemcitabine or combination gemcitabine paclitaxel treatment despite marked enhancement of glycolysis and slight reduction on proliferation (Table 1 and Figure 1). Further work is required to probe the effects of EMT-related metabolic reprogramming on chemo-sensitivity in additional PDAC cell lines and with other first-line chemotherapeutic agents.

#### **5. Conclusions**

Collectively, Snail overexpression in the PDAC cell line Panc1 and in non-tumorigenic HPDE cells resulted in the induction of EMT and a range of accompanying metabolic changes. In both cell lines Snail overexpression resulted in increased glucose uptake and lactate production, as well as reductions in mitochondrial ETC protein content. Additionally, Snail overexpression caused decreased carbon flux from glucose in the TCA cycle in Panc1 cells only, with no change in Snail-overexpressing HPDE cells. Despite the induction of EMT status and detection of metabolic reprogramming of glucose metabolism in Panc1 cells Snail-overexpression did not result in enhanced resistance to gemcitabine or combination

gemcitabine/paclitaxel therapy. This work highlights the role that Snail plays as an effector of EMT and its role in the induction of metabolic reprogramming in PDAC. Further research to uncover specific changes in metabolic enzymes, pathways and energetic profiles that are essential to EMT in PDAC is required to allow therapeutic interventions from a metabolic angle.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2077-0383/8/6/822/s1, Table S1: 13C metabolic flux analysis results.

**Author Contributions:** Conceptualization, M.L., L.-E.Q., N.T.; methodology, M.L., S.E.H., G.S., B.P.W., E.D., B.O., and L.-E.Q.; writing—original draft preparation, M.L., S.E.H. and L.-E.Q.; writing—review and editing, S.E.H. and N.T.; funding acquisition, L.-E.Q. and N.T.

**Funding:** This work was supported by an Australian Research Council Future Fellowship to N.T. M.L. was supported by an UNSW International Postgraduate Award, and L.-E.Q. by the Judith and Coffey Fund.

**Acknowledgments:** The authors wish to thank Phoebe Phillips for kind donation of the human pancreatic ductal epithelial cells. We also thank the staff at the Biomedical Imaging Facility at UNSW for their technical assistance. We gratefully acknowledge subsidized access to the Bioanalytical Mass Spectrometry Facility at UNSW, supported by the National Collaborative Research Infrastructure Scheme

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Phosphoproteomic Profiling Identifies Aberrant Activation of Integrin Signaling in Aggressive Non-Type Bladder Carcinoma**

**Barnali Deb 1,2, Vinuth N. Puttamallesh 1,3, Kirti Gondkar 1,3, Jean P. Thiery 4,5,6,\*, Harsha Gowda 1,\* and Prashant Kumar 1,2,\***

Received: 27 March 2019; Accepted: 23 April 2019; Published: 17 May 2019

**Abstract:** Bladder carcinoma is highly heterogeneous and its complex molecular landscape; thus, poses a significant challenge for resolving an effective treatment in metastatic tumors. We computed the epithelial-mesenchymal transition (EMT) scores of three bladder carcinoma subtypes—luminal, basal, and non-type. The EMT score of the non-type indicated a "mesenchymal-like" phenotype, which correlates with a relatively more aggressive form of carcinoma, typified by an increased migration and invasion. To identify the altered signaling pathways potentially regulating this EMT phenotype in bladder cancer cell lines, we utilized liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based phosphoproteomic approach. Bioinformatics analyses were carried out to determine the activated pathways, networks, and functions in bladder carcinoma cell lines. A total of 3125 proteins were identified, with 289 signature proteins noted to be differentially phosphorylated (*p* ≤ 0.05) in the non-type cell lines. The integrin pathway was significantly enriched and five major proteins (TLN1, CTTN, CRKL, ZYX and BCAR3) regulating cell motility and invasion were hyperphosphorylated. Our study reveals GSK3A/B and CDK1 as promising druggable targets for the non-type molecular subtype, which could improve the treatment outcomes for aggressive bladder carcinoma.

**Keywords:** urothelial cancer; phosphoproteomics; activated pathways; ingenuity pathway analysis; molecular subtypes

#### **1. Introduction**

Advanced cancer therapeutics demands a detailed understanding of the alteredmechanisms operating in malignant disease. Notably, in bladder carcinoma, the current treatment regimens and interventions are mostly determined by the clinicopathological characteristics of the tumor. Bladder carcinoma is categorized into two clinicopathologically distinct subgroups: non-muscle-invasive bladder cancer (NMIBC) and muscle-invasive bladder cancer (MIBC). At diagnosis, 75% of bladder carcinoma is NMIBC, whereas 20% to 25% is MIBC. About 50% to 70% of NMIBC, including Ta- and T1-stage tumors, frequently recur [1], whereas 10% to 15% progress to MIBC (T2, T3, and T4 stages) [2]. Unlike most other carcinomas, carcinoma in situ (Tis) tumors, a superficial tumor, progress very rapidly to MIBC. Emerging evidence suggests

that tumor heterogeneity alters treatment response and confers resistance, leading to recurrence. Thus, a detailed understanding of the clinicopathological subgroups and their molecular alterations is of paramount interest, especially for targeted therapy. Identification of distinct molecular subtypes of NMIBC and MIBC could highlight clinically relevant signaling pathways.

Clinical staging is still not precisely defined for bladder carcinoma because of its highly heterogeneous nature. The advent of whole-genomic and transcriptomic approaches has transformed our understanding of cancer heterogeneity. Whole-genome sequencing and expression profiling, have uncovered molecular subtypes that relate to distinct clinical entities in various cancers. For instance, several studies have identified intrinsic subtypes of MIBC that resemble the intrinsic subtypes of breast cancer [3–5]. Comprehensive studies have defined the molecular subtypes in bladder carcinoma, mainly by analyzing gene expression data, specific mutations, and copy-number changes with recent analyses suggesting that specific genetic alterations in fibroblast growth factor receptor 3 (*FGFR3*), tumor protein P53 (*TP53*), retinoblastoma 1 (*RB1*), Erb-B2 receptor tyrosine kinase 2 (*ERBB2*), nuclear-factor, erythroid 2 like 2 (*NFE2L2*), and lysine demethylase 6A (*KDM6A*/*UTX*) are enriched in different molecular subtypes and are clinically actionable [6–8].

So far, the intrinsic subtypes of bladder carcinoma have been categorized into two to seven classes by different groups. A broad stratification of high-grade bladder carcinoma was described by the University of North Carolina (UNC), where the subtypes (luminal and basal-like) reflected the hallmarks of breast cancer [4]. This study classified luminal and basal subtypes of bladder carcinoma using a 47-gene signature. Another similar molecular classification (MD Anderson) using gene expression profiling, restricted to MIBCs, found three subgroups: basal, luminal, and p53-like MIBCs [5]; the inclusion of the "p53-like" subtype was based on the mRNA expression of p53, which can predict MIBC chemoresistance. The Cancer Genome Atlas (TCGA) research network also categorized the intrinsic subtypes of bladder carcinoma into four clusters (Clusters I–IV) based on an integrated analysis of the mRNA, miRNA, and proteome expression profiling of 129 tumor samples. However, the most extensive classification by LUND University was conducted on gene expression profiling data for 308 tumor tissue samples, and resulted in five intrinsic types of carcinoma, designated Urobasal A, Genomically unstable, Urobasal B, Squamous cell carcinoma-like, and Infiltrated [9]. The LUND molecular subtypes were further characterized into seven gene signature (subdividing Urobasal A and Genomically instable subtypes into two groups each), based on the differential expression of genes associated with biological processes [10]. Using a 117-gene signature, Hedegaard et al. [11] classified only NMIBCs into three subclasses (Class I–III) with basal- and luminal-like characteristics as well as clinical features. A later meta-analysis on the molecular subtypes (established by UNC, MDA, TCGA and LUND) revealed some consensus patterning among the molecular subtyping strategies [10]. However, there remain inaccuracies in the current staging systems, and the clinicopathological features may also influence treatment selection [12]. Such inadequacies are to be expected due to the low sample size in these studies and a large, collaborative study is needed to identify the molecular signatures before assigning treatment with conventional chemotherapy.

In a study by Warrick et al., bladder carcinoma cell lines were used to study the molecular subtypes in a relevant in vitro model. A detailed analysis was carried out on the publicly available data for the 27 bladder carcinoma cell lines (CCLE database) [13] as per the TCGA study. The panel of bladder carcinoma cell lines was classified into three subtypes based on gene expression clusters: luminal, basal, and "non-type" [14]. The non-type subtype displayed low expression of both luminal and basal markers. We further included these three subtypes for our deeper understanding of the molecular events involved in bladder carcinoma.

We employed liquid chromatography tandem-mass spectrometry (LC-MS/MS)-based approach to identify differentially expressed proteins in all three distinct subtypes. We also sought to evaluate pathways enriched in luminal, basal, and non-type subtypes. We undertook an extensive quantitative phosphoproteomic analysis of six bladder carcinoma cell lines (RT112, SW780, VMCUB-1, T24, J82, and UMUC3) and compared the phosphoproteomic data of bladder carcinoma cell lines with a non-neoplastic bladder cell line (TERT-NHUC) [15]. The six cell lines were previously assigned to luminal (RT112 and SW780), basal (VMCUB-1), and non-type (UMUC3, J82 and T24) molecular subtypes [14]. To our knowledge, this is the first phosphoproteomic profiling of bladder carcinoma and it may provide a global platform to identify the complex kinase-driven signaling events in bladder carcinoma.

#### **2. Experimental Section**

#### *2.1. Cell Culture*

Bladder carcinoma cell lines, SW780, RT112, VMCUB-1, T24, J82, and UMUC3 were cultured in dulbecco's modified eagle's medium (DMEM; with high glucose, HIMEDIA), supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin. Cells were maintained in a humidified incubator at 37 ◦C and 5% CO2. The non-neoplastic bladder cell line (TERT-NHUC) was cultured in KGM-Gold Media (Lonza Group, Allendale, NJ, USA) supplemented with hydrocortisone, transferrin, epinephrine, gentamicin, amphotericin B (GA-1000), bovine pituitary extract (BPE), recombinant human epidermal growth factor (rhEGF), and insulin in a humidified incubator at 37 ◦C and 5% CO2.

TERT-NHUC cells were kindly provided by Prof. M.A. Knowles, University of Leeds, UK [15].

#### *2.2. Cell Lysis, Protein Extraction, and Digestion*

Cell lines were grown to 70% confluence, starved in serum-free medium for 12 h, and then lysed in cell lysis buffer (2% SDS, 5 mM sodium fluoride, 1 mM β-glycerophosphate, 1 mM sodium orthovanadate in 50 mM triethyl ammonium bicarbonate (TEABC)). Protein concentration was estimated using the BCA method (Pierce; Waltham, MA, USA). Equal amount of protein from each cell line was used for protein digestion. Cell lysates were reduced and alkylated using 5 mM dithiothreitol (DTT) and 20 mM iodoacetamide (IAA), respectively. Samples were then digested overnight at 37 ◦C using trypsin (1:20) (Promega; San Luis Obispo, CA, USA).

#### *2.3. TMT Labeling*

Before labeling, peptides were reconstituted in 50 mM TEABC (pH 8.0). Equal amounts of peptides from each condition were used for 10-plex Tandem Mass Tag (TMT) labeling (Thermo Fisher Scientific, Rockford, IL, USA). Labeling was carried out as per the manufacturer's protocol. TERT-NHUC was labeled with the channel 126, and SW780, RT112, VMCUB-1, T24, J82, and UMUC3 were labeled with 130N, 130C, 129C, 131, 128C, and 128N, respectively. The reaction was quenched by the addition of 8 μl of 5% hydroxylamine for 15 min at room temperature, and the samples were dried.

#### *2.4. Basic pH RPLC (bRPLC) and TiO2-Based Phosphopeptide Enrichment*

TMT-labeled peptides were pooled and fractionated using an Agilent 1100 high-pH reverse-phase liquid chromatography (RPLC) system with a flow rate of 1 mL/min, as described earlier [16]. Briefly, labeled peptides, resuspended in 1 mL of bRPLC solvent (10 mM TEABC pH 8.4), were fractionated using high-pH, reverse-phase XBridge C18 columns (5 μm, 250 × 4.6 mm; Waters Corporation, Milford, MA, USA), with an increasing gradient of bRPLC solvent B (10 mM TEABC in 90% ACN, pH 8.4). A total of 96 fractions were collected in a 96-well plate containing 0.1% formic acid. The fractions were then concatenated into 12 pools, vacuum-dried, and subjected to TiO2-based phosphopeptide enrichment [17]. The enriched phosphopeptides were eluted thrice into microfuge tubes with 40 μL of 2% ammonia solution containing 10 μL of 20% trifluoroacetic acid. The peptides were then dried, resuspended in 30 μL of 0.1% TFA, and desalted using C18 Stage Tips (Thermo Fisher Scientific). The eluted peptides were then subjected to LC-MS/MS analysis.

#### *2.5. LC-MS*/*MS Analysis*

Phosphoproteomic analyses were performed on an Orbitrap Fusion Tribrid mass spectrometer (Thermo Fisher Scientific) interfaced with an Easy-nLC II nanoflow liquid chromatography system (Thermo Fisher Scientific), as described earlier [18]. Briefly, each fraction was reconstituted in Solvent A (0.1% Formic acid) and loaded on trap column (75 μm × 2 cm) packed with Magic C18 AQ (Michrom Bioresources, Inc., Auburn, CA, USA). Peptides were resolved on an analytical column (75 μm × 15 cm) at a flow rate of 350 nL/min using a linear gradient from 5% to 60% ACN in a 120 min run. MS data were acquired using scan range of 400–1,600 m/z at mass resolution of 120,000, and MS/MS data were acquired using resolution of 30,000 at m/z of 400. HCD fragmentation for MS/MS analysis was carried out with isolation width of 2 m/z and normalized collision energy of 34%. Data-dependent acquisition was carried out where the most intense precursor ions were detected.

#### *2.6. Data Analysis*

The MS-derived data were searched using Mascot and Sequest HT search engines with Proteome Discoverer 2.0 (Thermo Fisher Scientific). Phosphopeptide-enriched fractions from each replicate were searched against the RefSeq protein database (version 70), National Institutes of Health, USA with carbamidomethylation of cysteine residues as a fixed modification. Oxidation of methionine; the phosphorylations of serine, threonine, and tyrosine; and the deamidation of asparagine and glutamine were selected as dynamic modifications. Trypsin was set as the protease and a maximum of one missed cleavage was allowed. Precursor mass tolerance was set to 20 ppm, and a fragment mass tolerance of 0.05 Da was allowed. All peptide-spectrum matches (PSMs) were identified at a 1% false-discovery rate (FDR). Posterior error probability was calculated for individual PSMs using a percolator, providing statistical confidence for each spectral match. The probability of phosphorylation for each site was calculated by the phosphoRS node in Proteome Discoverer. Only phosphopeptides with >75% site localization probability were considered for further analysis. A 2-fold cut-off was used for dysregulated phosphopeptides compared to TERT-NHUC cell line.

#### *2.7. Scratch Wound Assay*

For each cell line, 3 × 105 cells were seeded into the wells of a 6-well plate. The experiment was conducted in triplicate. When the cells had reached 90% to 95% confluency, a scratch was made with a 20-μl pipette tip. Cells were washed gently with media, incubated with PBS, and then imaged. Cells were maintained in DMEM supplemented with 2% FBS. Cells were observed and imaged at 0 and 24 h. The rate of migration was calculated using ImageJ macros (v1.50i; National Institutes of Health, Maryland, USA). A *t*-test was used to calculate significance and the area covered by the migrating cells was plotted in a bar graph.

#### *2.8. Invasion Assay*

The transwell system was used to measure invasion, as reported earlier [19]. Briefly, 20,000 cells per 500μl of serum free media were seeded onto the Matrigel-coated PET membrane (BD Bio Coat Matrigel Invasion Chamber; BD Biosciences, CA, USA) in the upper compartment, while the lower compartment was filled with complete growth media. Plates were maintained at 37 ◦C for 48 h. After the incubation period, the upper membrane surface was wiped with a cotton-tip applicator to remove non-migratory cells. Cells that migrated to underside of membrane were fixed and stained using 4% methylene blue in 50% methanol. The number of cells that penetrated was counted for randomly selected viewing fields at 10× magnification using CX 41 Olympus Polarizing Trinocular microscope. The counting of cells was done using ImageJ multi-point cell counter (v1.50i; National Institutes of Health, Maryland, USA). The experiment was conducted in replicates.

#### *2.9. Kinome Map*

The kinome map was built using the KinMap tool (http://www.kinhub.org/kinmap/index.html). The list of identified kinases was searched and the relevant kinases are highlighted on the kinome map. Dysregulated kinases (hyper- or hypo-phosphorylated in any cell line) are also depicted.

#### *2.10. Clustering the Molecular Subtypes*

All quantified peptides (across all cell lines and in both replicates) were considered to identify differentially phosphorylated peptides. Fold-change was calculated by taking the intensity ratios with TERT-NHUC as the control. A *t*-test was used to determine the differentially phosphorylated peptides in the non-type subtype as compared with the luminal/basal subtype and the control cells. Differentially phosphorylated peptides (*p* ≤ 0.05) were further considered for supervised clustering using Perseus data analysis tool v<IP_ADDRESS>, Martinsried, Germany (http://www.biochem.mpg.de/5111810/perseuss) [20].

#### *2.11. Ingenuity Pathway Analysis*

Pathway analysis was completed for significantly dysregulated phosphopeptides in the non-type subtype using the Ingenuity Pathway Analysis tool (IPA Build: 460209M; Version: 39480507; IPA, Qiagen, Redwood City, CA, USA). The Core Analysis module in the IPA Ingenuity Knowledge Base reference repository was used to predict relationships of differentially phosphorylated molecules in our dataset. IPA was used to overlay our input dataset with that of the knowledge base. Canonical pathway analyses identified the top canonical pathways significantly enriched in our datasets. The "integrin pathway" was enriched in the total quantified dataset across all tested cell lines and with a confident phosphosite assignment. Significantly enriched functions and networks were also obtained. The IPA analysis results were schematically replicated using Adobe Illustrator (vCS5.1; Adobe Systems, San Jose, CA, USA).

#### *2.12. Kinase-Substrate Enrichment Analysis*

Kinase-substrate enrichment analysis was done using the online KSEA tool (https://casecpb. shinyapps.io/ksea/) [21]. Proteins with differentially phosphorylated sites in the non-type subtype were used for the input file, and analyzed using PhosphoSitePlus and NetworKIN as the background datasets. The *p*-value cutoff (for plot) and number of substrates cutoff were set to 0.05 and 5, respectively.

#### *2.13. Motif Analysis*

Motifs enriched in the non-type subtype were analyzed using the motif-x algorithm v<IP_ADDRESS> (http://motif-x.med.harvard.edu). The "phospho-window" (7 amino acid residues on either side of the phosphorylated residue) was extracted from the RefSeq database. Significance threshold was set to *p* < 0.001. The minimum occurrence of the motif was set to 20 for pSer peptides and 10 for pThr against an IPI (International Proteome Index) human proteome background with the central character as "s" and "t", respectively.

#### *2.14. Immunohistochemistry*

High-grade and low-grade bladder cancer FFPE tissue sections were obtained from Kidwai Institute of Molecular Oncology after informed patient consent. IHC was carried out on both the cases. Briefly, sections were deparaffinized and antigen retrieval was carried by incubating sections in antigen retrieval buffer (0.01 M Trisodium citrate buffer, pH 6) for 20 minutes. Endogenous peroxidases were quenched using (1:1) methanol: chloroform solution followed by washes with PBS plus 0.05% Tween-20. The sections were blocked using 5% goat serum to avoid non-specific binding of primary antibody for 30 minutes. Further sections were incubated with primary anti-Talin1 (S425) (Abcam, Cambridge, United Kingdom) antibody at 1:200 dilutions overnight at 4 ◦C in a humidified chamber. Next day, the sections were washed thrice with wash buffer and incubated with appropriate horseradish peroxidase conjugated rabbit secondary antibody for 30 minutes at room temperature. Excess secondary antibody was removed using wash buffer followed by addition of DAB substrate. The signal was developed using DAB chromogen (DAKO, Glostrup, Denmark) and counterstained by hematoxylin. The immunohistochemical labeling was assessed by an experienced pathologist. Images were taken at 10× on Olympus DP-21 microscope.

#### **3. Results**

#### *3.1. EMT Scores of the Molecular Subtypes of Bladder Carcinoma Cell Lines*

EMT (−1.0 to +1.0) scores were computed to estimate the EMT phenotype of the molecular subtypes of bladder carcinoma cell lines (Supplementary Table S1) [3]. To compute the EMT score in bladder carcinoma cell lines, we adopted a similar approach to that used in ssGSEA [22]. Empirical cumulative distribution function (ECDF) was estimated for Epithelial and Mesenchymal gene sets. The 2KS test was employed to compute the difference between the Mesenchymal ECDF (ECDFMes) and the epithelial ECDF (ECDFEpi). The EMT signature specific to bladder was curated and applied to single sample gene set enrichment analysis (ssGSEA) to provide a gross assessment for the EMT phenotype for each cell line. The BinReg EMT signature was used to predict the EMT phenotype in the cell lines [3,23]. A cell line with a positive EMT score exhibits a more "mesenchymal-like" phenotype, whereas a negative EMT score reflected a more "epithelial-like" phenotype. EMT scores suggested that non-type cell lines were "mesenchymal-like", whereas the luminal and basal cell lines had an "epithelial-like" characteristic (Figure 1a).

#### *3.2. Non-Type Bladder Carcinoma Has Increased Migration and Invasive Ability*

We performed functional assays to determine and compare the tumorigenic properties of the non-type cell lines with the luminal/basal subtype. Using matrigel invasion assay, we estimated that the non-type were significantly more invasive as compared to the luminal (*p* = 6.9 × 10−14) and basal subtypes (*p* = 0.03) (Figure 1b,c). We found increased migration of non-type bladder carcinoma cell lines as compared to luminal (*p* = 0.002) and basal (*p* = 0.007) subtypes (Figure 1d,e).

#### *3.3. Phosphoproteomic Analysis of Bladder Carcinoma Cell Lines*

Phosphoproteomic profiling was conducted for six bladder carcinoma cell lines (RT112, SW780, VMCUB-1, J82, T24, and UMUC3) and the non-neoplastic bladder cell line (TERT-NHUC) using an established workflow for the enrichment of phosphopeptides (Figure 2). To increase the reliability of our phosphoproteomic analyses, we have included technical replicates in our study. The MS data were processed and searched against databases using SEQUEST-HT and MASCOT algorithms using the Proteome Discoverer 2.0 platform. Using an FDR cutoff of 1%, 10,979 phosphopeptides were identified, corresponding to 3125 proteins (Supplementary Table S2). The probability of Ser/Thr/Tyr phosphorylation sites on each peptide was calculated by the phosphoRS algorithm using a cut-off of >75%. We identified 4846 unique phosphopeptides corresponding to 2270 proteins (Supplementary Table S3), with a total of 4878 phosphosites identified: 4271 Ser residues, 574 Thr residues, and 33 Tyr residues (Figure 3a).

**Figure 1.** (**a**) Epithelial-mesenchymal transition score of the molecular subtypes of bladder carcinoma cell lines. (**b**) Invasion assay. Cells which invaded through the Matrigel were stained using methylene blue and imaged at 10× magnification. (**c**) Quantitative analysis of the number of cells which have invaded in each molecular subtype. (**d**) Scratch wound assay. Migration rate of the molecular subtypes of bladder carcinoma cell lines (magnification: 10×). (**e**) Quantitative analysis of scratch wound assay of the molecular subtype of bladder carcinoma.

**Figure 2.** Workflow illustrating the quantitative phosphoproteomics analysis of bladder carcinoma cell lines. For sample processing, proteins were extracted from the bladder carcinoma and non-neoplastic cell lines and digested using trypsin. Each cell line was tagged using the tandem mass tags TMT labeling kit and lyophilized and enriched using the phosphopeptide enrichment protocol (titanium dioxide enrichment). The samples were run on Orbitrap Fusion Tribrid Mass Spectrometer and MS2-based quantitation was achieved. The files were searched against Mascot and Sequest HT search engines. PhosphoRS node was used for the phosphosite assignment. Data was acquired in technical replicates.

**Figure 3.** (**a**) Pie chart depicting the phosphosites identified in the study. Identification of the unique phosphosites (percentage of unique serine, threonine, and tyrosine sites identified in the study is mentioned). (**b**) Kinome map depicting the identified kinases in the dataset. Kinases are highlighted as in the inset legend. The map was built using the KinMap tool.

#### *3.4. Kinases Enriched in the Phosphoproteomics Dataset*

A total of 151 kinases were enriched in the phosphoproteomic analysis of the bladder carcinoma cell lines. Of these, 35 were CMGC (which refers to the CDK, MAPK, GSK3, and CLK set of families); 26 were of the homologs of yeast Sterile 7, Sterile 11, and Sterile 20 kinases (STE); 24 were from the protein kinase A, G, and C families (AGC); 20 were calmodulin/calcium regulated kinases (CAMK); 16 were tyrosine kinases (TK); 12 were tyrosine-like kinases (TLK); seven were casein kinases (CK1); and 11 were atypical kinases. In our dataset, 16 kinases were hyperphosphorylated and 14 were hypophosphorylated (in at least one cell line) (Figure 3b).

#### *3.5. Unique Phosphorylation Signature Identified for the Non-Type Subtype of Bladder Carcinoma*

In a recent detailed analysis, Warrick et al. used agglomerative methods (using the expression of a subtype-specific gene list from the TCGA study) to characterize molecular subtypes of bladder carcinoma cell lines [14]. Based on that report, we selected six bladder carcinoma cell lines to identify differential changes in the signaling pathways of the non-type molecular subtype (Table 1).

**Table 1.** General characterization and molecular subtypes of the bladder carcinoma cell lines.

a UCC, urothelial cell carcinoma; EC, epithelial carcinoma; b Molecular subtypes are from Warrick et al., 2016 [14].

Unsupervised hierarchical clustering was employed to cluster quantified phosphopeptides across normal and malignant bladder cell lines (Supplementary Table S4). The luminal subtype showed distinct clusters (RT112 and SW780), whereas the non-type cell lines (J82, UMUC3 and T24) were clustered along with the basal type (VMCUB-1) (Supplementary Figure S1). Thus, we sought to use the phosphorylation pattern or signature in the non-type subtype to explain what drives the mesenchymal phenotype. 375 differentially phosphorylated peptides corresponding to 289 proteins were identified in the non-type cell lines as compared to the luminal/basal subtype (*p* ≤ 0.05) (Supplementary Table S5). Supervised clustering of the 375 differentially phosphorylated peptides revealed distinct phosphopeptide signature. Two distinct clusters were observed: the non-type cell lines were clustered away from the luminal/basal cells (Figure 4).

**Figure 4.** (**a**) Supervised clustering of the molecular subtypes of bladder carcinoma cell lines. A *t*-test conducted on the quantified data (across all cell lines) indicated that 375 peptides (corresponding to 289 proteins) were differentially phosphorylated in the non-type cell lines (T24, J82, and UMUC3) as compared with the luminal/basal subtype (SW780, RT112, and VMCUB-1) (*p* ≤ 0.05). (**b**) Canonical pathways enriched in the non-type subtype of bladder carcinoma. IPA analysis identified the "integrin pathway" as the most significantly enriched pathway in the aggressive non-type bladder carcinoma subtype.

#### *3.6. Ingenuity Pathway Analysis (IPA) Identifies Aberrant Activation of Pathways in the Molecular Subtypes of Bladder Carcinoma*

IPA of 46 hyperphosphorylated peptides identified significant enrichment of "integrin pathway" components (*p* = 6.29E-04; *z*-score = 2.00) (Figure 4b). Five proteins involved in cell motility had significantly increased phosphorylation in the non-type subtype as compared with the basal and luminal subtypes as well as the TERT-NHUC (control) cells (Supplementary Table S6). TLN1 was 3-fold hyperphosphorylated in the non-type and 2-fold in the basal. We further checked the hyperphosphorylation of TLN1 (S425) in bladder tumor sections. Immunohistochemical staining showed hyperphosphorylation of TLN1 (S425) in high-grade tumors as compared to low-grade tumors (Supplementary Figure S2). TLN1 is a major mediator of the integrin activation and crosstalk. It also stabilizes extracellular matrix-actin linkage [24]. CTTN and CRKL were 1.5-fold hyperphosphorylated in the non-type subtype. CTTN phosphorylation is reported to mediate and accompany integrin mediated cell adhesion to the extracellular matrix [25]. The translocation of CRKL to the focal adhesion activates integrin induced downstream signaling through Src family of kinases and further mediates cell migration [26]. ZYX and BCAR3 were 1.5-fold hyperphosphorylated in the non-type subtype whereas 0.5-fold hypophosphorylated in the luminal subtype. ZYX regulates cellular movement and binds strongly to the focal adhesions. Integrins bind its cytoplasmic tail to ZYX through alpha-actinin [27]. BCAR3 along with its binding partner BCAR1 leads to changes in cellular morphology, motility, and adhesion by the transduction of integrin signaling [28]. To check the global phosphorylation status in aggressive non-type we carried out a pathway analysis with the set of all the phosphopeptides quantified across all cell lines (2806 phosphopeptides). The analysis displayed integrin pathway enrichment which comprises global proposition of the other proteins (the other proteins which are identified/dysregulated in our study) (*p* = 5.06E-04; *z*-score = 5.303) with the total inclusion of 49 proteins in the signaling network (Figure 5).

**Figure 5.** Schematic diagram of enriched integrin signaling pathway in non-type subtype of bladder carcinoma cell lines. Ingenuity pathway analysis lead to the identification of integrin signaling to be most enriched in the aggressive molecular subtype. The dysregulated proteins, the number of phosphopeptides, and the phosphosites identified are highlighted in the pathway.

For the luminal subtype, there was an enrichment of "phospholipase C signaling" components (*p* = 6.29E-03; *z*-score = 2.00), including four differentially phosphorylated proteins—Rho guanine nucleotide exchange factor 16 (AREGEF16), actin-related protein 2/3 complex subunit 1B (ARPC1B), cdc42 effector protein 1 (CDC42EP1), and stathmin1 (STMN1). The basal subtype showed enrichment of "signaling by Rho family GTPase" components (*p* = 3.68E-03; *z*-score = 2.00), and included the differential phosphorylation of neuroblast differentiation-associated protein (AHNAK), AREGEF16, myristoylated alanine-rich C-kinase substrate (MARCKS), and proto-oncogene tyrosine-protein kinase Src (SRC) (Supplementary Figure S3).

#### *3.7. Regulatory Interaction Network Enriched in the Non-Type Subtype*

Network analysis appeared to link key molecules not identified in our study which may be important regulators of the network along with the non-type subtype specific phosphoproteins. The most enriched network in the non-type subtype contained 25 proteins from our dataset (score = 46) (Supplementary Figure S4a). The primary hub proteins included protein kinase A catalytic subunit (PRKACA), Rac-alpha serine/threonine protein kinase (AKT) and nuclear factor NF-kappa-B (NFκB). PRKACA interacts with protein kinase A regulatory subunit 1A (PRKAR1A), ADP ribosylation factor guanine nucleotide exchange factor 2 (ARFGEF2) and A-kinase anchor protein (AKAP11). AKT interacts which TBC2 domain family 4 (TBC2D4) and interferon-induced, double-stranded RNA activated protein (EIF2AK2), desmoplakin (DSP), SRC, and integrin beta-4 (ITGB4). NFκB interacts with TGFβ-activated kinase 1 (TGFB1), MAP3K7-binding protein 2 (TAB2), TBC2D4, casein kinase II subunit beta (CSNK2B), EIF2AK2, CDC42, CDC42EP1, DNA polymerase alpha subunit B (PLOA2), sentrin-specific protease 1 (SENP1), PKA, and immunoglobulin G (IgG). Cell motility was the most enriched function determined by the IPA analysis. 21 proteins related to cell motility were enriched in the non-type subtype of bladder carcinoma (Supplementary Figure S4b).

#### *3.8. CDK1 (Cyclin-Dependent Kinase 1) and GSK3A*/*GSK3B (Glycogen Synthase Kinase 3A and 3B) Are the Predicted Activated Kinases in Non-Type Subtype of Bladder Carcinoma Cell Lines*

Eight kinases—GSK3A, GSK3B, CDK1, ribosomal protein S6 kinase alpha-1 (RPS6KA1), ribosomal protein S6 kinase alpha-1 (RPS6KB1), serine/threonine protein kinase Sgk 1 (SGK1), TGFβ receptor type-2 (TGFBR2), and 5 -AMP-activated protein kinase catalytic subunit alpha-2 (PRKAA2)—were predicted to be significantly enriched (Supplementary Table S7). GSK3A (*p* = 0.003; *z*-score = 2.7) and GSK3B (*p* = 0.01; *z*-score = 2.3) were predicted to be activated and responsible for the phosphorylation of four downstream proteins (tuberin (TSC2), AKAP11, MAP1B, and TBC1 domain family 4 (TBC1D4)). GSK3A was predicted to independently phosphorylate translocon-associated protein subunit alpha (SSR1), whereas GSK3B was predicted to be upstream of neurogenic locus notch homolog protein 2 (NOTCH2) and transcription factor RelB (RELB). CDK1 was also predicted as the upstream activated kinase (*p* = 0.01; *z*-score = 2.06) for 23 proteins at single, double, or triple phosphorylation sites (Supplementary Table S8). RPS6KA1, RPS6KB1, and SGK1 were predicted as the negatively regulated kinases (*p* = 0.01; *z*-score = −2.10) (Figure 6a–c).

**Figure 6.** (**a**) Predicted upstream kinases enriched in the non-type subtype of bladder carcinoma cell lines. Graph showing the positively regulated upstream kinases (red bars) predicted to be activated in the non-type subtype. (**b**) Substrates of GSK3A/3B and (**c**) CDK1 enriched in non-type molecular subtype depicted by a schematic diagram. The respective phosphosites of the substrates identified are also highlighted. (**d**) Motif analysis of differentially phosphorylated peptides of non-type subtype bladder carcinoma. Serine and threonine motifs identified in the non-type subtype of bladder carcinoma.

#### *3.9. Active Proline-Directed Motifs Identified in the Non-Type Molecular Subtype*

In the non-type subtype, five serine and one threonine phosphorylated motifs were significantly differentially phosphorylated (among the 289 signature proteins). A maximum fold change of 36 was calculated for the pSDxE serine phosphorylated motif. The consensus motifs PxpSP and pSP motifs also showed a high fold change in phosphorylation of 7.5-fold. Supplementary Table S9 enlists the differentially phosphorylated protein sites matching to the consensus motifs. Only one motif for the threonine sites (pTP) was identified in the dataset (10-fold enriched among 24 phosphopeptides) (Figure 6d).

#### **4. Discussion**

A molecular subtype is an "intrinsic" feature that defines both clinical and biological stratification of a given tumor. The molecular subtypes of bladder carcinoma have not been extensively studied. In addition, switching between molecular subtypes and the multifocal characteristics of these tumors are two phenomena that preclude our understanding of the biological and clinical properties of these tumors. Most of the studies have defined the molecular subtypes by analyzing gene expression data based on the enrichment of specific genetic alterations and gene expression profiles. However, a comprehensive and quantitative phosphoproteomics study on bladder carcinoma has still not been reported. Here we provide compelling data that suggests the cellular signaling pathways that contribute to the mesenchymal phenotype in the non-type molecular subtype of bladder carcinoma. Our study identified 10,979 phosphopeptides corresponding to 3125 proteins from a panel of bladder carcinoma cell lines. This study offers new biological insight and identifies potential kinases that could be targeted in bladder carcinoma. We combined bioinformatics approaches with phosphoproteomic analyses of bladder carcinoma cell lines to generate highly enriched cellular signaling pathways. This provided the groundwork for defining the molecular subtypes in bladder carcinoma cell lines, which may help to identify and develop new therapeutic approaches within these model systems. Importantly, the study needs to be extended on large cohort of bladder tumor sections to validate a prognostic and predictive role of the subtype specific signature. This would aid in the clinical assessment of bladder carcinoma patients for appropriate treatment regimen based on the molecular subtypes. The computed EMT score suggests that the non-type molecular subtype is more "mesenchymal-like," whereas the luminal/basal subtypes are "epithelial-like." The non-type subtype cell lines show an increased migratory and invasive phenotype, reflecting typical characteristics of a mesenchymal-like phenotype. Non-type carcinomas may; thus, more readily invade the bladder wall, and this knowledge may provide the necessary evidence to spur a change in the field of muscle invasive bladder carcinoma, particularly in terms of their recurrent nature. Recently, existence of hybrid EMT state has also been reported by various groups [29,30]. A study by Jolly et al. showed that hybrid epithelial/mesenchymal cells adhere higher tumor-initiation capabilities and metastatic potential [31]. Another study by Yadavalli et al. computed the EMT score using transcriptomic datasets. The EMT scores distributed between −0.3 to +0.5, for cancer CTCs suggestive of intermediate phenotypes [32]. In our study, we observed the distribution of EMT scores in the aggressive non-type subtype were distributed between +0.35 to +0.65.

We identified the enrichment of the integrin signaling pathway in the non-type molecular subtype, with differential phosphorylation in five major proteins: TLN1, CTTN, BCAR3, CRKL, and ZYX. Integrins are a family of cell adhesion receptors that mediate cell–matrix and cell–cell interactions [33–35]. A major function of the integrins is to regulate the intracellular signaling cascades that lead to cell proliferation, survival, motility, and migration [36–38]. TLN1 is an adaptor protein critical to integrin signaling. In the non-type molecular subtype, TLN1 is significantly hyperphosphorylated at S425 in the head domain (band 4.1, ezrin, radixin, moesin homology domain) (Supplementary Figure S5), a site also known to be hyperphosphorylated in prostate cancer [39]. TLN1 S425 is phosphorylated by cyclin-dependent kinase 5 (CDK5), which controls the metastatic potential in prostate cancer [40]. S425 acts on integrin β with an "inside-out" mechanism; however, we did not identify a change in integrin β1. Cell migration is stimulated by TLN1 coupled to the integrin cytoplasmic domain: TLN1 connects ligand-bound integrins with the actin cytoskeleton, and this is needed for the catalysis of the focal adhesion-dependent pathways [41] to initiate cell movement. The phosphorylation of TLN1 at S425 limits focal-adhesion turnover, which, in turn, stabilizes the lamellipodia contributing to the sustained cell migration [42]. TLN1 overexpression is also correlated with advanced and aggressive oral squamous cell carcinoma [43]. In concordance with the reported studies, we also observed a hyperphosphorylation of TLN1 (S425) in high grade tumor sections. In our unpublished global proteomics data of bladder carcinoma cell lines, TLN1 is >1.5-fold overexpressed in the non-type cell lines; however, its role in bladder carcinoma and whether it contributes to the "mesenchymal-like" phenotype is unclear. CTTN is an actin-binding protein significantly hyperphosphorylated at S418, S405, and T401 in the non-type molecular subtype. CTTN promotes cell migration by lengthening lamellipodia, increasing the number of filopodia, and increasing protrusion time [44]. CTTN is phosphorylated at S418 and S405 by the serine/threonine extracellular signal-regulated kinases (ERK1/2) [45], which, in turn, stimulates binding of the Arp2/3 complex to N-WASP via its SH3 domain [46]. CTTN phosphorylation is required for actin regulation and is responsible for lamellipodia formation and associated with enhanced cancer cell migration [47]. In bladder carcinoma, CTTN may be involved in the regulation of actin-based extravasation through invadopodia formation [48]. CRKL is another adaptor protein that is overexpressed in many cancers, including bladder carcinoma [49]. It regulates cyclin D1 and modulates ERK1/2 phosphorylation. In bladder carcinoma cell lines, CRKL is hyperphosphorylated at S107; however, the function of this phosphorylation is unclear. The two other differentially phosphorylated proteins identified in our study are ZYX and BCAR3. ZYX is a LIM domain protein with distinct actin polymerization activity. It also modulates cell adhesion and the expression of integrins, controlling cell motility in lung carcinoma. ZYX may also regulate EMT during lung cancer development and may regulate cell–cell adhesion, integrin α5β1 expression, and ECM adhesion [50]. ZYX is overexpressed in breast cancer and positively correlated with breast carcinoma metastasis [51]. Its depletion reduces cell proliferation, motility, and in vivo tumorigenic activity. BCAR3, on the other hand, binds to the C-terminal of its partner p130Cas (BCAR1) and co-ordinates the activation of SRC, a non-receptor tyrosine kinase [52]. The BCAR3–SRC activation axis is governed by PTPα, binding to its SH2 domain to recruit the BCAR3/Cas complex. T789 phosphorylation of PTPα is required for integrin induced-adhesion signaling [53]. The functional roles of ZYX (S259) and BCAR3 (T50) phosphorylation have not been studied in bladder carcinoma, but these proteins may contribute to invasion, as both proteins are involved in cell motility. However, EMT is regulated at several levels through transcriptional controls, epigenetic modifications, translational modifications, and alternative splicing [54–58].

Phospholipase C has a decisive role in skin cancer development [59], and its overexpression in gastric-mucosa cells perhaps offers an opportunity to differentiate gastric cancer and inflammatory lesions [60]. Thus, its overexpression in luminal bladder carcinoma is not unexpected, and may be the result of early pro-inflammatory signaling. Our pathway analysis also identified an enrichment of Rho signaling in the basal cell line. Activated Rho proteins can alter cell behavior as well as cell morphology. In breast cancer, Rho family GTPases regulates cell motility by cytoskeleton remodeling and altering focal adhesions [61]. High RhoA and RhoC expression is correlated with poor tumor differentiation, muscle invasion, and lymph node metastasis in bladder carcinoma [62]. The basal subtype of bladder carcinoma is relatively more aggressive than the luminal subtype, and activated Rho protein signaling may contribute to this aggressive clinical outcome [63].

In our IPA analysis, NFκB was highlighted as the regulatory hub molecule of the most enriched network in the non-type subtype. This transcriptional regulator largely affects the expression of cytokines, chemokines, adhesion molecules, and controls various mechanisms, including inflammation, proliferation, transformation, angiogenesis, invasion, and chemoresistance [64]. Interestingly, β1-integrin contains a unique NFκB-binding site in its promoter region, suggesting a potential regulatory role of NFκB in integrin signaling [65]. Thus, NFκB may regulate the expression of integrins and enhance signaling downstream, affecting cell proliferation and motility.

The non-type significantly enriched in golgi vesicle (0.09%), cytoskeleton (5.1%) and in nuclear compartment (75.2%). In addition, ubiquitin-specific proteases (4.9%), cytoskeleton protein binding (4.5%), and signal transducers (1.1%) were significantly active in non-type bladder carcinoma cell line (Supplementary Figure S6). Through our kinase-substrate enrichment analysis, we predicted that the two most activated kinases were GSK3A and GSK3B, two highly conserved serine/threonine protein kinases. The roles of these two kinases in cancer remain controversial, with little evidence in the literature to support their role in the disease. The Y216-phosphorylated active form of GSK3B is significantly decreased in squamous cell carcinoma [66] and in lung cancer, inhibiting GSK3B increases the expression of the transcriptional repressor Snail, which suppresses E-cadherin and promotes EMT [67]. A potent GSK3 inhibitor LY2090314 has been tested in clinical trials against metastatic pancreatic cancer (NCT0163230) and acute leukemia (NCT01214603) [68]. Others have reported that pharmacological inhibition of GSK3B can induce apoptosis by sorafenib in melanoma cell lines and leukemic cells [69,70]. Yet, in ovarian cancer, GSK3B activity was linked with cell proliferation, and overexpression of its active form can induce CDK1 expression and facilitate cell proliferation [71]. CDK1 is a major mitotic regulating kinase involved in the G2/M phase of the cell cycle, and is upregulated in several cancers and also associated with clinicopathological factors [72–74]. Various CDK inhibitors have been trialed for the treatment of breast cancer and chronic lymphocytic leukemia [75–77]. Nevertheless, the role of CDK1 in bladder carcinoma has yet to be elucidated. GSK3A/B and CDK1 could be potential druggable targets for the aggressive non-type bladder carcinoma. Intriguingly, we identified enrichment of proline-directed motifs (pSP/ pTP) in our dataset, which are known to be active mostly in proliferating cells and phosphorylated by CDK1 during mitosis [78]. We also identified the pSer and pThr motifs that recognize acidic casein kinase II in our dataset [79], yet the importance of these motifs in bladder carcinoma remains unclear.

#### **5. Conclusions**

We present the first comprehensive phosphoproteome of bladder carcinoma cells and offer a potential starting point for further mechanistic studies. Our phosphoproteomic data and bioinformatic analyses suggest two key activated kinases (GSK3A/B and CDK1) as potential druggable targets for the aggressive subtype of bladder carcinoma. We provide a glimpse of the dynamic process of phosphorylation in bladder carcinoma and also the differentially regulated phosphoproteins in this aggressive molecular non-type subtype.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2077-0383/8/5/703/s1, Figure S1: Unsupervised clustering of the bladder carcinoma cell lines of the quantified data across all the cell lines; Figure S2: Immunohistochemistry showing TLN1 (S425) hyperphosphorylation in high-grade tumors as compared to low-grade tumors; Figure S3: Top five canonical pathways enriched by IPA in (a) luminal and (b) basal subtypes of bladder carcinoma; Figure S4: (a) A global network interaction map of the differentially phosphorylated proteins in the non-type subtype. Most enriched network in the aggressive subtype with an enrichment score of 46 (IPA). (b) Cell motility is an intrinsic characteristic of the non-type subtype of bladder carcinoma. Phosphoproteins involved in cell motility, as depicted from the IPA analysis; Figure S5: Protein domain structure for TLN1. The domains of TLN1 and the significantly differentially phosphorylated sites (S425; head domain) in the non-type cell lines; Figure S6: Gene ontology enrichment of the differentially phosphorylated proteins in the non-type molecular subtype. Most of the proteins are nuclear proteins (75.2%). Proteins with kinase activity and signal transducer activity are also significantly enriched (*p* ≤ 0.05). Data were analyzed using FunRich; Table S1: EMT score of bladder carcinoma cell lines molecular subtypes; Table S2: List of phosphopeptides identified in bladder carcinoma cell lines using SequestHT and Mascot search algorithms; Table S3: List of phosphopeptides identified in bladder carcinoma cell lines using SequestHT and Mascot search algorithms with a PhosphoRS probablity of ≥75%; Table S4: List of phosphopeptides identified in bladder carcinoma cell lines using SequestHT and Mascot search algorithms with a PhosphoRS probability of ≥75% and quantified across all the cell lines; Table S5: A total of 375 differentially phosphorylated peptides in the non-type bladder carcinoma molecular subtype; Table S6: Proteins differentially phosphorylated in non-type subtype in integrin signaling (*p*-value ≤ 0.05); Table S7: Kinase scores for the kinase-substrate enrichment analysis; Table S8: Kinase-substrate enrichment scores; and Table S9: Motifs identified by motif-x in the differentially phosphorylated peptides in non-type subtype.

**Author Contributions:** Conceptualization, P.K. and H.G.; methodology, B.D., V.N.P., and K.G.; software, B.D. and V.N.P.; validation, B.D. and K.G.; formal analysis, B.D., H.G., and P.K.; investigation, B.D., H.G., and P.K.; writing—original draft preparation, B.D. and P.K.; writing—review and editing, J.P.T., H.G., and P.K.; visualization, B.D., J.P.T., H.G., and P.K.; supervision, J.P.T., H.G., and P.K.; project administration, P.K.; funding acquisition, P.K. All authors have read and approved the final version of the manuscript.

**Funding:** This research was funded by Department of Science and Technology (DST), Ramanujan Fellowship, Government of India, grant number SB/S2/RJN-077/2015 and the APC was funded by SB/S2/RJN-077/2015, Department of Science and Technology (DST), Ramanujan Fellowship, Government of India.

**Acknowledgments:** We thank the Department of Biotechnology (DBT), Government of India for research support to the Institute of Bioinformatics, Bangalore. PK is a recipient of the Ramanujan Fellowship awarded by Department of Science and Technology (DST), Government of India. BD is a recipient of INSPIRE Fellowship from Department of Science and Technology (DST), Government of India. KG is a recipient of Senior Research Fellowship from University Grants Commission (UGC), Government of India. We also thank Rebecca A. Jackson, Scientific Editor, for English editing.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Clinical Stratification of High-Grade Ovarian Serous Carcinoma Using a Panel of Six Biomarkers**

**Swapnil C. Kamble 1,2, Arijit Sen 3, Rahul D. Dhake 4, Aparna N. Joshi 5, Divya Midha 6 and Sharmila A. Bapat 1,\***

Received: 31 January 2019; Accepted: 21 February 2019; Published: 8 March 2019

**Abstract:** Molecular stratification of high-grade serous ovarian carcinoma (HGSC) for targeted therapy is a pertinent approach in improving prognosis of this highly heterogeneous disease. Enabling the same necessitates identification of class-specific biomarkers and their robust detection in the clinic. We have earlier resolved three discrete molecular HGSC classes associated with distinct functional behavior based on their gene expression patterns, biological networks, and pathways. An important difference revealed was that Class 1 is likely to exhibit cooperative cell migration (CCM), Class 2 undergoes epithelial to mesenchymal transition (EMT), while Class 3 is possibly capable of both modes of migration. In the present study, we define clinical stratification of HGSC tumors through the establishment of standard operating procedures for immunohistochemistry and histochemistry based detection of a panel of biomarkers including TCF21, E-cadherin, PARP1, Slug, AnnexinA2, and hyaluronan. Further development and application of scoring guidelines based on expression of this panel in cell line-derived xenografts, commercial tissue microarrays, and patient tumors led to definitive stratification of samples. Biomarker expression was observed to vary significantly between primary and metastatic tumors suggesting class switching during disease progression. Another interesting feature in the study was of enhanced CCM-marker expression in tumors following disease progression and chemotherapy. These stratification principles and the new information thus generated is the first step towards class-specific personalized therapies in the disease.

**Keywords:** high-grade serous ovarian carcinoma; epithelial-to-mesenchymal transition; molecular stratification; biomarkers; scoring system; immunohistochemistry

#### **1. Introduction**

Personalized therapeutic decisions in cancer necessitate the development of accurate stratification schemes based on mutations and/or association of tumor sub-groups with specific biomarkers and biological functions, besides well-elucidated principles for their detection [1]. Recent resolution of four gastric cancer molecular groups identified predictive amplifications for subtype-specific treatment [2], including *PDCD1LG2* locus for use of nivolumab and pembrolizumab in the Epstein-Barr virus-associated group, EGFR for cetuximab, panitumumab, nimotuzumab, or matuzumab treatment in the chromosome instability group and Aurora kinase A/B inhibitors for treatment of the genomically stable (GS) subgroup [3–6]. Immunohistochemical (IHC) has become a significant tool in clinical diagnostics and is frequently utilized to classify malignant cells [7]. In gastric cancer, a panel of six biomarkers was used in tumor stratification [8,9]. In a similar approach, cancers of the endometrium [10,11], lung [12], triple-negative breast [13], esophagogastric junction carcinomas [14] were stratified into discrete molecular classes using tumor-specific IHC-based biomarkers. Multiplexed IHC for the concurrent detection of a number of biomarkers in lung cancer is increasingly becoming point-of care in treatment [15]. Such translation of molecular information implies the feasibility of similar applications in other tumor types. High-grade serous ovarian cancer (HGSC) represents aggressive tumors characterized by swift metastatic progression and poor patient prognosis [16]. Despite radical surgery and initial response to platinum and taxane based chemotherapy, most patients relapse following median progression-free survival of ~18 months [17,18]. Clinical outcomes vary considerably emphasizing an imminent need to improve therapeutic options. Large-scale molecular analyses have recently identified diverse molecular pathways, mutations, gene expression, morphologies, cell(s) of origin, etc. leading to a systematic understanding of HGSC despite its heterogeneity [19–23]. Our earlier analyses of gene expression datasets also resolved three classes in HGSC that were associated with discrete mechanisms of metastases [24]. Development of targeted therapies now necessitates the establishment of a robust diagnostic pipeline for HGSC stratification. As a first step towards this aim, the present study evaluates the application of six markers using immunohistochemistry (IHC) and histochemistry (HC), the establishment of standard operating procedures (SOPs) and development of a reference human tissue library for these markers along with scoring guidelines for interpretation of marker expression. Further evaluation and application were performed in xenografts and commercial tissue microarrays (TMAs), along with the determination of thresholds for clinical classification in resected primary tumors and secondary metastases and/or cell blocks prepared from ascitic fluid of chemo-naïve and chemo-treated patients were also achieved (Supplementary Figure S1). These efforts define the establishment of diagnostic principles for application in clinical practice.

#### **2. Materials and Methods**

#### *2.1. Sample Collection and Preparation*

Formalin-fixed and paraffin-embedded (FFPE) tissue collection and processing using routine methods following surgery, after obtaining informed consent, were approved by the respective Institutional Review Board of NCCS with project identification code IEC/22/12/2014. All subjects gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the National Centre for Cell Science IEC/22/12/2014. In all, retrospective 96 primary high-grade serous ovarian adenocarcinoma patient cases with information of name, age, grade, stage, and treatment status were selected, who had undergone surgery at the Armed Forces Medical College (Pune, India; 2008–2015), Tata Medical Centre (Kolkata, India; 2013–2014), Jehangir Hospital (Pune, India; 2003–2005), Command Hospital (Pune, India; 2010–2011) and Inlaks & Budhrani Hospital (Morbai Naraindas Budrani Cancer Institute, Pune, India; 2013–2015).

#### *2.2. Animal Studies*

Animal experimentation was in accordance with the rules and regulations of the National Centre for Cell Science (NCCS) Institutional Animal Ethics Committee. The study was approved with project number IAEC/2011/B-163. Xenografts were raised as described earlier [24]. In brief, 2.5 × 106 cells of cell lines OVCAR3, OV90, OVMZ6, A4, CP70, PEO14, and CAOV3 were injected subcutaneously in non-obese diabetic/severe combined immunodeficient (NOD/SCID) mice. Animals were maintained under pathogen-free conditions and assessed every 2 days until the tumor diameter was ~1 cm, whereupon animals were sacrificed and tumors harvested.

#### *2.3. Immunohistochemical (IHC) and Histochemical Staining (HC)*

IHC and HC were performed in 5 μm sections of FFPE blocks fixed by drying at 60 ◦C for at least 1 h in an oven using standard protocol, deparaffinized in xylene and hydrated in ethanol-distilled water gradient. Heat-induced epitope retrieval (HIER) was carried out for 30 min at pH = 9/pH = 6 (Himedia, India). For peroxidase inactivation, sections were incubated in 3% H2O2 for 30 min (Qualigens, MA, USA), followed by 1× Blocking Solution for 10 min (Biogenex, CA, USA) and overnight incubation in primary antibody (Abcam, MA, USA; E-cadherin ready-to-use, Biogenex, CA, USA; Santa Cruz Biotechnology Inc.Texas, USA; Abcam, MA, USA). Sections were washed and incubated with anti-rabbit HRP-conjugate (Jackson ImmunoResearch Laboratories, Inc., PA, USA) or anti-mouse HRP-conjugate (Jackson ImmunoResearch Laboratories, Inc., PA, USA) for 1 h, and color developed with DAB (Thermo Pierce, MA, USA); hematoxylin used as a counterstain. These sections were dehydrated and mounted in DPX (Qualigens, MA, USA). Negative controls were prepared in the absence of primary antibody. IHC methods were standardized for each marker as Standard Operating Procedures (SOPs; Supplementary Dataset 1). SOPs were developed considering the positive and negative expression tissue controls, and secondary antibody control for each batch of IHC run. For HC-based HA detection, test sections were exposed to freshly prepared hyaluronidase (1 mg/mL; Sigma-Aldrich, MA, USA); control slides were incubated in phosphate buffer for 1 h at 37 ◦C. Sections were washed in running water for 10 min and stained with Alcian blue for 30 min (pH = 2.5, Fluka, MA, USA), counterstained with Nuclear Fast Red Solution for 2 min (Sigma-Aldrich, MA, United States), and dehydrated and mounted in DPX (Qualigens, MA, USA). Positive experimental controls included testis (TCF21, PARP1), liver (E-cadherin), lymphocytes (Slug), gall bladder (ANXA2), and small intestine (hyaluronan); negative controls included heart (TCF21, E-cadherin, Slug, ANXA2) and mucosa of the small intestine (PARP1). Slides for human tissues, xenografts, and TMA were reviewed independently; a consensus was reached to establish tissues for reference score.

#### *2.4. Statistical Analysis*

Each observer scored the biomarkers for frequency, intensity, and localization. Computation of these scores led to the derivation of biomarker and class indices, which compared between groups by Pearson correlation using SPSS (version 20, SPSS Inc., Chicago, IL, USA) for Windows to delineate classes. Student's t-test and ANOVA were determined in Microsoft Excel 2016; *p* < 0.05 was considered significant.

#### **3. Results**

#### *3.1. Selection of Class-Specific Biomarkers, Development of SOPs for Detection, a Reference Human Tissue Library and Guidelines for Scoring*

A strong correlation of the transcription factors TCF21 and Slug with Class 1 (Cooperative Cell Migration/CCM-class) and Class 2 (Epithelial-to-Mesenchymal/EMT-class) tumors respectively [15] lent consideration to their inclusion in this study. E-cadherin was selected as a feature of cell-cell adhesion to substantiate CCM class-specific purported biological functions and PARP1 for defects in homologous recombination. Known associations of EMT with AnnexinA2 (ANXA2) and its interactions with Slug, and extra-cellular matrix components including hyaluronan (HA) and its synthesizing genes (*HAS1* and *HAS2*) suggested HA as a candidate marker (Supplementary Table S1, Supplementary Figure S2). Class 3 tumors lacked any unique biological features, hence no specific markers were assigned for their identification. An inability to correlate variations in Vimentin protein levels with either of the transcription factor present in these groups refrained its inclusion. The final

screening biomarker panel thus comprised of TCF21, E-cadherin, PARP1, Slug, ANXA2, and HA (Figure 1A).

**Figure 1.** Class associations and marker scoring guidelines; (**A**) rationale for class-specific biological function based putative marker selection, MET: Mesenchymal-to-epithelial transition, HR: Homologous Recombination mediated DNA Damage Repair, EMT: Epithelial-to-mesenchymal transition; (**B**) schematic of scoring guidelines for IHC based staining of nuclear markers (TCF21, PARP1, Slug), A: Absent, W: Weak, I: Intermediate, S: Strong, Mis: Mislocalised, N: normal localization. A similar approach was used for scoring of membrane markers (E-cadherin, ANXA2) except that sub-cellular location was scored either 1 (cytoplasm) or 2 (cell membrane), while extracellular expression of hyaluronan fibers (evaluated as blue color developed by Alcian blue staining that is lost on hyaluronidase) was scored 1 in distant tumor stroma, and 2 in tumor epithelial cell nests. Scoring and analyses of marker expression in xenografts and TMAs; (**C**) Tissues and markers for Scoring of Frequency–0: A4 (TCF21), 1: OV90 (Slug), 2: OVCAR3 (PARP1), 3: A4 (Slug); scoring of intensity and localization-CAOV3 (TCF21), different regions representing scores of 0–3 (intensity) and 0–2 (localization), Scale bar is 50 μm; (**D**) Representative micrographs of HGSC xenografts for: Row 1-H&E (hematoxylin and eosin) stained section while Rows 2,3,4,7,8 represent IHC-based identification of TCF21, Slug, E-cadherin, PARP1, and ANXA2; Rows 5 and 6 represent HC-based identification of HA fibers in untreated and hyaluronidase digested xenograft sections respectively, scale bar is 50 μm; (**E**) class Indices of xenografts; (**F**,**G**) scatter plots of CICCM vs. CIEMT in xenografts and high-grade serous ovarian carcinoma (HGSC) cases in TMA respectively. DP-double positive; CCM-cooperative cell migration.

Evaluation of any novel marker in tumor stratification necessitates the establishment of standard operating protocols (SOPs) to address pre-analytic (slide coating, tissue selection, fixation, processing), analytic (clone and antibody selection, buffers and instruments for antigen retrieval, antibody/enzyme concentration, duration of incubations at each step, etc.), and post-analytic parameters (interpretation, analyses and reporting of expression in the reference and control tissues). These were established for our panel (Supplementary Dataset 1), along with the development of a reference library based on the Human Protein Atlas (HPA) [25] using appropriate normal human tissues. Three specific metrics associated with IHC/HC detection viz. frequency, intensity, and localization were applied in developing universal guidelines for marker scoring (Figure 1B depicts a schematic for transcription factor marker scoring, while a reference score sheet is provided in Supplementary Figure S3). The subjectivity of analyses due to inter-personal observation variation was minimized by collecting independent scores from five observers followed by a comprehensive pathology review to arrive at a consensus in case of difference in opinions. Specific scoring guidelines for each marker that were thereby agreed on and corresponding healthy tissue included the following:

	- a. TCF21: cardiac myocytes, ovarian stromal cells, and germinal cells of testis represented SFreq 0, 1, and 3 respectively; SFreq = 2 could not be identified in healthy tissues.
	- b. E-cadherin: cardiac myocytes, liver hepatocytes, and prostate epithelial cells represented SFreq 0, 2, and 3 respectively; healthy tissues representing SFreq = 1 could not be identified.
	- c. PARP1: mucosa of the small intestine, cardiac myocytes, germinal basal cells of testis represented SFreq as 0, 1, and 3 respectively; healthy tissues representing SFreq = 2 could not be identified.
	- d. Slug: cardiac myocytes, smooth muscles of the appendix, lymphocytes of the small intestine represented SFreq 0, 1, and 2 respectively; healthy tissues representing SFreq = 3 could not be identified.
	- e. HA: cartilage and sub-mucosa of the small intestine represented SFreq as 2 and 3 respectively; healthy tissues representing SFreq = 0 or 1 could not be identified.
	- f. ANXA2: cardiac myocytes, the somatic muscle of the small intestine, epithelial cells of the gall bladder represented SFreq 0, 1, and 3 respectively; healthy tissues representing SFreq = 2 could not be identified.
	- a. TCF21: cardiac myocytes, ovarian stromal cells, germinal basal cells of testis represented SInt 0, 1, and 2 respectively; SInt = 3 could not be identified in healthy tissues.
	- b. E-cadherin: cardiac myocytes, epithelial cells of the small intestine, epithelial cells of prostate represented SInt 0, 2, and 3 respectively; healthy tissues representing SInt = 1 could not be identified.
	- c. PARP1: mucosa of the small intestine, cardiac myocytes, and germinal basal cells of testis represented SInt 0, 1, and 2 respectively; healthy tissues representing SFreq = 3 could not be identified.
	- d. Slug: cardiac myocytes, smooth muscle of the appendix, and lymphocytes of the small intestine represented SInt 0, 1, and 2 respectively; healthy tissues representing SInt = 3 could not be identified.
	- e. HA: Intensity for hyaluronan was measured as blue color intensity developed by Alcian blue in comparison to hyaluronidase digested tissue section. Sub-mucosa of the small

intestine and cartilage tissues represented SInt 1 and 2 respectively; healthy tissues representing SInt = 0 or 3 could not be identified.

	- a. TCF21: cardiac myocytes, liver hepatocytes, germinal basal cells of testis represented SLoc 0, 1, and 2 respectively.
	- b. E-cadherin: cardiac myocytes, prostate epithelial cells represented SLoc 0 and 2 respectively; healthy tissues representing SLoc = 1 could not be identified.
	- c. PARP1: mucosa of the small intestine, germinal basal cells of testis represented SLoc 0 and 2 respectively; healthy tissues representing SLoc = 1 could not be identified.
	- d. Slug: cardiac myocytes, the somatic muscle of the appendix, lymphocytes of the small intestine, represented SLoc 0, 1, and 2 respectively.
	- e. HA: cartilage represented SLoc of score 2; healthy tissues representing SLoc = 1 could not be identified. A further consensus was reached in the pathology review to consider extracellular staining in tumor nests that is eliminated following hyaluronidase treatment as a proper localization, while distant stroma-associated HA was considered as mislocalization.
	- f. ANXA2: cardiac myocytes, stromal cells of the gall bladder, epithelial cells of the gall bladder represented SLoc as 0, 1, and 2 respectively.

#### *3.2. Establishment of Scoring Guidelines for Stratification Using a Panel of Xenograft*

Initial validation of the biomarker expression and scoring scheme was achieved using HGSC cell line derived xenografts generated in NOD/SCID mice (Figure 1C,D). TCF21 localization in xenograft sections was either dominantly nuclear (CAOV3 and PEO14), cytoplasmic (OVMZ6, OV90, and OVCAR3), or negligible (CP70 and A4); similarly, Slug was nuclear (OVMZ6, OV90, and A4), cytoplasmic (CAOV3 and CP70), or absent (OVCAR3 and PEO14). Moderate intensity of E-cadherin at the cell membrane was observed in ~50% tumor cells in CAOV3 and OVCAR3 xenografts but was lower in OVMZ6, CP70, OV90, A4, and PEO14 xenografts. Significantly, high-intensity expression of nuclear PARP1 was evident only in OVCAR3 xenografts; while other xenografts had significantly lower expression. High frequency, moderate intensity of HA was observed in CAOV3, OV90, and A4 xenografts, while OVCAR3 and PEO14 expressed HA at low to moderate frequency with weak intensity. Significant expression of ANXA2 was evident only in OVMZ6 and A4 xenografts. Consensus marker scores consolidated by the pathologist panel for each marker and xenograft (SFreq, SInt, SLoc; Table 1) were used to derive specific Biomarker Indices (BI; Equation (1)). Class-indices representing class-specific metrics of consolidated marker expression were derived from class-specific BI (CICCM and CIEMT; Equations (2) and (3) respectively; Table 1).

$$BI = \frac{1}{3} \left( \frac{\text{observed } S\_{Freq}}{\max \, S\_{Freq}} \right) + \frac{1}{3} \left( \frac{\text{observed } S\_{Int}}{\max \, S\_{Int}} \right) + \frac{1}{3} \left( \frac{\text{observed } S\_{Lac}}{\max \, S\_{Lac}} \right) \tag{1}$$

$$CI\_{\rm{CCM}} = \frac{BI\_{\rm{TCF21}} + BI\_{\rm{E-cadherin}} + BI\_{\rm{PARP1}}}{3} \tag{2}$$

$$CI\_{EMT} = \frac{BI\_{Slug} + BI\_{HA} + BI\_{ANXA2}}{3} \tag{3}$$

**Table 1.** Scores, biomarker and class indices (BI and CI respectively) for cooperative cell migration (CCM) and epithelial-to-mesenchymal transition (EMT) markers in xenografts.

Class indices represent class-specific metrics of consolidated markers expression; CICCM and CIEMT computed for xenografts ranged from (0–0.76) and (0.22–0.89) respectively (Table 1; Figure 1D). The distribution of median CIEMT vs. CICCM (±10%) values were further applied in class identification. Thus, OVCAR3 represents CCM-class; A4, OVMZ6, and OV90 the EMT-class; CAOV3, PEO14, and CP70 being double positive (DP; Figure 1E). Such inclusiveness of expression of the six markers quantifies molecular heterogeneity in mixed/unclassified tumors and assigns biological functions to the ambiguous Class 3 through relative marker expression.

#### *3.3. Evaluation of Stratification Guidelines in TMAs*

The above biomarker scoring and class identification guidelines were applied to commercial TMAs (duplicate cores per sample) which included two normal ovary and 13 HGSC tumor cases among other ovarian cancer subtypes. Availability of limited consecutive slides (*n* = 5) led to screening of only four of the six biomarkers (TCF21, E-cadherin, Slug, and HA) and equations (ii) and (iii) were appropriately modified in consideration of a four marker-based class identification (Supplementary Figure S4; Supplementary Table S2). Biomarker score averages for clinical cases represented on TMAs were considered for computation of BI and CI values (of each core in duplicate showed a near-similar expression for all biomarkers). We observed that a majority of HGSC TMA-cores expressed high-intensity cytoplasmic TCF21, moderate nuclear expression of which was present in ~5–10% of tumor cells. Except in four cases, Slug expression was weak to moderate cytoplasmic and nuclear localization was evident in only 5–10% tumor cells. Moderate expression of E-cadherin at cell membranes was observed in nearly 50% of tumor cells, and extracellular HA fibers also stained at a moderate intensity in tumor cell nests. Consolidation of biomarker scores of each case, computation of BI and CI values followed by plotting the distribution of CICCM vs. CIEMT values of TMA cores indicated three cases to represent CCM-Class while the remaining belonged to DP-Class; and EMT-Class remained unrepresented (Figure 1F).

#### *3.4. Evaluation of Clinical Samples Associates CCM-Markers with Metastases and Chemotherapy*

The variance in frequencies of class profiles between xenografts and TMAs emphasized the pertinence of screening larger numbers of clinical samples in validation. Towards assessing clinical representation, we obtained and stratified 160 tumor samples pathologically diagnosed as HGSC from 96 patients. These included primary (ovary (T), fallopian tube (FT)) and metastatic tumors (omentum (O), peritoneal ascites-derived cell blocks (A); Supplementary Table S3). BI and CI scores with CICCM and CIEMT values for each tumor were computed from marker scores (Supplementary Tables S4–S8) followed by evaluation of CI distribution towards class assignment as performed for xenografts and TMA. Data analysis of this clinical cohort was conducted in tumor groups as given Table 2.

Examination of 'Group A' tumors representing different sites of metastases and stages of tumor progression in six chemo-naïve cases revealed consistent class-associations across different sites in two cases (B/2774/12 and B/3136/09), while suggesting class switching in the remaining four cases wherein marker expression was altered in following metastases (Figure 2A). One of the latter four cases expressed CCM-markers in ovarian tumors; FT and omental tumors were DP-Class (B/1627/13). The ovarian tumor of case B/825/10 expressed EMT-markers that switched to DP-Class in FT and omental tumors. Ovarian and omental tumors of the remaining two cases segregated into DP-Class, while FT tumors expressed either CCM (B/749/13) or EMT (B/1716/09) markers. Overall, the three tumor sites predominantly segregated into CCM or DP Class; only one case of ovarian and FT tumors was represented as EMT-Class (Figure 2B; Supplementary Tables S9 and S10).

CI-based class-assignment in 'Group B' tumors stratified chemo-naïve ovarian tumors into DP, CCM or EMT classes (33.3%, 29.4%, and 29.4% respectively; Figure 2C,E), while treated ovarian tumors exhibited dominant representation of CCM-Class with lower frequencies of EMT and DP classes (48.0%, 23.0%, and 19.2% respectively). Chemo-naïve FT tumors stratified into CCM, EMT or DP Class (*n* = 2, 1, and 8 respectively; Figure 2C,E), while those after treatment belonged to either CCM or DP Class. Omental tumor deposits were predominantly DP-Class, with a marginal increase in the frequency of CCM and EMT markers in treated samples (Figure 2D,E), while chemo-naïve as well as treated ascites cell blocks presented more frequently with the CCM-class (Figure 2D,E). Overall, 'Group B' patient tumors show predominantly CCM or DP expression (Supplementary Table S11) as compared with either xenografts (mostly EMT-class) or TMAs (CCM and DP class). These findings support metastasis/chemotherapy-induced HGSC tumor expression towards a CCM subtype.

*J. Clin. Med.* **2019**, *8*, 330

**Figure 2.** (**A**) Scatter plot of CICCM vs. CIEMT distribution for chemo-naïve cases with tumors detected in ovary (Ov), fallopian tube (FT), and omentum (O) (left panel), and a reference case-chart (right panel), (**B**) graphical representation of class-specification of Group A tumors, (**C**) scatter plots of CICCM vs. CIEMT distribution for single chemo-naïve or -treated (red and blue shapes respectively) tumors from-ovary & fallopian tube -, (**D**) omentum ◦ and ascites cell blocks -, (**E**) graphical representation of Group B tumors (chemo-naïve-N; chemo-treated-T). EMT-epithelial to mesenchymal transition; DP-double positive; CCM-cooperative cell migration.

#### *3.5. HGSC Tumors at Different Sites Exhibit Molecular Heterogeneity and Class-Switching*

Class switching during tumor progression from ovarian to omental sites was further examined in 'Group C' samples that comprised of tumors from either chemo-naïve (*n* = 17) or -treated cases (*n =* 16). Almost half of the cases in both groups did indeed exhibit metastases associated class-switching (Figure 3(Ai,Aii); Supplementary Figure S5A), although lack of a specific direction to the switch possibly suggests the involvement of other factors in the determining marker expression.

The last analytical set of Group D tumors comprised of six ovarian, ascites, and/or omental tumor pairs before and after therapy (Figure 3B; Supplementary Figure S5B). Considering the limited cases that represent unique behavior, we have discussed them on a case to case basis. Cases 1 and 2 strongly conformed to CCM-class after disease progression as well as treatment, while the remaining four cases exhibited class-switching. Chemotherapy in Case 3 resulted in enhanced expression of CCM-markers over a DP-profile in untreated ovarian tumors, while Case 4 was associated with heterogeneity of marker expression following treatment. Case 5 was the most complex of all six and showed considerable marker heterogeneity between different tumor sites. Case 6 exhibited progression and therapy-associated class-shift towards the CCM class from a DP ovarian tumor. These findings further support the class-switching of HGSC tumors upon metastases and/or chemotherapy.

**Figure 3.** Class switching detected in cases with (**A**) ovarian and omental tumors tissues of (i) chemo-naïve (*n* = 17) and (ii) chemo-treated (*n* = 16) patient cohort; (**B**) class switching in tumors collected from patients (*n* = 6) prior to chemo-treatment (filled shapes) and post chemo-therapy (empty shapes) as determined through ascites (diamond), primary tumor (triangle), and omentum (circle). EMT-epithelial to mesenchymal transition; DP-double positive; DN-double negative; CCM-cooperative cell migration.

#### *3.6. Disease Progression is Inclined Towards Enrichment of CCM-Markers*

To evaluate the effects of therapy and disease progression vis-à-vis metastases and stage advancement on stratification, we further compared the means of CICCM and CIEMT (M-CICCM and M-CIEMT respectively) of ovarian, FT, and omental tumors from the same patient in chemo-naïve (CN) vs. chemo-treated (CT) groups. M-CICCM of CN ovarian-omental (T-O) tumors (Group 'C') was lower than that of CN ovarian-FT-omental (T-F-O) tumors (Group 'A'), while M-CICCM of CT T-O tumors was enhanced (Figure 4A, Supplementary Table S12). Interestingly, M-CIEMT of ovarian as well as omental tumors were similar between CN T-O vs. T-F-O tumors, and between CN T-O vs. CT T-O groups suggesting minimal effects of either disease progression or therapy on expression of EMT-markers. Towards elucidation of individual biomarker contribution likely to contribute to these differences, we analyzed their expression through means of BI (M-BI) during stage progression in CN HGSC cases represented on TMAs (T1 to T2) and in patient tumors. Nearly steady M-BI of TCF21, increased Slug, E-cadherin, and HA levels were revealed in the TMAs, with Slug expression being maximal at the T2 stage (*p* < 0.007; Figure 4B). Similar analyses of Group 'B' cases with available tumor stage information in CN (31-T, 1-FT, 17-O) and CT cases (35-T, 7-FT, and 9-O) were performed. CN T samples expressed comparable M-BI levels of TCF21 and HA, higher E-cadherin and PARP1 and lower Slug and ANXA2 at T3 over T1 stage; CT T samples had almost comparable M-BI profiles as CN tumors (Figure 4C, Supplementary Figure S6). Limited FT and O tumors at stages T1 and T2 restricted their analyses; CN- and CT-O tumors at stage T3 expressed comparable levels of markers. Group C CN O samples at stage T3 were associated with decreased TCF21 and E-cadherin concurrent with increased PARP1,

ANXA2, and HA M-BI profiles (Figure 4D). CT T tumors had low M-BI scores for TCF21, PARP1, and HA in comparison to O and similar levels of E-cadherin, Slug, and ANXA2. CT T samples in this group displayed steady levels of TCF21 and PARP1, lower E-cadherin, Slug, and ANXA2 and increased HA levels over CN tumors; while treatment was associated with higher BI levels of TCF21, PARP1, and HA and lower levels of E-cadherin, Slug, and ANXA2 in the O tumors. This suggests that selected biomarkers play a significant role during disease progression.

**Figure 4.** HGSC progression-associated marker expression. (**A**) Plot comparing CICCM with CIEMT of tumors of Groups B and C respectively; (**B**) plot comparing biomarker (BI) scores for TCF21, E-cadherin, Slug, and HA in chemo-naïve ovarian tumors present in TMA for stages T1 and T2; (**C**) plot for chemo-naïve (left) and treated (right) paired ovarian (T)-omental (O) tumors; (**D**) plot for chemo-naïve (CN) and treated (CT) tumors at stages T1, T2, and T3 in ovarian and omental ◦ tumors.

Class-switching in paired samples led us to examine similar effects of a chemotherapeutic challenge (paclitaxel) in the three classes representing cell lines derived xenografts viz. CCM-Class (OVCAR3), EMT-Class (A4, OVMZ6), and DP-Class (CAOV3, PEO14). Distribution of CICCM and CIEMT values revealed a DP to CCM (CAOV3, PEO14) and EMT to DP (A4, OVMZ6) class switch. An outlier was the OVCAR3 xenograft (CCM class) that despite an increased CI score, provided no evidence of class switching (Figure 5A,B).

**Figure 5.** (**A**). Paclitaxel exposure alters scoring marker panel in HGSC cell line derived xenografts. Representative images of HGSC xenograft (control and paclitaxel treated) sections for Row 1-HE (hematoxylin and eosin), Rows 2,3,4,7,8 represent IHC-based detection of TCF21, Slug, E cadherin, PARP1, and ANXA2, Rows 5 and 6 represent HC-based identification of HA fibers in untreated and hyaluronidase digested sections respectively, scale bars-50 μm; (**B**) scatter plots of CICCM vs. CIEMT derived from xenograft (control-grey and paclitaxel treated-red) scoring. EMT-epithelial to mesenchymal transition; DP-double positive; CCM-cooperative cell migration.

#### *3.7. Correlation Between Transcript- and Protein-Based Stratification*

Proteome-based profiling of TCGA samples has been reported earlier (Figure 6A) [26,27]. Comparison of tumor samples common to our earlier study and the five proteomic subtypes (*n* = 61; Refs. 24 vs. 26) indicated correlation between CCM and proliferative groups, while EMT tumors were dominantly mesenchymal, with a few being either immunoreactive or differentiated (Figure 6B); surprisingly, the stromal class had negligible associations. A similar comparison of tumor samples common to the transcript study and two proteomic subtypes (*n* = 34) [24,27] correlated the CCM-Class with TCGA-A/Epithelial cluster and EMT-Class with TCGA-B/mesenchymal cluster (Figure 6C). These observations indicate some degree of variation that could arise from differences between transcriptomic and protein abundances.

**Figure 6.** Comparison of HGSC stratification approaches. (**A**) Proteomics-based TCGA HGSC tumors stratification in 5 sub-groups by Zhang et al. 2016 (169 cases) and Coscia et al. 2016 (84 cases) and transcriptomics based HGSC stratification into three classes (Gardi et al. 2014); (**B**) comparison of Gardi et al. 2014 vs. Zhang et al., 2016; (**C**) comparison of Gardi et al. 2014 vs. Coscia et al. 2016.

#### **4. Discussion**