Quantitative Analysis of Differential Proteome Expression in Bladder Cancer vs. Normal Bladder Cells Using SILAC Method

The best way to increase patient survival rate is to identify patients who are likely to progress to muscle-invasive or metastatic disease upfront and treat them more aggressively. The human cell lines HCV29 (normal bladder epithelia), KK47 (low grade nonmuscle invasive bladder cancer, NMIBC), and YTS1 (metastatic bladder cancer) have been widely used in studies of molecular mechanisms and cell signaling during bladder cancer (BC) progression. However, little attention has been paid to global quantitative proteome analysis of these three cell lines. We labeled HCV29, KK47, and YTS1 cells by the SILAC method using three stable isotopes each of arginine and lysine. Labeled proteins were analyzed by 2D ultrahigh-resolution liquid chromatography LTQ Orbitrap mass spectrometry. Among 3721 unique identified and annotated proteins in KK47 and YTS1 cells, 36 were significantly upregulated and 74 were significantly downregulated with >95% confidence. Differential expression of these proteins was confirmed by western blotting, quantitative RT-PCR, and cell staining with specific antibodies. Gene ontology (GO) term and pathway analysis indicated that the differentially regulated proteins were involved in DNA replication and molecular transport, cell growth and proliferation, cellular movement, immune cell trafficking, and cell death and survival. These proteins and the advanced proteome techniques described here will be useful for further elucidation of molecular mechanisms in BC and other types of cancer.


Introduction
Bladder cancer (BC) is the fifth most common type of human cancer. There were an estimated 74,690 newly diagnosed cases and 15,580 deaths from this disease in the United States in 2013 [1]. Of total BC patients, >70% have nonmuscle-invasive disease and~25% present initially with muscle invasion. Patients with the muscle-invasive form have a 50% risk of distant metastases and a poor prognosis [2]. The recurrence of superficial bladder tumors is a major reason for the worldwide prevalence of BC [3]. The majority (90%) of BCs are classified histologically as urothelial carcinomas (UCs), derived from the bladder urothelium [4]. Bladder epithelial tissues have a clear hierarchical organization consisting of three morphologically distinct cell types: basal, intermediate, and umbrella cells, corresponding respectively to early, middle, and late differentiation states [5]. Malignant transformation may occur in each of these cell types, resulting in a diversity of tumor phenotypes [6]. According to the latest report of the American Cancer Society, the relative 5-year survival rate for BC with early detection (stage I, (T1, N0, M0)) is~88% [7]. Therefore, identification of novel early-stage molecular markers is desirable for improved risk stratification.
Candidate biomarkers for BC detection evaluated to date include telomerase, bladder tumor antigen (BTA), nuclear matrix protein 22 , and fibrin degradation product (FDP). The reliability of tests based on these biomarkers is poor because of low sensitivity and high false-positive rates [8][9][10][11]. Proteins can potentially be identified specific to aggressive or nonaggressive types of cancer. Proteome analysis is challenging because of the limited amount of available clinical sample [12]. Monitoring of the proteome of BC cells could provide additional information for clinical diagnostic purposes.
Recent advances in mass spectrometry (MS) for protein identification and quantification facilitate in-depth analysis of large numbers of proteins, and have been used for examination of the whole proteome in several systems. Such methods include 2D difference gel electrophoresis (2D DIGE), the similar iTRAQ (isobaric tag for relative and absolute quantitation), isotopecoded affinity tagging (ICAT), and stable isotope labeling by amino acids in cell culture (SILAC) [13][14][15]. In comparison with peptide-based absolute quantitation methods, SILAC has the advantages of mixing samples at the very beginning, and reduced sample-to-sample variability. Metabolic labeling with stable isotopes has been described as the "gold standard" in protein quantification [16]. Arginine (Arg) and lysine (Lys) are the stable isotope-labeled amino acids most frequently used in SILAC-based studies, because subsequent trypsin digestion of isolated proteins (which cleaves at basic residues) for MS analysis generates peptides with a single labeled amino acid, simplifying analysis and quantification [17]. In the present study, three stable isotopes each of Arg (R0, R6, R10) and Lys (K0, K4, K8)in three separate cultures ("light" (L), "medium" (M), and "heavy" (H)) were used to analyze proteome differences at various stages of BC. Distinctive L, M, and H forms of each peptide as detected by MS reflected relative amounts of the corresponding protein in three isotopically encoded BC cell stages.
Three human cell lines were studied: normal bladder epithelial HCV29, low grade nonmuscle invasive bladder cancer (NMIBC) KK47, and metastatic muscle invasive bladder cancer YTS1. Each of the three cell lines was cultured in media added with three combinations of unlabeled Lys and Arg ("light"), D 4 -Lys and 13 C 6 -Arg ("medium"), and 13 C 6 15 N 2 -Lys and 13 C 6 15 N 4 -Arg ("heavy"). Proteins with >98% label incorporation were analyzed and quantified by 2D-HPLC-LTQ Orbitrap MS (Fig 1). Differential expression of the identified proteins, which  Total peptides were concentrated and desalted using a 10KD size-exclusion spin ultrafiltration unit and dried using a SpeedVac concentrator.

Data analysis
Raw MS data were analyzed using the MaxQuant software program (V. 1.2.2.5) [26,27]. A false discovery rate (FDR) of 0.01 for proteins and peptides and a minimum peptide length of 6 amino acids were required. MS/MS spectra were searched by Andromeda [28] against the IPI human database (V. 3.85). The MaxQuant program determined the SILAC state of peptides from mass differences between SILAC peptide pairs, and this information was used to perform searches with fixed Arg6 and Lys4 or Arg10 and Lys8 modifications, as appropriate. Quantification in MaxQuant was performed as described previously [26]. Differential regulation within each experimental M/L ("medium/ light") ratio and H/L ("heavy/ light") ratio of the identified proteins was normalized using z-score analysis, as described previously [29,30]. In brief, M/L and H/L ratios were converted into log2 space, and average ratios and SD (standard deviations) were calculated for each data set. The log2 M/L and H/L ratio of each protein were converted into a z-score, using the following formula: z À scoresðsÞ of ½b ¼ log 2 XðM or HÞ L ½b À Average of ðlog 2 of each number; a . . . nÞ Standard deviation of ðlog 2 of each number; a . . . nÞ where b were deemed as a single protein in a data set population (a. . ..n). The z-score was a measure of how many SD units (σ) of the log2 M/L or H/L ratio of the protein was away from the population mean. A z-score !1.960σ represented that differential expression of the protein lied outside the 95% confidence interval, a score !2.576σ represented expression outside the 99% confidence interval, and a score !3.291σ represented 99.9% confidence. Z-scores !1.960σ were considered to be significant [29].

Functional annotation and Ingenuity Pathways Analysis
Identified proteins were further analyzed using the SWISS-PROT database to classify their biological process, cellular component, and molecular function [31]. Significant over-represented gene ontology (GO) terms were identified using the Database for Annotation, Visualization and Integrated Discovery (DAVID) gene bioinformatic resources [32,33]. Proteins determined to be differentially regulated as described in the preceding section were tabulated in Excel and their International Protein Index (IPI) numbers were uploaded into DAVID (http://david.abcc.ncifcrf. gov/home.jsp) for functional annotation analysis. Data sets containing gene identifiers and corresponding expression values were then uploaded into the Ingenuity Systems application. Each IPI number was mapped to its corresponding gene object in the Ingenuity Pathways Knowledge Base. Networks of the proteins were generated algorithmically based on their connectivity. Fisher's exact test was used to calculate a p-value indicating the probability that a particular biological function and/or disease assigned to that network was due to chance alone.

Quantitative real-time PCR
Cells (1×10 5 per well in a 6-well plate) were cultured and treated as described above. Total RNA was isolated using an RNApure Tissue Kit (CWBiotech; Beijing, China) according to the manufacturer's instructions. Primers were designed using the DNAMAN program (V. 6.0.3; Lynnon Biosoft, Canada). First-strand cDNA was synthesized from total RNA using ReverTra Ace-α (Toyobo; Osaka, Japan). Quantitative real-time PCR was performed by LightCyclerbased SYBR Green I dye detection with UltraSYBR Mixture (CWBiotech). Gene expression was quantified by the 2 -ΔΔCT method [35].

Cell staining
Cells were cultured on sterilized coverslips in 24-well plates until 70-80% confluence, washed, immobilized, permeabilized with 0.2% Triton X-100 for 10 min at room temperature, and blocked with 5% nonfat milk overnight at 4°C. Fixed cells were incubated with diluted primary antibodies for 12 hr, incubated with FITC-labeled secondary antibodies at 4°C for 6 hr in the dark, washed, stained with 4 μg/mL DAPI at room temp for 10 min, washed with 1×PBS, and photographed with a fluorescence microscope (Eclipse E600, Nikon; Tokyo, Japan).

Determination of isotope incorporation efficiency
To analyze dynamic changes in BC oncogenesis at the proteome level, the SILAC method was applied to three bladder cell lines to obtain labeled cell populations. Sufficient labeling is a prerequisite for reliable quantification using this method. In the case of incomplete labeling of proteins, particularly for labeling efficiency <95%, quantitation of low-abundance proteins would be masked by contaminated "light" peptides. To determine incorporation efficiency of labeled Lys and Arg, "light" (HCV29), "medium" (KK47), and "heavy" (YTS1) proteins were separated, and the high-abundance protein was in-gel digested. MALDI-TOF/TOF-MS results for peptide GVVDSEDLPLNISR in heat shock protein 90 (P08238) indicated that complete incorporation of isotopically labeled Arg and Lys was achieved in KK47 and YTS1, and no Arg-to-Pro conversion occurred (Fig 2A). LC-ESI-MS/MS analysis of the doubly charged peptide VNQIGSVTESLQACK of alpha-enolase (P06733) and GGPEVQQVPAGER of fatty acid synthase (P49327) showed that these doublets of actual peak clusters were from HCV29, KK47, and YTS1 cells, respectively ( Fig 2B).

SILAC cell model for quantification of proteome in BC progression
Proteins isolated from the three cell lines were mixed (1:1:1) and digested using a 10 KD filter (Millipore; Billerica, MA, USA). Peptides were analyzed by ultrahigh-resolution liquid chromatography-tandem MS (nLC-ESI-MS/MS) on a hybrid linear ion trap LTQ Orbitrap instrument. A total of 3721 unique proteins were identified in two independent replicate experiments (Fig 3A and S1 Table). Of these, 1766 proteins (47.5% of the total) that were identified in both experiments and satisfied the criteria established for protein quantitation were subjected to further bioinformatic analysis. The distribution histograms of log ratios for both M/L and H/L fit a Gaussian distribution. Most of the identified proteins were within the ±1 range of log ratios (Fig 3C and 3D). Using 1 as the threshold log ratio, expression of 255 proteins was higher in both KK47 and HCV29, and expression of 434 proteins was lower in both KK47 and HCV29 Quantitative Proteome Analysis of Bladder Cells ( Fig 3B). Population distribution-based z-scores allowed direct comparison of proteins from different experiments. Differing confidence level cutoffs were applied to the data by z-score analysis to determine which proteins were significantly differentially regulated. The cutoffs applied were 95%, 99% and 99.9%, corresponding to z-scores of ±1.960, ±2.576, and ±3.291, respectively. Using a 95% cutoff, significant differential regulation was observed for 110 proteins in KK47 vs. HCV29 (36 upregulated, 74 downregulated) and for 87 proteins in YTS1 vs. HCV29 (17 up, 70 down). Differential regulation was observed for 35 proteins in the two BC lines vs. HCV29 using a 99% cutoff (2 up, 33 down), but for only five of these proteins using a 99.9% cutoff (2 up, 3 down) ( Table 1). Tables 2 and 3 list the upregulated and downregulated proteins determined in the two experiments, with their average SILAC ratios and z-scores. All proteins differentially regulated with >95% confidence had a >5-fold alteration of SILAC ratios, and most proteins differentially regulated with >99% confidence had a >10-fold alteration of SILAC ratios.
Hierarchical clustering analysis of samples was performed to examine correlations of proteome patterns among the three cell lines. The Cluster graph ("heat map") shows that samples of the same cell line cluster together ( Fig 3E). Some proteomes were distinctive among the three cell lines based on significant alterations in metastatic vs. low grade nonmuscle invasive, whereas other proteomes were moderate and consistent. Proteins in the former group may be involved in BC development.

Functional classification and pathway analysis of identified proteins
Functional interpretation is a crucial step in data analysis when extensive functional annotation of the data sets is not available. Taking into account their nonexclusive localization in GO, the identified proteins were linked to at least one annotation term each within the GO molecular function, biological process, and molecular component categories. The most common molecular functions were binding (47.2%), and catalytic activity (30.9%) (Fig 4A). The major biological process categories were cellular (16.3%), single-organism (14.2%), and metabolic (13.8%) ( Fig 4B). The major cellular component categories were cell (17.6%), cell part (17.6%), and organelle (15.8%) (Fig 4C).
To identify enrichment terms associated with the upregulated and downregulated groups of proteins after averaging of z-scores using the 95% cutoff, lists of proteins were uploaded to the DAVID website using the complete human proteome as background. To help clarify which molecular functions and biological processes were most affected during BC maturation, overrepresented GO terms were identified based on threshold count ! 2 and Expression Analysis Systematic Explorer (EASE) < 0.1. The over-represented molecular functions, biological processes, and cellular components in the significant enriched GO terms of upregulated proteins were analyzed. The most highly ranked molecular function was neutral amino acid transmembrane transporter activity (2 proteins). The most highly ranked biological processes were cellular amino acid derivative metabolic process (4 proteins), ribonucleoprotein complex biogenesis (4 proteins), response to extracellular stimulus (4 proteins), and nitrogen compound biosynthetic process (4 proteins). The most highly ranked cellular components were intracellular non-membrane-bounded organelle (14 proteins) and non-membrane-bounded organelle (14 proteins).
Next, over-represented molecular functions, biological processes, and cellular components in the significant enriched GO terms of downregulated proteins were analyzed. The most highly ranked molecular functions were calcium ion binding (16 proteins), structural molecule activity (11 proteins), and identical protein binding (10 proteins). The most highly ranked biological processes were cell adhesion (14 proteins), biological adhesion (14 proteins), response to wounding (13 proteins), and immune response (13 proteins). The most highly ranked cellular components were plasma membrane (40 proteins), plasma membrane part (30 proteins), and non-membrane-bounded organelle (24 proteins) (S2 Table).
Proteins were further analyzed, and metabolic and canonical pathways and interconnecting proteins were generated, using Ingenuity Pathways Analysis (IngenuityH Systems, www. ingenuity.com). The top network functions identified as upregulated proteins in BC cells were involved in DNA replication, amino acid metabolism, molecular transport (52 proteins; Fig  5A), gene expression and hereditary disorders (33 proteins), cell growth and proliferation (20 proteins; Fig 5B), and post-translational modification and cancer (4 proteins). The top network functions identified as downregulated proteins in BC cells were cellular movement and immune cell trafficking (97 proteins; Fig 5C), lipid metabolism (31 proteins; Fig 5D), cellular development, growth and proliferation, and cell death and survival (6 proteins). These findings indicate that BC cell proteomes were continuously shifting depending on the stage of cell metastasis.

Confirmation of MS results by western blotting
Variations of the differential proteins described above were confirmed by western blotting. THY1, MAGEA4, IGF2BP1 and SFN were detected at higher levels in BC KK47 and YTS1 cells than in normal bladder epithelial HCV29 cells, whereas VIM, CTNNB1, FN1 and CD70 were detected at lower levels in KK47 and YTS1 than in HCV29 (Fig 6B and 6C). In general, the western blotting results were consistent with the variables from MS analysis (Fig 6A; S1

Confirmation of SILAC results by qRT-PCR
The expression of six responding genes at the transcriptional level was evaluated by qRT-PCR. In BC KK47 and YTS1 cells, expression of MAGEA4, THY1, and IGF2BP1 was significantly increased, whereas that of VIM, CTNNB1, and FN1 was greatly reduced (Fig 6D). These findings are consistent with SILAC results.

Confirmation of SILAC and western blotting results by cell staining with antibodies
Antibodies recognizing MAGEA4, THY1, IGF2BP1, VIM, CTNNB1, and FN1 proteins, whose expression differed significantly in BC cells vs. HCV29 cells, were used to confirm the results of previous analyses and to assess protein distributions. Fluorescence signal intensities in the two BC cell lines were significantly higher for MAGEA4, THY1, and IGF2BP1 and significantly lower for VIM, CTNNB1, and FN1, in agreement with SILAC, western blotting, and RT-PCR results. Preferential localization was observed for MAGEA4 in central cytoplasm (including mitochondria and centrosome) and the nuclear region, for THY1 and VIM in cytoplasmic membrane and cytoplasm, for IGF2BP1 in the nuclear region and cytoplasm, and for CTNNB1 and FN1 in cytoplasmic membrane (Fig 7).

Discussion
Urothelial carcinoma of the bladder is unique among epithelial carcinomas in its divergent pathways of tumorigenesis. At the time of diagnosis for transitional cell BC,~80% of patients are in the low-grade (grade 1-2) non-muscle invasive (cTa or cT1) stage. In patients with low grade Ta disease, the 15-year progression-free survival is 95% with no cancer-specific mortality [36]. Therefore, biomarkers are needed for predicting a risk of stage progression from noninvasive to invasive, or non-metastatic to metastatic and predicting responsiveness to systemic therapies. The human cell lines HCV29 (normal bladder epithelia), KK47 (low grade nonmuscle invasive bladder cancer), and YTS1 (metastatic bladder cancer) have been widely used in studies of molecular mechanisms and cell signaling during the progression of bladder cancer to muscle or metastatic states [37]. However, little attention has been paid to global quantitative proteome analysis of these three cell lines. Quantitative analysis of proteins at different stages of BC progression is a challenging but important task for understanding of disease mechanisms. SILAC, a differential isotope labeling strategy that involves metabolic labeling of proteins in vivo, has been widely applied in cell biology and studies of model organisms such as yeast, bacteria, nematodes, plants, and mice [38,39]. In SILAC, natural "light" isotopes of carbon, nitrogen, and hydrogen in amino acids incorporated during protein translation are substituted with "heavy" isotopes such as 13 C, 15 N, and 2 H. No study to date has quantified differences in protein abundance at various stages of BC. Previous studies have focused on comparative protein levels in BC patients vs. control subjects, but not on alterations of protein levels during BC progression [8,9].
Identification and characterization of protein levels at various steps of differentiation are essential for our understanding of normal tissue development and malignant transformation. In the present study, we applied the SILAC method to analysis of HCV29, KK47, and YTS1. This strategy provided protein information for all three cell lines during one MS experiment. Following normalization by the z-score method, differential regulation was observed for 110 proteins in KK47 vs. HCV29 and for 87 proteins in YTS1 vs. HCV29. These differentially regulated proteins, which may play important roles in BC development, include SFN (14-3-3 protein sigma), SELENBP1 (selenium-binding protein 1), COL6A3 (collagen α3 (VI) chain), and CD70. COL6A2 and COL6A3 protein levels are reduced in urine of BC patients [40]. SFN is downregulated in invasive bladder transitional cell carcinomas undergoing epithelial-to-mesenchymal conversion and highly upregulated in pure squamous cell carcinomas [41]. Some of the proteins that are altered in BC are also altered in other cancers. In particular, expression of FOLR1 (folate receptor alpha) is increased relative to normal tissue in many types of epithelial cancer, including non-mucinous ovarian, endometrial, non-small cell lung, colorectal, and breast cancers [42].
The DAVID bioinformatic resources provide a set of powerful tools to explore their large gene lists in depth from many different biological perspectives in order to fully extract associated biological meanings. Enrichment analysis of GO terms in our identified proteins indicated that neutral amino acid transmembrane transporter activity and metabolic processes were significantly upregulated and that protein binding, cell adhesion, biological adhesion, and immune response were significantly downregulated in the BC cell lines. The observed changes in GO terms tended to promote cellular proliferation, tumor development, cancer cell progression and metastasis, and escape from immune system surveillance.
A major challenge in cancer biology is the formulation of biological hypotheses regarding biomarker candidates [43]. Pathway analysis of differentially expressed proteins indicates that DNA replication and molecular transport, cell growth, and cell proliferation are requirements for cancer cell metastasis. Intercellular adhesion molecule-1 (iCAM-1) is a transmembrane glycoprotein present at basal levels in a wide variety of cell types and is upregulated in response to a number of inflammatory mediators [44]. The biological significance of iCAM-1 expression in cancer remains controversial; it is elevated in gastric, breast, oral, and thyroid cancer tissues [44][45][46][47] but reduced in some ovarian adenocarcinoma cell lines and primary tumors [48]. Also, iCAM-1 expression is up-regulated in squamous cell types associated with inflammation, such as schistosomal bladder cancer [49]. However, in the present study, iCAM-1 expression was reduced in metastatic bladder cancer cells. Moreover, CIT, a novel tissue-specific Ser/Thr kinase that encompasses the Rho-Rac-binding protein Citron, plays a role in cytokinesis and in Proteins were transferred to a PVDF membrane, probed with their primary antibodies, and incubated with HRP-conjugated rabbit anti-mouse or goat anti-rabbit secondary antibodies. (C) Densitometric quantitation of the protein levels. The protein was normalized by the tublin, and then compared to HCV29, which were arbitrarily set at 1.0. (D) Gene expression for the proteins was analyzed by qRT-PCR. Relative expression in comparison to control samples was analyzed by the 2 −ΔΔCt method and represented as Log 2 . Expression of genes above Log 2 (2) or below Log 2 (1/2) was significantly upregulated or downregulated, respectively.  Rho signaling that modulates myosin phosphorylation and cell adhesion [50]. Upregulation of CIT promotes DNA replication and cancer cell proliferation.
In conclusion, we successfully applied the SILAC method to identify and quantify proteins whose level is significantly up-or downregulated during BC development or progression. Similar advanced proteome techniques will be useful for further elucidation of biomarkers and molecular mechanisms in BC and other types of cancer.  Table. Gene ontology statistics of significantly upregulated and downregulated proteins. (XLSX)