Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

A Unique Four-Hub Protein Cluster Associates to Glioblastoma Progression

  • Pasquale Simeone,

    Affiliation Unit of Cancer Pathology, Ce.S.I., Foundation University “G. d'Annunzio,” Chieti, Italy

  • Marco Trerotola,

    Affiliation Unit of Cancer Pathology, Ce.S.I., Foundation University “G. d'Annunzio,” Chieti, Italy

  • Andrea Urbanella,

    Affiliation Unit of Cancer Pathology, Ce.S.I., Foundation University “G. d'Annunzio,” Chieti, Italy

  • Rossano Lattanzio,

    Affiliations Unit of Cancer Pathology, Ce.S.I., Foundation University “G. d'Annunzio,” Chieti, Italy, Department of Experimental and Clinical Sciences, School of Medicine and Health Science, University “G. d'Annunzio,” Chieti, Italy

  • Domenico Ciavardelli,

    Affiliations School of Human and Social Science, University “Kore” of Enna, Enna, Italy, Molecular Neurology Unit, Ce.S.I., University “G. d'Annunzio,” Chieti, Italy

  • Fabrizio Di Giuseppe,

    Affiliations Aging Research Center, Ce.S.I., University “G. d'Annunzio” Foundation, Chieti, Italy, Department of Experimental and Clinical Sciences, School of Medicine and Health Science, University “G. d'Annunzio,” Chieti, Italy, StemTeCh Group, Chieti, Italy

  • Enrica Eleuterio,

    Affiliations Aging Research Center, Ce.S.I., University “G. d'Annunzio” Foundation, Chieti, Italy, Department of Experimental and Clinical Sciences, School of Medicine and Health Science, University “G. d'Annunzio,” Chieti, Italy, StemTeCh Group, Chieti, Italy

  • Marilisa Sulpizio,

    Affiliations Aging Research Center, Ce.S.I., University “G. d'Annunzio” Foundation, Chieti, Italy, Department of Experimental and Clinical Sciences, School of Medicine and Health Science, University “G. d'Annunzio,” Chieti, Italy, StemTeCh Group, Chieti, Italy

  • Vincenzo Eusebi,

    Affiliation Department of “Tutela Salute Donna, Vita nascente, Bambino e Adolescente,” Catholic University of the Sacred Heart, Policlinico Universitario “Agostino Gemelli,” Roma, Italy

  • Annalisa Pession,

    Affiliation Section of Surgical Pathology, “M. Malpighi,” Bellaria Hospital, Bologna, Italy

  • Mauro Piantelli,

    Affiliations Unit of Cancer Pathology, Ce.S.I., Foundation University “G. d'Annunzio,” Chieti, Italy, Department of Experimental and Clinical Sciences, School of Medicine and Health Science, University “G. d'Annunzio,” Chieti, Italy

  • Saverio Alberti

    s.alberti@unich.it

    Affiliations Unit of Cancer Pathology, Ce.S.I., Foundation University “G. d'Annunzio,” Chieti, Italy, Department of Neuroscience, Imaging and Clinical Sciences, University “G. d'Annunzio,” Chieti, Italy

Abstract

Gliomas are the most frequent brain tumors. Among them, glioblastomas are malignant and largely resistant to available treatments. Histopathology is the gold standard for classification and grading of brain tumors. However, brain tumor heterogeneity is remarkable and histopathology procedures for glioma classification remain unsatisfactory for predicting disease course as well as response to treatment. Proteins that tightly associate with cancer differentiation and progression, can bear important prognostic information. Here, we describe the identification of protein clusters differentially expressed in high-grade versus low-grade gliomas. Tissue samples from 25 high-grade tumors, 10 low-grade tumors and 5 normal brain cortices were analyzed by 2D-PAGE and proteomic profiling by mass spectrometry. This led to identify 48 differentially expressed protein markers between tumors and normal samples. Protein clustering by multivariate analyses (PCA and PLS-DA) provided discrimination between pathological samples to an unprecedented extent, and revealed a unique network of deranged proteins. We discovered a novel glioblastoma control module centered on four major network hubs: Huntingtin, HNF4α, c-Myc and 14-3-3ζ. Immunohistochemistry, western blotting and unbiased proteome-wide meta-analysis revealed altered expression of this glioblastoma control module in human glioma samples as compared with normal controls. Moreover, the four-hub network was found to cross-talk with both p53 and EGFR pathways. In summary, the findings of this study indicate the existence of a unifying signaling module controlling glioblastoma pathogenesis and malignant progression, and suggest novel targets for development of diagnostic and therapeutic procedures.

Introduction

In 2013, more than 23,000 individuals were expected to be diagnosed with primary tumors of brain and central nervous system and more than 14,000 deaths were expected in the US alone [1]. The World Health Organization defines pilocytic (Grade I) and diffuse (Grade II) astrocytomas as low-grade brain tumors; anaplastic astrocytomas (Grade III) and glioblastomas (Grade IV; also designated as glioblastoma multiforme, GBM) are high-grade malignant tumors [2],[3]. With an annual incidence of 2–3 per 100,000 in Europe and US, GBM is the most frequent and aggressive form of brain tumor (60–70% of total malignant gliomas), and is essentially incurable [3], [4]. GBM consists of poorly differentiated, highly invasive neoplastic astrocytes; histopathological features include cellular polymorphism, nuclear atypia, mitotic activity, vascular thrombosis, microvascular proliferation and necrosis [5]. Regional heterogeneity of GBM frequently causes diagnostic discrepancies (≥20% of cases). Moreover, a high percentage of gliomas, such as mixed oligoastrocytomas and lower-grade gliomas, remain difficult to categorize reproducibly due to considerable histological overlap. These factors can compromise choice as well as effectiveness of therapeutic options [6]. Histopathologic diagnosis can be further compromised when only small biopsies are available. Additional molecular markers are thus urgently needed to efficiently discriminate among patients with distinct outcomes.

Loss of PTEN, amplification of EGFR and alterations of TP53, PDGFRA and CDKN2A/P16 are frequently found to be associated with GBM pathogenesis [5], [7]. Primary GBMs develop de novo after a short clinical history and without evidence of precursor lesions, whereas “secondary” GBMs arise from pre-existing diffuse or anaplastic astrocytomas. The signaling pathways responsible for development and growth of primary versus secondary GBM appeared as profoundly diverse, suggesting these two types of GBM to be different disease entities. Rather diverse genetic signatures were further proposed in the attempt of explaining GBM pathogenesis and heterogeneity [5], [8][12]. However, the actual impact of genetic signatures for GBM diagnosis and prognosis remains to be defined.

Genomic and transcriptomic data have provided key insight in GBM pathophysiology [11][14]. Corresponding insight into GBM proteomics [15] has not yet been achieved. Proteomic analysis of low- and high-grade tumors has tried to fill this gap [16][20], but analysis of specific protein markers has largely failed to provide a comprehensive view of GBM pathology.

In this study, we set to identify significantly modulated protein clusters that may bear functional impact and robustly explain distinct, relevant GBM pathology components. Proteomic analysis of human high-grade tumors, low-grade tumors and control tissue samples from normal brain cortex was thus systematically intersected through multivariate statistical procedures (principal component analysis, PCA and partial least square-discriminant analysis, PLS-DA). Using this approach, we were able to identify protein clusters discriminating tumors from normal tissues as well as high-grade from low-grade gliomas. Connectivity network analysis then allowed to discover a GBM control module that encompassed four major signaling hubs centered on Huntingtin, HNF4α, 14-3-3ζ and c-Myc. This proteomic signature was shown to underlie p53 and EGFR signaling, as an interconnected network. The GBM control module is candidate to be used as diagnostic biomarker and as target for therapeutic intervention. It may also help drafting a unifying model for glioblastoma appearance.

Materials and Methods

Patients and tissue specimens

Bioptic samples from low-grade and high-grade glioma patients were frozen in liquid nitrogen and stored at −80°C at the Section of Pathology “M. Malpighi” of the Bellaria Hospital, University of Bologna, between 1990 and 2002. Corresponding formalin-fixed paraffin embedded (FFPE) samples were stained with haematoxylin-eosin for routine histological diagnosis. The protocol of this study was approved by the board of the Ministry of the University and Research (“Novel technologies for glioblastoma assessment”, FISR Neurobiotechnologies, Grant N 481). Informed consent was previously obtained as indicated in Marucci et al. [21].

A total of 10 low-grade glial tumors (4 oligodendrogliomas, OL, 4 pilocytic astrocytomas, PA and 2 fibrillary astrocytoma, FA) and 25 GBMs were collected. All samples were re-staged and graded by expert pathologists according to 2007 WHO central nervous system tumor classification [2]. Control samples were tissues from five normal cortices from different brain regions, as obtained at autopsy from individuals deceased from diseases not involving the brain.

Reagents and chemicals

All reagents and chemicals are purchased from Sigma-Adrich (St. Louis, MO, USA), Bio-Rad Laboratories (Hercules, CA, USA) and GE Healthcare (Little Chalfont, UK).

Brain tumor lysates

Frozen brain tumor specimens were thawed on ice and resuspended in 2-DE lysis buffer (8 M urea, 40 mM Tris base, 65 mM DTT). Tumor lysates were briefly sonicated in Eppendorf (Hamburg, Germany) tubes with three 10-sec bursts, in 4% CHAPS. The lysates were centrifuged for 15 min at 12000 rpm to remove cell debris. Lysate supernatants were then processed for 2D PAGE analysis (Two dimensional polyacrylamide gel electrophoresis).

2D Electrophoresis

The first dimension was run over non-linear immobilized pH gradients (3.5–10.0 NL IPG 18 cm) (Pharmacia-Hoeffer Biotechnology AB, CA, USA). Hydration was achieved overnight in the reswelling cassette with 25 ml of a solution containing 8 M urea, 2% CHAPS (w/v), 10 mM DTE, 2% (v/v) pH 3.5–10 Ampholites, bromophenol blue and 200 µg of protein extract [22]. Run strips were equilibrated in 50 mM Tris-HCl pH 8.4, 6 M urea, 2% (w/v) DTE, 2% (w/v) SDS, 30% (v/v) glycerol for 12 min. Sulphydrilic groups were blocked in 2.5% (w/v) iodoacetamide, 50 mM Tris-HCl pH 6.8, 6 M urea, 2% (w/v) SDS, 30% (v/v) glycerol, bromophenol blue for 5 min.

The SDS–PAGE (Sodium dodecyl sulfate polyacrylamide gel electrophoresis) dimension was run in a vertical gradient acrylamide/PDA (9–16% T/2.6% C) slab gel. Sodium thiosulfate was used as an additive to reduce background in the silver staining. A constant current of 40 mA/gel was applied [23]. Gels were removed from glass plates, washed in deionized water for 5 minutes, and stained with ammoniacal silver as described by [22]. Preparative gels were stained with the Protea silver stain kit compatible with mass spectrometry analysis (Protea Bioscience, Morgantown, WV, USA).

Image analysis

The GS-700 Densitometer Gel Doc (Bio-Rad Laboratories, Hercules, CA, USA) was used as scanning device. Protein spots were detected using ImageMaster 2-D Platinum software, version 6.0 (GE Healthcare, Little Chalfont, UK). Spot borders were visually inspected and misidentification caused by confluent spots, artifacts and low signal to noise ratio, were manually corrected. Parameters like “saliency” (a measure of spot curvature) and “min area” (lowest area threshold under which spots are considered artefacts) were used to identify ‘true’ and ‘false’ protein spots. Manual contour drawing was then applied in all cases of sub-optimal spot auto-detection. This procedure was validated by assessing total spot numbers and spot volume ratios before and after background subtraction.

In order to optimize quantitative analysis of protein spots, the volume of each candidate spot was normalized using 4 surrounding landmark spots localized in areas closeby to the spot of interest, i.e. within corresponding background grey and with analogous signal staining exposure. Landmarks ratios were used for first-level normalization. Sums of landmark signals were then used for target spot normalized quantification. Using these criteria, spot detection and quantification were obtained, which minimized intensity variations among the gels, as assessed by Mann-Whitney test (p<0.05). Spots with more than 50% of data missing were not included in subsequent analytical steps.

Expression data of each identified spot were plotted into a frequency histogram to highlight main differences between analyzed sets, to visualize subtypes within the same histopathological group and to assess valued distribution shapes. The Shapiro-Wilk test was utilized to assess for normal (Gaussian) distributions.

Mass spectrometry analysis

After tryptic in gel-digestion, overnight at 37°C [24], the differentially expressed spots were excised from preparative gels and analyzed by mass spectrometry (MS) to identify amino acid sequences, using a Bruker Ultraflex III (Bruker, Bremen, Germany) operating in reflectron mode. This instrument was equipped with a Nd:YAG smartbeam laser to acquire positive-ion MALDI mass spectra over a mass range of m/z 800–4000. Spectral processing and peak list generation were implemented by Bruker flexAnalysis software (version 3.3, Bruker Daltonics) for MS and MS/MS spectra. For each protein spot, the most intense precursor ion signals in each MS spectrum were analyzed by MS/MS fragmentation in LIFT mode. α-cyano-4-hydroxycinnamic acid was used as matrix. Spot identifications were performed by querying the Mascot database. Trypsin cut, carbamidomethyl (C) as fixed, oxidation (M) as variable, a maximum of one missed cleavage allowed, were imposed as modifications in the search parameters. Peptide tolerance and MS/MS tolerance were set at 250 ppm and 0.5 Da respectively.

Univariate and multivariate statistical analysis

Univariate statistical analyses were performed with GraphPad Prism (GraphPad Software Inc., La Jolla, Ca) and XLStat (Addinsoft, Paris, France). Spearman's correlation analysis was performed using MetaboAnalyst 2.0 software (www.metaboanalyst.ca) [25][27]. Multivariate statistical analysis and data modeling were performed using MetaboAnalyst 2.0 (www.metaboanalyst.ca) [25][27] and SIMCA 13 (Umetrics, Umea, Sweden) [28] software packages.

Column-wise normalization was applied to provide Gaussian-like distributions [25], [29]. Analyses were then performed on autoscaled data (mean-centered and divided by the standard deviation of each variable) [30]. A diagnostic plot was utilized to represent normalization procedures for normal distribution assessments [29]. As examples, value intensities for e.g. APOA1, PRDX2, ALDOC, CRYAB_b, TTHY are ≥3- fold higher than others (e.g. NDUS1, QCR1, NFM, ACTB), thus inducing a skewed distribution. After autoscaling normalization, box plots have nearly same mean, standard deviation, and their distribution better matches a Gaussian curve (Figure S1B right) (Kernel density plot right-bottom) [31]. PCA was used as an unsupervised method in order to find the directions of maximum covariance among our protein spots without referring to class labels (tissue samples). This allowed to visualize differences among samples, to detect clustering and pick-up outliers.

PCA condenses datasets to obtain optimal dimensions that best capture signal covariance. However, it fails providing working hypotheses for some causal relations among data subsets [32]. Hence, we went on performing histopathology classification-guided PLS-DA. As many supervised classification algorithms tend to overfit the data [26], [33], PLS-DA model validation was performed as previously described [34]. Briefly, to define the optimal number of PCs (principal components), “7-fold cross-validation” (CV) was applied [35], [36]. Using CV, the predictive power of the model was verified. Two parameters were calculated for evaluating the models: R2 (goodness of fit) and Q2 (goodness of prediction). A model with Q2>0.5 was considered good, Q2>0.9 excellent [37], [38]. As cross-validation only assesses the predictive power without a statistical validation, the performance of PLS-DA models was also validated by a permutation test (200 times).

To help interpreting results from PLS-DA, we considered the variable importance in the projection scores (VIP score) and regression coefficients (CoeffCS). This allowed to evaluate protein influence (including prediction performance) on the model and identify the best descriptors of the differences among the three groups. The VIP score is a weighted sum of squares of the PLS loading weights taking into account the amount of explained Y-variation in each dimension [25], [26], [39]. Since the average of squared VIP scores equals 1, the “greater than 1” rule is generally used as a criterion to identify the most significant variables [37], [39]. PLS-DA CoeffCS express the relation between the Y variables (classes) and all the terms in the model and are used for interpreting the influence of the X variables (proteins) on Y. VIP and CoeffCS values are cumulatively calculated from all extracted PLS components. The coefficients express how strongly the Y groups are correlated to the systematic part of each X variable considering all three components.

Gene ontology, networks and functional analyses

Gene Ontology (GO) analysis was performed using PANTHER 7.2 software (www.pantherdb.org/). The signaling hubs and connectivity networks were obtained using Ingenuity Pathway analysis (IPA, Ingenuity Systems, www.ingenuity.com) and STRING 9.1 (string-db.org) package.

Western blotting

Expression levels of randomly selected proteins were analyzed by Western blotting in order to validate the dataset identified by proteomic analysis. Blots were incubated with the following primary antibodies (Santa Cruz, Santa Cruz, CA): LDH-B (Q-21) sc-133731 rabbit polyclonal; SOD-1 (V-17) sc-34015 goat polyclonal; APOA-I (FL-267) sc-30089 rabbit polyclonal; Aldolase C Antibody (N-14) sc-12065 goat polyclonal and PRXII (9A1) sc-59660 mouse monoclonal. Signal intensities of the bands were quantified with Image JA 1.46b, using a Kodak grey-scale standards power curve (www.kodak.com) as reference. Band intensity values were normalized versus red Ponceau signals of transferred proteins on Western blot filters. Normalized densitometry values between proteomic gel spots and Western blot bands were correlated with the Spearman's rank correlation analysis, to obtain rho coefficients and corresponding p values.

The expression profiles of hub proteins (Huntingtin, HNF4α, c-Myc, 14-3-3ζ) were also determined on glioma tissue lysates.

Immunohistochemistry staining

Five µm sections from FFPE samples were stained overnight using antibodies to Huntingtin (Millipore MAB2166; mouse monoclonal, clone 1HU-4C8, 1∶150), HNF4α (AbCam ab41898; mouse monoclonal, clone K9218, 1∶70), 14-3-3ζ (AbCam ab51129; rabbit polyclonal, 1∶50) and c-Myc (AbCam ab32072; rabbit monoclonal, clone Y69, 1∶100). Antigen retrieval was performed using hot citrate buffer pH 6.0 (Huntingtin, HNF4α, and 14-3-3ζ) or 1 mmol/L EDTA pH 8.0 (c-Myc). Antigen-antibody reactions were visualized using a polymer-based detection system (EnVision Kit, Dako), using diaminobenzidine as chromogen.

Immunohistochemistry data bank meta-analysis

Publicly available databases [40], [41], containing high-resolution IHC (Immunohistochemistry) images extending proteome-wide were analyzed for patterns of expression of GBM-driving engines. The Human Protein Atlas (v. 12, www.proteinatlas.org) provides spatial distribution and expression data from 16621 proteins/21984 antibodies in different normal human tissues and different cancer types. The expression profiles of hub proteins were generated for antibody staining parameters, intensity, and fraction of positive cells in control versus glioma arrays. EGFR and p53 IHC staining were used as internal benchmark for performance assessment and quantification standards.

Results

Proteomic analysis

Twenty-five high-grade GBMs, 10 low-grade gliomas and 5 tissue samples from normal brain cortex were analyzed by 2D PAGE (Figure 1A). Clinical data from brain tumor patients are summarized in Table 1 [21].

thumbnail
Figure 1. Gel analysis, spot detection and edge improvement.

(A) Representative 2-DE gels from normal brain, low-grade and high-grade gliomas. Proteins from normal and tumor brain tissues were processed as described in materials and methods. Protein spots were visualized by ammoniacal silver staining. (B) 3-D representation of protein volume spots and comparison between automatic and manual procedures for definition of protein spots and edges. Manual strategy is time-consuming but allows reducing the number of misidentifications and improving the quantification of the same protein markers between different gels. (C) Normalization of density values of differentially expressed spots (green) by using surrounding landmarks (red). (D) High-grade histotype tumors (green) showing a bimodal distribution of triosephosphate isomerase (TPIS1) expression values. In low-grade tumors (blue) and control samples (red) the values appear to follow a typical Gaussian (normal) distribution.

https://doi.org/10.1371/journal.pone.0103030.g001

To improve quantification accuracy and allow robust statistical analysis of 2D gel data [42], we optimized image processing, spot detection and signal quantification procedures by applying operator-guided background subtraction and spot contour optimization (Figure 1B and S1A). These procedures allowed to obtain 49±15% increase of the detected signal (spot volume upon spot contour optimization) (Table S1), and better spot detection (+16±2%), as compared with the basic/automatic procedure. Most spot variation was not detectable across all samples; therefore, protein spots were marked as “differential” if expression changes could be observed in at least one tumor sample versus all control samples. The statistical distributions were found to be not-normal and often bimodal for most protein spots (Figure 1D). Forty-eight protein markers were identified by MS (Table 2, Figure S2, Supplemental Material S1), and their expression levels were quantified upon the image optimization procedures described above (Figure S3A and Table S2A,B). For enhancing quantification robustness, the density value of each candidate spot was normalized versus 4 surrounding landmarks (Figure 1C, Table S2A) as described in Materials and Methods section. This allowed to obtain robust quantification, through compensating for local staining dishomogeneities.

GO analysis was performed, revealing the largest protein classes to have metabolic (26 proteins), structural (9) molecule binding (6) and antioxidant (4) activities (Figure 2A, top). Dual/multiple activities accounted for redundant/over-represented functions. Most detected proteins were hydrolases (10), oxidoreductases (10), components/interactors of the cytoskeleton (8), transfer/carrier proteins (5) and transferases (4). (Figure 2A, bottom).

thumbnail
Figure 2. Classification of identified spots and correlation analysis.

(A) GO pie charts show PANTHER classifications made according to the associated Molecular function (top) and Protein class (bottom). (B) Graphical representation of Spearman's correlation matrix. Heatmap shows Spearman's correlation between differentially expressed protein spots. Each column and row defines an individual variable. Positive correlation values are in red, and negative correlation values are in blue. Hierarchical clustering was applied to both dimensions. (C) Positive (rho ≥0.5) and negative (rho ≤−0.5) correlations are listed in red and blue, respectively.

https://doi.org/10.1371/journal.pone.0103030.g002

Spearman's correlation analysis

Negative or positive correlations were globally revealed by Spearman's correlation analysis (Figure 2B). Highest positive correlations (Figure 2C, red) were found to occur between HBA and HBD; CN37 and TBA1B; CN37 and LDHB; ALDOC and ESTD; CN37 and DEST; CN37 and RAN; ALDOC and NFM; TBA1B and DEST. Highest negative correlations (Figure 2C, blue) were observed between UCHL1 and HBD, and between UCHL1 and HCD2 (Table S3).

PCA analysis

The proteomic matrix was processed by scaling protein expression values in order to reduce potential systematic bias and make the variables comparable in magnitude to each other [25], [27], as indicated. The data scaling results and normalization procedures, are summarized graphically in the Figure S1B. The horizontal box plots represent the distributions of individual variables, the bottom curves show the global data distribution based on kernel density estimation (Figure S1B).

Then, we went on to utilize PCA as an unsupervised multivariate method for analyzing the dataset and identifying the best discriminators among sample classes. PCA score plots were generated (Figure 3A), where each axis represented a PC identifying linear combinations of the most tightly interconnected proteins/signaling networks [32], [43]. Samples with similar protein expression profiles/PC scores clustered together with striking fitness. PC1 (score vector t1) was found to discriminate controls from tumors; PC2 (score vector t2) separated low-grade from high-grade gliomas (Figure 3A, left). Seven major discriminators between control and tumor samples were found: expression levels of APOA1, PRDX3_a and CLIC1 were higher in tumors than in normal brain cortex, whereas significantly lower levels of NFM, CN37, NDUS1 and MDHC were found in tumors as compared with normal tissues (component PC1, Figure 3A, right). Thirteen major discriminators between low- and high-grade tumor samples were found: expression levels of HCD2, HBA and HBD were strongly up-regulated in high-grade gliomas, whereas CRYAB_b, IPYR, TPIS, PEA15, PSD13, GFAP, PHP14, 6PGL, KCRB, IDH3A had higher expression in low-grade than high-grade tumors (component PC2, Figure 3A, right). PCA analysis also revealed that UCHL1 have a high discriminating power of this marker when control/low-grade samples were compared with high-grade tumors. Global sets of protein marker with higher (Figure S3B,C) and lower discriminating power (Figure S3D) were identified.

thumbnail
Figure 3. PCA and PLS-DA models.

(A, left) PCA score plot showing separation between control samples, low-grade and high-grade tumors. (A, right) PCA loading plot showing the proteins (variables) responsible for discrimination between the groups. (B) PLS-DA score plot showing separation between control samples, low-grade and high-grade tumors. (C) Overlapping of PCA correlation groups (as shown in Figure 3A, right panel) to the PLS-DA weight plot in order to enhance the discriminating power of identified markers. Colored ovals and solid circles represent PLS-DA and PCA protein clusters, respectively. Color code: red, controls; blue, low-grade tumors; green, high-grade tumors; magenta, low and high-grade tumors; yellow, low-grade tumors and controls.

https://doi.org/10.1371/journal.pone.0103030.g003

PLS-DA analysis

To verify the strength of the unsupervised PCA analysis, and to further build on it, we analyzed the dataset on the basis of known classes (controls vs high-grade tumors vs low-grade tumors) using a supervised PLS-DA method [37], [38], [44], [45]. This model was found to have strong goodness of fit (cumulative R2Y = 0.890) and prediction power (cumulative Q2 = 0.813) (Figures 4 A,B). The separation between normal brain tissues, low-grade and high-grade gliomas yielded a staggering clear discrimination (Figure 3B). Most significant separations were explained by a three-component model, where principal components PC1, PC2 and PC3 represented 15.1%, 13.7% and 10% of the total variance in the protein spot-matrix, respectively (Figure 4). Permutation tests were carried out in order to validate the PLS-DA model [37]: as shown in Figure S4A, the original model was found to have higher R2 and Q2 values than the permuted models, and negative Q2 values were obtained for all three permuted groups tested.

thumbnail
Figure 4. PLS-DA cross-validation, performance and protein VIP scores.

(A–B) Bar plot showing the performance measures (R2Ycum and Q2cum) using different number of components. The selected performance measure Q2 shows the three-component model performs as the best one. R2X: portion of the variation of X explained by specified PC; R2X(cum) Cumulative explained portion of X set variation; Eigenvalue: number of variables (K) times R2X; R2Y: portion of the Y set variation modeled by the PC; R2Y(cum): cumulative modeled variation of Y set; Q2: overall cross-validated R2 for the specific PC; Limit: threshold cross-validation for the specific PC; Q2(cum): cumulative Q2 up to the specified component, is a model predictive power according to cross validation. Unlike R2X(cum), Q2(cum) is not additive. (C) 3-D score plot. The samples are represented in the 3-D score plot, the first three components (PC1, PC2, PC3) are reported accounting 15.1%, 13.7% and 10% of total variation respectively. (D) Proteins able to discriminate between controls, low-grade tumors and high-grade tumors, ordered by VIP score. VIP scores ≥1 are significant (above the red line) and indicate important X variables (proteins) that predict Y responses (tissue samples).

https://doi.org/10.1371/journal.pone.0103030.g004

A PLS-DA loading plot was generated in order to find major discriminants between the groups analyzed (Figures 3C). Eleven major discriminators between control and tumor samples were identified: APOA1, CLIC1 and PRDX3_a were over-expressed, whereas NFM, CN37, NDUS1, MDHC, ALDOC, STMN1, PEBP1 and DDAH1 were down-regulated in tumors as compared with control samples. Twelve major discriminators between low- and high-grade gliomas were identified: HCD2, HBA and HBD were up-regulated in high-grade gliomas, whereas expression levels of CRYAB_b, IPYR, TPIS, PEA15, PSD13, GFAP, IDH3A, 6PGL and PHP14 were found to be higher in low-grade than in high-grade tumors. The regression coefficients calculated for PLS-DA outcomes confirmed the power of the identified clusters in distinguishing among sample classes (Figure S4B). Next, we calculated the VIP score for each protein in our dataset. Out of 48 variables analyzed as potential predictors, the following 18 descriptors were found to significantly contribute to the classification model (VIP score ≥1): HBD, UCHL1, HBA, STMN1, HCD2, ALDOC, NFM, IPYR, NDUS1, MDHC, DDAH1, PSD13, APOA1, 6PGL, PEBP1, TTHY, ACTY, IDH3A (Figure 4D).

In order to improve the discriminating power of the identified protein markers, we went on to intersect PCA and PLS-DA loading plots and find the shared proteins/best discriminators for each condition (Figure 3C). HBA, HBD and HCD2 positively correlated with high-grade gliomas; GFAP, PHP14, 6PGL, PSD13, PEA15, TPIS, CRYAB_b, IPYR and IDH3A correlated with low-grade tumors; on the other hand, NFM, CN37, NDUS1 and MDHC negatively correlated with tumor samples. APOA1, PRDX3_a and CLIC1 were the best discriminators between tumors and negative controls.

Validation of proteomic profiles

The proteomic landscape of human gliomas was then validated by protein immunoblotting, quantifying the expression levels of randomly-selected spots. Protein markers (APO-A1, SOD1, PRDXII, LDHB, ALDOC) were randomly selected (∼10%) among the 48 differentially expressed proteins and analyzed. Western blot chemiluminescence images were acquired at sub-saturation levels and quantified with ImageJ. Silver staining density quantified as above and Western blot signals were then subjected to Spearman's correlation analysis. Paired signal analysis supported the accuracy of proteomic profiles: APOA1 (ρ = 0.550, pvalue = 0.015), SOD1 (ρ = 0.517, pvalue = 0.025), LDHB (ρ = 0.517, pvalue = 0.044), PRDXII upper band (ρ = −0.086, pvalue = 0.387), PRDXII lower band (ρ = 0.030; pvalue = 0.458), ALDOC (ρ = 0.771, pvalue = 0.051). Table S4 presents in extenso Western blotting data, densitometry, normalization procedure details, scatter plots and elliptic confidence intervals.

Pathway analysis

The discriminating performance of the protein clusters identified by unsupervised analysis, and the tight correspondence between PCA and PLS-DA clusters suggested deep, intersected biological relevance. Hence, we went on to explore the existence of a “GBM control module” connecting the discriminating protein clusters through physical and functional interactions. In order to reveal these cross-talks we carried out bioinformatic network detection followed by data meta-analysis. The top-score network (Fisher's exact test: p = 1x10−104) was found to contain 46 out of the 48 proteins identified by proteomics (Figure 5 and Table S5A). Most members were found to play key roles in neurological diseases (28), genetic disorders (34), skeletal and muscle diseases (22), and cancer (25). The five top score pathways of this network included proteins involved in mitochondrial function (HCD2, PRDX3, NDUS1, PRDX5 and QCR1), pentose phosphate pathway (6PGL, ALDOC, TKT); glycolysis/gluconeogenesis (HCD2, TPIS, LDHB and ALDOC), inositol metabolism (TPIS and ALDOC) and oxidative phosphorylation (NDUS1, ATP5H, QCR1 and IPYR).

thumbnail
Figure 5. Pathway analysis.

Graphical representation of the protein network retrieved using the Ingenuity Pathway Analysis Tool. Proteins are represented as nodes, the biological relationships between the nodes are represented as lines. Proteins identified by PCA and PLS-DA analysis and indicated as discriminant among controls, low-grade and high-grade tumors are in highlighted blue. Most proteins differentially expressed in gliomas are closely connected in the network through four major hubs: Huntingtin, HNF4α, 14-3-3ζ (YWHAZ) and c-Myc. Four external edges: differentially expressed proteins identified by MS (red) and direct interactors with the four major hubs (blue bridges and strengthened contours).

https://doi.org/10.1371/journal.pone.0103030.g005

Strikingly, most proteins of the identified network were then found to converge on four major hubs: Hungtintin (HTT, 16 edges), Hepatocyte nuclear factor 4α (HNF4A, 10 edges), 14-3-3ζ (YWHAZ, 9 edges) and c-Myc (MYC, 9 edges) (Figure 5). Remarkably, major discriminators identified by PCA and PLS-DA analyses were found to interact with Huntingtin (10), HNF4α (5), c-Myc (4) and 14-3-3ζ (3) (Table S6).

We then went on to assess the relevance of the GBM control module versus known molecular pathways involved in GBM biogenesis, via IPA and STRING platforms algorithms (Tables S5B). To our surprise, the four hubs we had identified were found to converge on both major players of glioma development and progression: epidermal growth factor receptor (EGFR) and p53. Bridging proteins between network hubs and EGFR/p53 were as follows: UCHL1, TPI1 and SH3BGRL to EGFR; ACTB, CRYAB, STMN1, NME1, Tubulin, GFAP, UBE2N, PPA1 and UCHL1 to p53 (Table S5B).

Contrary to a wide-held belief, proteins and mRNA levels correlate poorly in most cellular systems [46][48], differential protein/mRNA stability playing a major role in this discordant scenario [49], [50]. Nevertheless, identification of transcription factor-driven differential gene expression landscapes provides insight into tumor-driving gene networks [51]. Hence, we went on to identify upstream transcription factors potentially involved in a coordinated regulation of proteins taking part to the GBM control module. Using stringent criteria for the analysis (P values <0.005, interactions ≥5), we discovered that 9 transcription factors (HTT, MYC, HNF4A, TP53, ESRRA, NFE2L2, PPARGC1A, MYCN, ESR1) interacted with 33 out of 48 differentially expressed proteins. Importantly, these transcription factors were also found to interact with/regulate expression of 18 of the best discriminators identified by PCA and PLS-DA (Table S5C). Of major relevance, all transcription factors identified as central hubs (Huntingtin, c-Myc, HNF4α) of the GBM control module, together with p53, stood-up as major drivers of the expression of the vast majority of the components of the module.

Validation of hub expression in tumors

A prediction of our model was that the four hubs of the GBM control module should be broadly expressed. Hence, we assessed their expression in human glioma samples by protein immunoblotting. 14-3-3ζ, HNF4α and Huntingtin were widely expressed in glioma samples. Expression of HNF4α was higher in tumors samples than controls; Huntingtin and c-Myc were found to be overexpressed in high-grade gliomas (Figure 6). Moreover, we analyzed the expression of the four hubs in tissue samples by IHC (Figure 7 and S5). In high-grade gliomas 14-3-3ζ and HNF4α were strongly expressed in nuclear and cytoplasmatic compartments (Figure 7A, Figure 7C). On the other hand, specific nuclear accumulation was observed in low-grade samples (Figure 7B, Figure 7D). c-Myc was specifically expressed in the nucleus of the glioma samples, with a trend toward an up-regulation from low- to high-grade tumors(Figure 7G, Figure 7H).

thumbnail
Figure 6. Network hubs expression in glioma samples - Western blot analysis.

14-3-3ζ, HNF4α Huntingtin and c-Myc protein expression levels in tumor and control samples, as determined by Western blotting. GB: glioblastoma mutiforme; OL: oligodendroglioma; PA: pilocytic astrocytoma; FA: fibrillary astrocytoma. Control sample (CTR) was tissue from normal cortex. HG: high-grade tumors. LG: low-grade tumors.

https://doi.org/10.1371/journal.pone.0103030.g006

thumbnail
Figure 7. Network hubs expression in glioma samples - IHC analysis.

Expression of 14-3-3ζ (A, B), HNF4α (C, D), Huntingtin (E, F) and c-Myc (G, H), as determined by IHC staining of glioblastoma samples (left column) and low-grade tumors (right column). Representative samples are shown. Scale bars  = 20 µm. Nuclei were counterstained with hematoxylin (in blue).

https://doi.org/10.1371/journal.pone.0103030.g007

Four hubs expression meta-analysis

The findings above suggested broad expression of the four hubs of the GBM control module. We verified this prediction by performing a proteome-wide profiling of IHC expression patterns (Human Protein Atlas; www.proteinatlas.org) for Huntingtin, HNF4α, 14-3-3ζ and c-Myc. HNF4α level in normal glial cell showed low levels of staining, with moderate intensity in <25% of the cells. In analyzed gliomas 5/10 had a corresponding expression profile as compared to controls; 2/10 presented an increase of expression (medium staining, moderate intensity and percent reactive cells of 75–25%), 3/10 presented lower expression (<25% of cells) compared with control samples. Hungtintin level in normal glial cell showed low expression (low staining, moderate intensity and percentage <25%) in IHC array stained with mouse mAb. In glioma tissue arrays 0/12 have the same expression profiles compared to controls. Strikingly, 12/12 presented a global increase of expression or a substantial increase of positive cells. A second IHC array set (12 samples) stained with rabbit polyclonal antibody was analyzed confirming this evidence. Consistent with our findings, c-Myc expression in normal glial cells was not detectable by IHC. Rabbit polyclonal antibody targeting the C-terminal portion of c-Myc led to positive staining on 4/11 tumors samples. The mouse mAb gave positive staining in 11/12 astrocytoma samples. 14-3-3ζ presented strong staining levels in normal glial cell (high staining, strong intensity, percentages between 75%–25%). Three out of 10 array samples presented the same staining patterns as controls. The remaining 7/10 samples showed a prevalence of positive cells of >75%. p53 and EGFR expression patterns were analyzed as internal benchmarks of the robustness of analysis and were shown to possess expected expression profiles and prevalence of expression findings (Table S7).

Discussion

Proteomic profiling of human GBM allowed to discover differentially expressed protein clusters, that were shown to craft a tightly interconnected control network. This was recapitulated into a four-hub control module, as centered on Huntingtin, HNF4α, c-Myc and 14-3-3ζ. This was able to stringently discriminate between high-grade GBMs, low-grade tumors and normal tissues. The proteomic clusters included tumor upregulated (PRDX3, APOA1, CLIC1) and downregulated (NFM, NDUS1, MDHC, ALDOC, STMN1, PEBP1, DDAH1, CN37) proteins. Major discriminator between high-grade and low-grade tumors included CRYAB, IPYR, TPIS, PEA15, PSD13, GFAP, IDH3A, 6PGL, PHP14, KCRB as overexpressed in low-grade gliomas; HCD2, HBA, HBD as overexpressed in high-grade GBM. UCHL1 expression showed a positive correlation with normal brain tissue and low-grade tumor, and a negative correlation with high-grade tumors.

Huntingtin, whose mutations are responsible for the neurodegenerative disorders of Huntington's disease, is found in neurites and at synapses, has anti-apoptotic functions and is neuroprotective in brain cells exposed to apoptotic stimuli, such as serum deprivation, mitochondrial toxins or death-inducing genes [52]. Notably, pathogenic Huntingtin affects the expression, redox state, disulfide bonding of antioxidant proteins identified here, among them SODC, and PRDX2, together with PRDXI [53], thus supporting a shared functional link. Taken together, our findings provide first evidence of function of Huntingtin in brain tumors, thus paving novel avenues of investigation on GBM pathophysiology.

HNF4α is a modulator of cell proliferation [54][56] through the cell cycle inhibitor p21 [57] and the transmembrane glycoprotein Trop-2 [51]. A tight interplay/feedback loop occurs between HNF4α and c-Myc. Both HNF4α and c-Myc proteins compete for control of the P21/CDKN1A gene transcription [57], and deletion of HNF4A in hepatocellular carcinoma cells results in significant up-regulation of c-Myc and enhanced cell proliferation rates [58]. Essentially no evidence for expression and function of HNF4α in brain tumors was available before this study, again opening novel avenues for investigation on GBM pathophysiology.

Deregulation of MYC is a frequent driver of cancer [59]. c-Myc has been reported to bind a large number of genes [60] and regulates cell proliferation by affecting cell-cycle checkpoint genes, CDK inhibitors and cyclins [61]. c-Myc also plays a major role in regulating metabolic genes required for energy production [62], [63] and ribosomal biogenesis. mTORC2 controls glycolytic metabolism by regulating c-Myc cellular levels and ultimately determines overall survival of GBM patients [64]. This evidence is consistent with our GO analysis showing that major targets of our analysis, such as ALDOC, TPIS, IPYR, 6PGL, HCD2, IDH3A, NDUS1, MDHC are involved in metabolism (e.g. glycolysis, inositol metabolism and oxidative phosphorylation) and cancer-related metabolic reprogramming, including the Warburg effect [65][69].

Our findings support a major involvement of 14-3-3ζ in the progression of GBM [70], in agreement with previous studies showing that 14-3-3ζ expression levels were a prognostic factor in GBM [71]. 14-3-3ζ is involved in oral squamous cell [72], stomach [73], breast [74] and papillomavirus-induced carcinomas [75]. Major targets identified in our analysis, such as GFAP, CRYAB and TPIS, are major interactors of 14-3-3ζ and are powerful discriminators of low-grade astrocytoma.

Glioma development is frequently associated with mutations of the isocitrate dehydrogenase IDH1 and IDH2 genes [76], [77], whereas mutations of IDH3 have never been observed in GBM [78]. Our analysis discovered IDH3A quantitative variations in low-grade samples versus high-grade tumors. This provided support for a novel model of interference of IDH proteins in GBM progression, through differential expression of a wild-type protein. Of interest, IDH3A was found to be a specific target for p53-dependent phosphorylation [79], further supporting the functional relevance of the GBM control module.

Remarkably, three of the four hubs of the GBM control module (Huntingtin, HNF4α, c-Myc) are transcription factors. Transcription factor network analysis then highlighted all three of them as major regulators of the expression of most proteins of the GBM control module, supporting a joint driving role in GBM development. A functional role of the newly discovered four-hub control module in GBM appearance and progression further required vast, coordinate expression in tumor cases.

Strikingly, analysis of the transcription factors steering the GBM control module led to the discovery that these tightly interrelate with p53. Notably, p53 can regulate Huntingtin's expression at the transcriptional level [80], thus suggesting a cooperation of these signaling pathways not only in neurological diseases, but also in development and progression of brain tumors. p53 plays a critical role as modulator of the HNF4α/c-Myc feedback loop, since it binds c-Myc [81] as well as HNF4α [82], and inhibits the activity of HNF4α via recruitment of histone deacetylase [82]. Additionally, ATM-dependent activation of p53 involves dephosphorylation and association with 14-3-3 [83].

GBM type II are linked to mutations of TP53, whereas GBM type I are thought to be driven by EGFR amplification/disregulation [5], [7]. Notably, overexpression of Huntingtin interacting protein 1 (HIP1) has been shown to correlate with increased EGFR levels [84]. HIP1 physically associates with EGFR and maintains its levels in brain tumors [84]. It was recently reported that EGFR induces expression of the oncogenic miRNA miR-7 through a Ras/ERK/Myc pathway, and that c-Myc binds to the miR-7 promoter, enhancing its activity [85]. The 14-3-3ζ protein was reported to directly bind EGFR upon stimulation with EGF [86]. Recently, downregulation of PEPB1/RKIP (Raf kinase inhibitor protein) was shown to be associated with poor outcome and malignant progression [87]. PEPB1 inhibits RTKs signaling blocking Raf/MEK/ERK cascade. We found PEPB1/RKIP downregulated in low-grade and high-grade tumors, extending previous indications [20]. Taken together, our findings indicate deep intertwining of the GBM control module also with EGFR signaling pathway.

In summary, using multitiered proteomic profiling, we discovered previously unidentified hubs (centered on Huntingtin, HNF4α, c-Myc and 14-3-3ζ; these then attract p53 and EGFR) as major signaling drivers in the pathogenesis of brain tumors. Our findings thus support the unexpected existence of a unique GBM control module, helping providing a much needed unifying model for GBM appearance and progression. Future studies should be undertaken to validate a diagnostic/prognostic role of the GBM control module. This may also provide better tools for classification and clinical evaluation of GBM, for more effective procedures for tumor diagnosis, prognosis and patients cure.

Supporting Information

Figure S1.

Signal versus noise in spot detection. (A) The gel images were subjected to automatic spot detection setting the same parameters: the number of detected spots was increased in the gray adjusted image (right) as compared with the original one (left). (B) Data normalization view. Box plots and kernel density plots show the distribution of protein concentration before (left) and after (right) autoscaling (mean-centered and divided by the standard deviation of each variable) as described in the text.

https://doi.org/10.1371/journal.pone.0103030.s001

(TIF)

Figure S2.

Protein Identification by mass spectrometry. (A) Red circles indicate spots excised on preparative gels and subjected to in-gel tryptic digestion, followed by MS and MS/MS spectrometry analysis for protein identification. Gels and samples were processed as described in materials and methods section. (B) Examples of mass spectra from identified proteins. Numbers on X axis represent precise m/z values of detected peptide ion signals. The peak masses were used to identify the proteins. For each protein spot the strongest peaks were analyzed by MS/MS fragmentation in LIFT mode. α-cyano-4-hydroxycinnamic acid was used as matrix. (top) P02511, Alpha-crystallin B chain (CRYAB). (middle) P50213 Isocitrate dehydrogenase NAD subunit alpha, mitochondrial (IDH3A). (bottom) P10768, S-formylglutathione hydrolase (ESTD).

https://doi.org/10.1371/journal.pone.0103030.s002

(TIF)

Figure S3.

Differential spot expression analysis. (A) Differentially expressed proteins spots as quantified by image analysis (Materials and Methods) (B–D) Examples of and proteins with multimodal distribution in glioblastomas, low-grade astrocytomas and control samples; protein with higher (APOA1, NFM) or lower discriminating power (TTHY) are shown.

https://doi.org/10.1371/journal.pone.0103030.s003

(TIF)

Figures S4.

PLS-DA model permutation test plots and coefficient scores of proteins from the PLS-DA analysis. Permutation tests for: High-grade tumors (left), low-grade tumors (middle) and controls (right). Permutation tests were performed by comparing goodness of fit and prediction (R2 and Q2 values) of the original model with the goodness of fit and prediction of several models based on data in which the order of the Y observations were randomly permuted. The two intercepts can be considered as measures of degrees of overfit and overprediction. The correlation coefficients of original and permuted data are reported on the x axis; 200 random permutations were carried out. The values of R2 and Q2 are reported on the y axis. The two circles on the in the upper right (ρ = 1) correspond to the values of R2 (green circles) and Q2 (blue circles) of the original data. The other circles represent permutation results. The low values of intercepts show that the model has a statistical significance (not over-fitting). (B) Coefficient scores were utilized to provide an estimate of the protein changes in the various groups. Larger coefficient scores (positive or negative) indicate stronger correlations with proteomic group profile classification. The highest positive (black boxes) or negative (red boxes) discriminating coefficient scores of high-grade tumors (left) were exemplified by translation to low-grade tumors (middle) and controls (right).

https://doi.org/10.1371/journal.pone.0103030.s004

(TIF)

Figure S5.

IHC staining in positive control samples. Expression of the 4-hub proteins in positive control tissue sections. (A) Expression of Huntingtin in small bowel. (B) Expression of HNF4α in small bowel. (C) Expression of 14-3-3ζ in breast cancer. (D) Expression of c-Myc in colon cancer. Scale bars  = 20 µm. Nuclei were counterstained with hematoxylin (in blue).

https://doi.org/10.1371/journal.pone.0103030.s005

(TIF)

Table S1.

Inter-gel spot volume quantification variance analysis. Three replica gels of a reference protein lysate from high-grade brain tumor were utilized for inter-gel variance analysis. Ten randomly chosen spot volumes were quantified using Image Master Platinum 6.0; basic/automated procedure was compared with operator-guided contour drawing. The bar graph shows the gain-of-signal (expressed as percent of variation) detected using operator-guided contour drawing versus basic/automated analysis. P = 0.037 (comparison among means in manual vs basic/automated procedures). SD, standard deviation. a: values from basic-automated quantification procedure. b: values from operator-guided quantification procedure. c: values from operator-guided quantification procedure have been normalized on values from basic-automated procedure, and expressed as percent of variation.

https://doi.org/10.1371/journal.pone.0103030.s006

(XLSX)

Table S2.

Table S2A: Differentially expressed proteins in tumor samples versus normal brain. 2D-gels (5 Control samples, 10 Low-Grade and 25 High-Grade tumors) were processed and quantified as described. Density values from differentially-expressed protein spots were determined. Density normalization on local landmarking was performed as described in the main text. Table S2B: Protein level distributions in normal brain cortex and tumor samples. Box and scatter plots of protein markers defined by proteomics analysis. The graphs show normalized density values. The boxes encompass values from the first quartile (bottom) to the third quartile (top) for the three category (CTR = control; LG = Low-Grade; HG = Hi-Grade). Red horizontal line, median value. Red cross, average value. Each black dot represents an individual sample.

https://doi.org/10.1371/journal.pone.0103030.s007

(XLSX)

Table S3.

Spearman correlation matrix. Spearman's correlation matrix of all marker proteins identified by MS analysis. Numeric values of Spearman's correlation coefficients (ρ) between variables are reported. Each column and row show individual variables. Global correlation analyses are presented in Figure 3.

https://doi.org/10.1371/journal.pone.0103030.s008

(XLSX)

Table S4.

Validation of proteomic target proteins by immunoblotting analysis. Immunoblot analysis (second column) versus silver normalized density values (first column). Five proteins (∼10%) were randomly selected among the 48 differentially expressed proteins and analyzed in tumor samples. Density values from blots were quantified as described in Materials and Methods, as normalized on red Ponceau signal. Silver staining density and Western blot signals were subjected to Spearman's correlation analysis; correlation coefficients (rho) and p-values are reported. Scatter plots for the two variables with confidence ellipses were generated. Representative blots from APOA1, SOD1, LDHB, PRDXII (upper and lower bands) ALDOC are shown.

https://doi.org/10.1371/journal.pone.0103030.s009

(XLSX)

Table S5.

Table S5A: Ingenuity Pathway Analysis. The significance values for canonical pathways and other biological functions were calculated using the right-tailed Fisher's exact test by comparing the number of user-specified proteins that participate in a given function or pathway, relative to the total number of occurrences of these proteins in all pathway or functional annotations stored in the Ingenuity pathway knowledge base (IPKB). a: The degree of interaction between differentially expressed markers was compared with that expected by chance. A p-value  = 1×10−104 was computed by a hypergeometric test. Table S5B: Supervised pathway analysis. Interaction of EGFR (A) and p53 (B) with network proteins, as determined by IPA analysis. (C) Pathway analysis, as performed by STRING 9.1, of the four major hubs (HTT, HNF4A, Myc, YWHAZ) cross-interacting with p53 and EGFR. Table S5C: Transcription factor pathway analysis. Transcription Factor Analysis, as performed by IPA Upstream Regulator Analysis Tool. Using stringent cut-offs for interaction significance (p value <0.005); a threshold value for interaction with ≥5 target proteins was applied. Nine transcription factors (HTT, MYC, HNF4A, TP53, ESRRA, NFE2L2, PPARGC1A, MYCN, ESR1) were shown to modulate 33 out of 48 differentially expressed proteins. Color codes correspond to those of discriminating proteins by PCA and PLS-DA analysis (Figure 4). Proteins that positively correlate with controls are in red; with low-grade tumors are in blue, with both high-grade and low-grade are in magenta. Correlation of UCHL1 with low-grade/control group is in yellow.

https://doi.org/10.1371/journal.pone.0103030.s010

(XLSX)

Table S6.

Table S6A: Discriminator proteins from PCA/PLS-DA clusters mapping on network hubs - HTT. Table S6B: Discriminator proteins from PCA/PLS-DA clusters mapping on network hubs - HNF4A. Table S6C: Discriminator proteins from PCA/PLS-DA clusters mapping on network hubs - Myc. Table S6D: Discriminator proteins from PCA/PLS-DA clusters mapping on network hubs - YWHAZ (14-3-3 zeta protein). Protein members from the four discriminator clusters, as defined by PCA and PLS-DA analysis, were mapped on the IPA network (highlighted in blue).

https://doi.org/10.1371/journal.pone.0103030.s011

(XLSX)

Table S7.

Table S7A: Immunohistochemistry proteome profiling meta-analysis - Huntingtin. Representative examples of IHC-stained sections of glioma samples. Specific staining for the huntingtin protein is identifiable as brown spots. Nuclei are counterstained with hematoxylin (in blue). Table S7B: Immunohistochemistry proteome profiling meta-analysi - HNF4A. Representative examples of IHC-stained sections of glioma samples. Specific staining for the HNF4α protein is identifiable as brown spots. Nuclei are counterstained with hematoxylin (in blue). Table S7C: Immunohistochemistry proteome profiling meta-analysis - c-Myc. Representative examples of IHC-stained sections of glioma samples. Specific staining for the c-Myc protein is identifiable as brown spots. Nuclei are counterstained with hematoxylin (in blue). Table S7D: Immunohistochemistry proteome profiling meta-analysis - YWHAZ (14-3-3ζ protein). Representative examples of IHC-stained sections of glioma samples. Specific staining for YWHAZ (14-3-3ζ protein) is identifiable as brown spots. Nuclei are counterstained with hematoxylin (in blue). Table S7E: Immunohistochemistry proteome profiling meta-analysis - EGFR. Representative examples of IHC-stained sections of glioma samples. Specific staining for the EGFR protein is identifiable as brown spots. Nuclei are counterstained with hematoxylin (in blue). Table S7F: Immunohistochemistry proteome profiling meta-analysis - p53. Representative examples of IHC-stained sections of glioma samples. Specific staining for the p53 protein is identifiable as brown spots. Nuclei are counterstained with hematoxylin (in blue).

https://doi.org/10.1371/journal.pone.0103030.s012

(XLSX)

Acknowledgments

We thank Gianluca Marucci and Stefania Angelucci for help during the course of this work.

Author Contributions

Conceived and designed the experiments: SA PS MT. Performed the experiments: PS MT FDG EE MS RL VE AP MP. Analyzed the data: PS AU DC. Wrote the paper: SA PS MT.

References

  1. 1. Siegel R, Naishadham D, Jemal A (2013) Cancer statistics, 2013. CA Cancer J Clin 63: 11–30.
  2. 2. Louis D, Ohgaki H, Wiestler O, Cavenee W (2007) World Health Organization Classification of Tumours of the Central Nervous System. Lyon: IARC Press.
  3. 3. Louis DN, Ohgaki H, Wiestler OD, Cavenee WK, Burger PC, et al. (2007) The 2007 WHO classification of tumours of the central nervous system. Acta Neuropathol 114: 97–109.
  4. 4. De Vita VT, Lawrence TS, Rosenberg SA (2008) De Vita, Hellman & Rosenberg's Cancer: Principles & Practice of Oncology, 8th Edition; De Vita VT, Hellman S, Rosenberg SA, editors. Philadelphia: Lippincott Williams & Wilkins.
  5. 5. Kleihues P, Louis DN, Scheithauer BW, Rorke LB, Reifenberger G, et al. (2002) The WHO classification of tumors of the nervous system. J Neuropathol Exp Neurol 61: 215–225 discussion 226–219.
  6. 6. Aldape K, Simmons ML, Davis RL, Miike R, Wiencke J, et al. (2000) Discrepancies in diagnoses of neuroepithelial neoplasms: the San Francisco Bay Area Adult Glioma Study. Cancer 88: 2342–2349.
  7. 7. Radner H, Blumcke I, Reifenberger G, Wiestler OD (2002) [The new WHO classification of tumors of the nervous system 2000. Pathology and genetics]. Pathologe 23: 260–283.
  8. 8. Ohgaki H, Dessen P, Jourde B, Horstmann S, Nishikawa T, et al. (2004) Genetic pathways to glioblastoma: a population-based study. Cancer Res 64: 6892–6899.
  9. 9. Rasheed BK, Wiltshire RN, Bigner SH, Bigner DD (1999) Molecular pathogenesis of malignant gliomas. CurrOpinOncol 11: 162–167.
  10. 10. Rich JN, Bigner DD (2004) Development of novel targeted therapies in the treatment of malignant glioma. Nat Rev Drug Discov 3: 430–446.
  11. 11. Yang L, Luquette LJ, Gehlenborg N, Xi R, Haseley PS, et al. (2013) Diverse mechanisms of somatic structural variations in human cancer genomes. Cell 153: 919–929.
  12. 12. Frattini V, Trifonov V, Chan JM, Castano A, Lia M, et al.. (2013) The integrated landscape of driver genomic alterations in glioblastoma. Nat Genet.
  13. 13. Chen J, McKay RM, Parada LF (2012) Malignant glioma: lessons from genomics, mouse models, and stem cells. Cell 149: 36–47.
  14. 14. Bleeker FE, Molenaar RJ, Leenstra S (2012) Recent advances in the molecular understanding of glioblastoma. J Neurooncol 108: 11–27.
  15. 15. Chumbalkar V, Sawaya R, Bogler O (2008) Proteomics: the new frontier also for brain tumor research. Curr Probl Cancer 32: 143–154.
  16. 16. Hanash SM, Bobek MP, Rickman DS, Williams T, Rouillard JM, et al. (2002) Integrating cancer genomics and proteomics in the post-genome era. Proteomics 2: 69–75.
  17. 17. Odreman F, Vindigni M, Gonzales ML, Niccolini B, Candiano G, et al. (2005) Proteomic studies on low- and high-grade human brain astrocytomas. J Proteome Res 4: 698–708.
  18. 18. Okamoto H, Li J, Glasker S, Vortmeyer AO, Jaffe H, et al. (2007) Proteomic Comparison of Oligodendrogliomas with and without 1pLOH. Cancer Biol Ther 6: 391–396.
  19. 19. Niclou SP, Fack F, Rajcevic U (2010) Glioma proteomics: status and perspectives. J Proteomics 73: 1823–1838.
  20. 20. Gimenez M, Souza VC, Izumi C, Barbieri MR, Chammas R, et al. (2010) Proteomic analysis of low- to high-grade astrocytomas reveals an alteration of the expression level of raf kinase inhibitor protein and nucleophosmin. Proteomics 10: 2812–2821.
  21. 21. Marucci G, Morandi L, Magrini E, Farnedi A, Franceschi E, et al. (2008) Gene expression profiling in glioblastoma and immunohistochemical evaluation of IGFBP-2 and CDC20. Virchows Arch 453: 599–609.
  22. 22. Bjellqvist B, Pasquali C, Ravier F, Sanchez JC, Hochstrasser D (1993) A nonlinear wide-range immobilized pH gradient for two-dimensional electrophoresis and its definition in a relevant pH scale. Electrophoresis 14: 1357–1365.
  23. 23. Hochstrasser D, Harrington M, Hochstrasser A, Miller M, Merril C (1988) Methods for increasing the resolution of two-dimensional protein electrophoresis. Anal Biochem 173: 424–435.
  24. 24. McLeod J, Curtis N, Lewis HD, Good MA, Fagan MJ, et al. (2009) Gamma-secretase-dependent cleavage of amyloid precursor protein regulates osteoblast behavior. FASEB J 23: 2942–2955.
  25. 25. Xia J, Psychogios N, Young N, Wishart DS (2009) MetaboAnalyst: a web server for metabolomic data analysis and interpretation. Nucleic Acids Res 37: W652–660.
  26. 26. Xia J, Wishart DS (2011) Web-based inference of biological patterns, functions and pathways from metabolomic data using MetaboAnalyst. Nat Protoc 6: 743–760.
  27. 27. Xia J, Mandal R, Sinelnikov IV, Broadhurst D, Wishart DS (2012) MetaboAnalyst 2.0–a comprehensive server for metabolomic data analysis. Nucleic Acids Res.
  28. 28. Eriksson L, Antti H, Gottfries J, Holmes E, Johansson E, et al. (2004) Using chemometrics for navigating in the large data sets of genomics, proteomics, and metabonomics (gpm). Anal Bioanal Chem 380: 419–429.
  29. 29. Arndt D, Xia J, Liu Y, Zhou Y, Guo AC, et al. (2012) METAGENassist: a comprehensive web server for comparative metagenomics. Nucleic Acids Res 40: W88–95.
  30. 30. van den Berg RA, Hoefsloot HC, Westerhuis JA, Smilde AK, van der Werf MJ (2006) Centering, scaling, and transformations: improving the biological information content of metabolomics data. BMC Genomics 7: 142.
  31. 31. Xia J, Wishart DS (2002) Metabolomic Data Processing, Analysis, and Interpretation Using MetaboAnalyst. Current Protocols in Bioinformatics: John Wiley & Sons, Inc.
  32. 32. Janes KA, Yaffe MB (2006) Data-driven modelling of signal-transduction networks. Nat Rev Mol Cell Biol 7: 820–828.
  33. 33. Westerhuis JA, Hoefsloot HCJ, Smit S, Vis DJ, Smilde AK, et al. (2008) Assessment of PLSDA cross validation. Metabolomics 4: 81–89.
  34. 34. Vergara D, Simeone P, del Boccio P, Toto C, Pieragostino D, et al. (2013) Comparative proteome profiling of breast tumor cell lines by gel electrophoresis and mass spectrometry reveals an epithelial mesenchymal transition associated protein signature. Mol Biosyst 9: 1127–1138.
  35. 35. Wiklund S, Nilsson D, Eriksson L, Sjöström M, Wold S, et al. (2007) A randomization test for PLS component selection. Journal of Chemometrics 21: 427–439.
  36. 36. Bjerrum JT, Nielsen OH, Hao F, Tang H, Nicholson JK, et al. (2010) Metabonomics in ulcerative colitis: diagnostics, biomarker identification, and insight into the pathophysiology. J Proteome Res 9: 954–962.
  37. 37. Chan EC, Pasikanti KK, Nicholson JK (2011) Global urinary metabolic profiling procedures using gas chromatography-mass spectrometry. Nat Protoc 6: 1483–1499.
  38. 38. Sun H (2004) A universal molecular descriptor system for prediction of logP, logS, logBB, and absorption. J Chem Inf Comput Sci 44: 748–757.
  39. 39. Eriksson L, Umetrics AB (2006) Multi- and Megavariate Data Analysis, Part 1, Basic Principles and Applications: Umetrics AB.
  40. 40. Berglund L, Bjorling E, Oksvold P, Fagerberg L, Asplund A, et al. (2008) A genecentric Human Protein Atlas for expression profiles based on antibodies. Mol Cell Proteomics 7: 2019–2027.
  41. 41. Uhlen M, Oksvold P, Fagerberg L, Lundberg E, Jonasson K, et al. (2010) Towards a knowledge-based Human Protein Atlas. Nat Biotechnol 28: 1248–1250.
  42. 42. Levanen B, Wheelock AM (2009) Troubleshooting image analysis in 2DE. Methods Mol Biol 519: 113–129.
  43. 43. Kleno TG, Leonardsen LR, Kjeldal HO, Laursen SM, Jensen ON, et al. (2004) Mechanisms of hydrazine toxicity in rat liver investigated by proteomics and multivariate data analysis. Proteomics 4: 868–880.
  44. 44. Norden B, Broberg P, Lindberg C, Plymoth A (2005) Analysis and understanding of high-dimensionality data by means of multivariate data analysis. Chem Biodivers 2: 1487–1494.
  45. 45. Marengo E, Robotti E, Bobba M, Milli A, Campostrini N, et al. (2008) Application of partial least squares discriminant analysis and variable selection procedures: a 2D-PAGE proteomic study. Anal Bioanal Chem 390: 1327–1342.
  46. 46. Varambally S, Yu J, Laxman B, Rhodes DR, Mehra R, et al. (2005) Integrative genomic and proteomic analysis of prostate cancer reveals signatures of metastatic progression. Cancer Cell 8: 393–406.
  47. 47. Minagawa H, Honda M, Miyazaki K, Tabuse Y, Teramoto R, et al. (2008) Comparative proteomic and transcriptomic profiling of the human hepatocellular carcinoma. Biochem Biophys Res Commun 366: 186–192.
  48. 48. Gygi SP, Rochon Y, Franza BR, Aebersold R (1999) Correlation between protein and mRNA abundance in yeast. Mol Cell Biol 19: 1720–1730.
  49. 49. Guerra E, Trerotola M, Dell' Arciprete R, Bonasera V, Palombo B, et al. (2008) A bi-cistronic CYCLIN D1-TROP2 mRNA chimera demonstrates a novel oncogenic mechanism in human cancer. Cancer Res 68: 8113–8121.
  50. 50. Trerotola M, Cantanelli P, Guerra E, Tripaldi R, Aloisi AL, et al.. (2013) Up-regulation of Trop-2 quantitatively stimulates human cancer growth. Oncogene 32: 222–233.
  51. 51. Guerra E, Trerotola M, Aloisi AL, Tripaldi R, Vacca G, et al. (2013) The Trop-2 signalling network in cancer growth. Oncogene 32: 1594–1600.
  52. 52. Rigamonti D, Bauer JH, De-Fraja C, Conti L, Sipione S, et al. (2000) Wild-type huntingtin protects from apoptosis upstream of caspase-3. J Neurosci 20: 3705–3713.
  53. 53. Pitts A, Dailey K, Newington JT, Chien A, Arseneault R, et al. (2012) Dithiol-based compounds maintain expression of antioxidant protein peroxiredoxin 1 that counteracts toxicity of mutant huntingtin. J Biol Chem 287: 22717–22729.
  54. 54. Erdmann S, Senkel S, Arndt T, Lucas B, Lausen J, et al. (2007) Tissue-specific transcription factor HNF4alpha inhibits cell proliferation and induces apoptosis in the pancreatic INS-1 beta-cell line. Biol Chem 388: 91–106.
  55. 55. Lucas B, Grigo K, Erdmann S, Lausen J, Klein-Hitpass L, et al. (2005) HNF4alpha reduces proliferation of kidney cells and affects genes deregulated in renal cell carcinoma. Oncogene 24: 6418–6431.
  56. 56. Spath GF, Weiss MC (1998) Hepatocyte nuclear factor 4 provokes expression of epithelial marker genes, acting as a morphogen in dedifferentiated hepatoma cells. J Cell Biol 140: 935–946.
  57. 57. Hwang-Verslues WW, Sladek FM (2008) Nuclear receptor hepatocyte nuclear factor 4alpha1 competes with oncoprotein c-Myc for control of the p21/WAF1 promoter. Mol Endocrinol 22: 78–90.
  58. 58. Walesky C, Edwards G, Borude P, Gunewardena S, O'Neil M, et al. (2013) Hepatocyte nuclear factor 4 alpha deletion promotes diethylnitrosamine-induced hepatocellular carcinoma in rodents. Hepatology 57: 2480–2490.
  59. 59. Kohn KW (1999) Molecular interaction map of the mammalian cell cycle control and DNA repair systems. Mol Biol Cell 10: 2703–2734.
  60. 60. Orian A, van Steensel B, Delrow J, Bussemaker HJ, Li L, et al. (2003) Genomic binding by the Drosophila Myc, Max, Mad/Mnt transcription factor network. Genes Dev 17: 1101–1114.
  61. 61. Meyer N, Penn LZ (2008) Reflecting on 25 years with MYC. Nat Rev Cancer 8: 976–990.
  62. 62. Sebastian C, Zwaans BM, Silberman DM, Gymrek M, Goren A, et al. (2013) The histone deacetylase SIRT6 is a tumor suppressor that controls cancer metabolism. Cell 151: 1185–1199.
  63. 63. Morrish F, Isern N, Sadilek M, Jeffrey M, Hockenbery DM (2009) c-Myc activates multiple metabolic networks to generate substrates for cell-cycle entry. Oncogene 28: 2485–2491.
  64. 64. Masui K, Tanaka K, Akhavan D, Babic I, Gini B, et al. (2013) mTOR Complex 2 Controls Glycolytic Metabolism in Glioblastoma through FoxO Acetylation and Upregulation of c-Myc. Cell Metab 18: 726–739.
  65. 65. Marie SK, Shinjo SM (2011) Metabolism and brain cancer. Clinics (Sao Paulo) 66 Suppl 133–43.
  66. 66. Furnari FB, Fenton T, Bachoo RM, Mukasa A, Stommel JM, et al. (2007) Malignant astrocytic glioma: genetics, biology, and paths to treatment. Genes Dev 21: 2683–2710.
  67. 67. Ziegler DS, Kung AL, Kieran MW (2008) Anti-apoptosis mechanisms in malignant gliomas. J Clin Oncol 26: 493–500.
  68. 68. Seyfried TN, Sanderson TM, El-Abbadi MM, McGowan R, Mukherjee P (2003) Role of glucose and ketone bodies in the metabolic control of experimental brain cancer. Br J Cancer 89: 1375–1382.
  69. 69. Galarraga J, Loreck DJ, Graham JF, DeLaPaz RL, Smith BH, et al. (1986) Glucose metabolism in human gliomas: correspondence of in situ and in vitro metabolic rates and altered energy metabolism. Metab Brain Dis 1: 279–291.
  70. 70. Liang S, Shen G, Liu Q, Xu Y, Zhou L, et al. (2009) Isoform-specific expression and characterization of 14-3-3 proteins in human glioma tissues discovered by stable isotope labeling with amino acids in cell culture-based proteomic analysis. Proteomics Clin Appl 3: 743–753.
  71. 71. Yang X, Cao W, Zhou J, Zhang W, Zhang X, et al. (2011) 14-3-3zeta positive expression is associated with a poor prognosis in patients with glioblastoma. Neurosurgery 68: 932–938 discussion 938.
  72. 72. Arora S, Matta A, Shukla NK, Deo SV, Ralhan R (2005) Identification of differentially expressed genes in oral squamous cell carcinoma. Mol Carcinog 42: 97–108.
  73. 73. Jang JS, Cho HY, Lee YJ, Ha WS, Kim HW (2004) The differential proteome profile of stomach cancer: identification of the biomarker candidates. Oncol Res 14: 491–499.
  74. 74. Somiari RI, Somiari S, Russell S, Shriver CD (2005) Proteomics of breast carcinoma. J Chromatogr B Analyt Technol Biomed Life Sci 815: 215–225.
  75. 75. Huber E, Vlasny D, Jeckel S, Stubenrauch F, Iftner T (2004) Gene profiling of cottontail rabbit papillomavirus-induced carcinomas identifies upregulated genes directly Involved in stroma invasion as shown by small interfering RNA-mediated gene silencing. J Virol 78: 7478–7489.
  76. 76. Parsons DW, Jones S, Zhang X, Lin JC, Leary RJ, et al. (2008) An integrated genomic analysis of human glioblastoma multiforme. Science 321: 1807–1812.
  77. 77. Yan H, Parsons DW, Jin G, McLendon R, Rasheed BA, et al. (2009) IDH1 and IDH2 mutations in gliomas. N Engl J Med 360: 765–773.
  78. 78. Krell D, Assoku M, Galloway M, Mulholland P, Tomlinson I, et al. (2011) Screen for IDH1, IDH2, IDH3, D2HGDH and L2HGDH mutations in glioblastoma. PLoS One 6: e19868.
  79. 79. Rahman-Roblick R, Hellman U, Becker S, Bader FG, Auer G, et al. (2008) Proteomic identification of p53-dependent protein phosphorylation. Oncogene 27: 4854–4859.
  80. 80. Feng Z, Jin S, Zupnick A, Hoh J, de Stanchina E, et al. (2006) p53 tumor suppressor protein regulates the levels of huntingtin gene expression. Oncogene 25: 1–7.
  81. 81. Agrawal P, Yu K, Salomon AR, Sedivy JM (2010) Proteomic profiling of Myc-associated proteins. Cell Cycle 9: 4908–4921.
  82. 82. Maeda Y, Seidel SD, Wei G, Liu X, Sladek FM (2002) Repression of hepatocyte nuclear factor 4alpha tumor suppressor p53: involvement of the ligand-binding domain and histone deacetylase activity. Mol Endocrinol 16: 402–410.
  83. 83. Waterman MJ, Stavridi ES, Waterman JL, Halazonetis TD (1998) ATM-dependent activation of p53 involves dephosphorylation and association with 14-3-3 proteins. Nat Genet 19: 175–178.
  84. 84. Bradley SV, Holland EC, Liu GY, Thomas D, Hyun TS, et al. (2007) Huntingtin interacting protein 1 is a novel brain tumor marker that associates with epidermal growth factor receptor. Cancer Res 67: 3609–3615.
  85. 85. Chou YT, Lin HH, Lien YC, Wang YH, Hong CF, et al. (2010) EGFR promotes lung tumorigenesis by activating miR-7 through a Ras/ERK/Myc pathway that targets the Ets2 transcriptional repressor ERF. Cancer Res 70: 8822–8831.
  86. 86. Oksvold MP, Huitfeldt HS, Langdon WY (2004) Identification of 14-3-3zeta as an EGF receptor interacting protein. FEBS Lett 569: 207–210.
  87. 87. Martinho O, Granja S, Jaraquemada T, Caeiro C, Miranda-Goncalves V, et al. (2012) Downregulation of RKIP is associated with poor outcome and malignant progression in gliomas. PLoS One 7: e30769.