Abstract
Purpose: There is an urgent need for therapies that will reduce the mortality of patients with bone metastasis. In this study, we profiled the protein signal pathway networks of the human bone metastasis microenvironment. The goal was to identify sets of interacting proteins that correlate with survival time following the first diagnosis of bone metastasis.
Experimental Design: Using Reverse Phase Protein Microarray technology, we measured the expression of 88 end points in the bone microenvironment of 159 bone metastasis tissue samples derived from patients with primary carcinomas and sarcomas.
Results: Metastases originating from different primary tumors showed similar levels of cell signaling across tissue types for the majority of proteins analyzed, suggesting that the bone microenvironment strongly influences the metastatic tumor signaling profiles. In a training set (72 samples), TNF receptor 1, alone (P = 0.0013) or combined with serotonin (P = 0.0004), TNFα (P = 0.0214), and RANK (P = 0.0226), was associated with poor survival, regardless of the primary tumor of origin. Results were confirmed by (i) analysis of an independent validation set (71 samples) and (ii) independent bioinformatic analysis using a support vector machine learning model. Spearman rho analysis revealed a highly significant number of interactions intersecting with ERα S118, serotonin, TNFα, RANKL, and matrix metalloproteinase in the bone metastasis signaling network, regardless of the primary tumor. The interaction network pattern was significantly different in the short versus long survivors.
Conclusions: TNF receptor 1 and neuroendocrine-regulated protein signal pathways seem to play an important role in bone metastasis and may constitute a novel drug-targetable mechanism of seed-soil cross talk in bone metastasis. Clin Cancer Res; 19(9); 2473–85. ©2013 AACR.
Translational Relevance
Current therapies do not significantly extend life for most patients suffering with bone metastasis. Following a diagnosis of bone metastasis, the molecular determinants of patient survival are largely unknown. Do prognostic indicators exist that predict survival independent of primary tumor type? We addressed these questions by applying protein array technology to map the network of 88 signal pathway proteins in bone metastasis specimens from a study set of 159 cases spanning eight types of primary tumors. We identified, and independently validated, dramatic differences in the signaling networks linked to TNF receptor 1 (TNFR1) and the neurohormone serotonin that strongly correlated with survival independent of primary tumor site. These insights provide strategies for treating bone metastasis by combination therapy targeting TNFR1 and serotonin uptake or production in bone.
Introduction
There is an urgent need to develop novel and reliable molecular targets for the prognosis and treatment of bone metastasis, a major cause of cancer-associated pain, fracture, spinal cord compression, and mortality. Although a great deal is known about how metastatic tumor cells influence osteoblasts and osteoclasts to induce bone destruction (1–3), this has not led to therapies that significantly extend life for patients suffering with bone metastasis (4–8). It is unknown why some patients with bone metastasis follow a steep downhill course, whereas other patients survive for many years. What are the key molecules determinants of bone metastasis clinical aggressiveness? Primary tumors of diverse histologic origin can be source of bone tumors. It is unclear whether primary carcinomas versus sarcomas use different bone colonization mechanisms. Are the bone metastasis aggressiveness determinants independent of the primary tumor origin? In our study, we have attempted to address these critical questions by analyzing a unique study set of 159 fresh-frozen bone metastasis specimens collected at the time of first diagnosis of bone metastasis for patients previously diagnosed with a variety of carcinomas and sarcomas. We conducted protein signal network analysis of 88 proteins/phosphoproteins known to be directly/indirectly influenced by the bone remodeling ligands tumor necrosis factor α (TNFα) and receptor activator of NF-κβ ligand (RANKL). TNFα and RANKL are abundant in sites of inflammatory bone erosion (9, 10). Within bone metastasis, the liberation of RANKL, which binds to RANK present on osteoclast precursors, is associated with osteolysis or bone matrix disorganization (11, 12). Using Reverse Phase Protein Microarray (RPMA) technology, we evaluated pathway interconnections in the bone metastasis microenvironment correlating with survival time after first diagnosis of bone metastasis, comparing both carcinoma and sarcoma primary tumors. Results revealed a previously unknown role for neuroendocrine signaling and TNF receptor 1 (TNFR1) as major determinants of survival.
Materials and Methods
The biospecimens used for this study were snap-frozen portions of bone metastasis collected over a 10-year period from a variety of primary tumors, including some rare categories of primary tumors (Table 1). Bone metastases specimens (159) from a spectrum of primary carcinomas and primary sarcomas (Table 1 and Supplementary Table S1) were collected at the Istituto Ortopedico Rizzoli (Bologna, Italy) under patient informed consent. The specimens were collected at the time of first work-up of the metastatic disease and before treatment of the metastatic disease. Tissue was excised in the surgical suite and transported refrigerated to the pathology laboratory. A board-certified pathologist carried out gross examination of each tissue sample and provided diseased tissue not required for diagnosis. Bone metastasis specimens were snap-frozen, maintained at −80°C, and thawed once for analysis.
Sample set summary.
Frozen bone metastasis samples were divided into 2 pieces. One piece was thawed/fixed in 10% formalin and decalcified in formic/nitric acid and paraffin embedded (FFPE) to evaluate the tumor cell percentage. The second piece was lysed using Adaptive Focused Acoustic (AFA) technology (Covaris) to quantify the levels of 88 proteins by RPMA.
We have recently developed a novel fixation method for decalcification of bone and preservation of proteins (13). This fixative was used for immunohistochemical examination of serotonin in example bone metastasis cases material from recently diagnosed cases that were not part of the 159 bone metastasis cases described above with 10-year follow-up.
Survival under standard of care therapy was defined as the number of months from the first diagnosis of bone metastasis (Supplementary Table S1).
Tissue homogenization and protein extraction
One part of each tissue sample was weighed while frozen, then submerged in liquid nitrogen, and completely pulverized using CryoPrep technology (Covaris). Pulverized samples were homogenized in protein extraction buffer—45% T-PER (Thermo Scientific), 45% Novex Tris-Glycine SDS Sample Buffer 2× (Invitrogen), and 10% TCEP (Thermo Scientific)—using AFA technology (Covaris; 20% duty factor, 275 pick incident power, 200 cycles per burst, continuous degassing, 90 seconds). After lysis, the samples were boiled for 7 minutes at 100°C.
Reverse Phase Protein Microarrays
RPMA was conducted as described previously (14–17). Whole tissue lysates were printed in triplicate on ONCYTE Avid nitrocellulose film-slides (GRACE Bio-Labs) using an Aushon 2470 arrayer equipped with 185 μm pins (Aushon Biosystems). Each lysate was printed in a 2-fold dilution curve representing undiluted lysate, 1:2, and 1:4 dilutions. Control lysates were printed in a 2-fold dilution curve. All RPMAs were baked for 2 hours at 80°C to allow DNA fixation and then stored with desiccant at −20°C.
RPMA immunostaining
RPMA slides were blocked (I-Block, Applied Biosystems) for 2 hours before immunostaining. Immunostaining was conducted on a Dako Autostainer [Catalyzed Signal Amplification (CSA) kit, Dako]. Each slide was incubated with a single primary antibody at room temperature for 30 minutes (Supplementary Table S2). The negative control slide was incubated with antibody diluent (Dako). For each immunostaining run, 1 slide was incubated with anti ss-DNA antibody (1:15,000; IBL International GmbH). Secondary antibody was goat antirabbit or rabbit antigoat IgG heavy+light (1:7,500; Vector Laboratories), or rabbit antimouse IgG (1:10; Dako), and amplified via horseradish peroxidase-mediated biotinyl tyramide with chromogenic detection (diaminobenzidine; Dako).
Antibody validation
Primary antibodies were validated before use by Western blotting with complex cellular lysates (cell line or human tissue lysates). Antibodies were valid if a single band at the correct molecular weight comprised greater than 80% of the lane (Supplementary Fig. S1A). Specificity of antibodies against analytes too small to be detected by Western blot analysis was verified by peptide competition (Supplementary Fig. S1C).
Image acquisition and data analysis
RPMA slides were scanned on a UMAX 2100XL flatbed scanner (white balance 255, black 0, middle tone 1.37, 1800 dpi, 14 bit). Spot intensity was analyzed by Image Quant v5.2 software (Molecular Dynamics). Data reduction was carried out with a VBA Excel macro, RPMA Analysis Suite (18). To normalize data, the relative intensity for each protein spot was divided by the ssDNA relative intensity for the corresponding spot (19).
Statistical methods
The Kruskal–Wallis 1-way analysis of variance (GraphPad Prism v5.03; GraphPad Software, Inc.) was applied to the 5 groups of bone metastasis specimens from different primary tumors of origin, using the median of the intensities for each group and to compare TNFR1 expression levels across the 4-patient groups identified by the survival quartiles. The Mann–Whitney t test on medians (GraphPad) was used to compare (i) protein expression levels in metastases colonized from primary sarcomas with metastases derived from primary carcinomas, (ii) TNFR1 expression levels in long-term and short-term survivors, and (iii) sample age in patients with high and low TNFR1 expression levels. P ≤ 0.05 was considered significant.
The Kaplan–Meier product limit estimator (GraphPad) was applied to bone metastasis survival data. Patients were divided into 2 groups as follows (i) above and below the median expression intensities for each endpoint; (ii) 25th/75th percentile of TNFR1 expression; (iii) the presence/absence of a specific metastatic site; (iv) males/females; (v) above/below the median of the sample age. Patients were divided into 8 groups based on the primary tumor origin (head–neck tumors, urogenital tumors, gastrointestinal tumors, breast cancer, lung cancer, melanoma, unknown carcinomas, and sarcomas) to investigate the influence of the primary tumor site on survival. Noninformative censoring was applied. The log-rank nonparametric test (GraphPad) was used to compare the survival distributions of the patient groups. The Spearman rank correlation coefficient, ρ, (JMP 5.1.2; SAS Institute Inc.) was calculated for each protein pair in the RPMA quantitative expression profiles of the 143 bone metastasis specimens with known survival data. P ≥ 0.75 with P ≤ 0.01 was considered significant. A Cox proportional hazards model (The SAS System; SAS Institute Inc.) was applied to calculate the HR of the continuous variables: age at diagnosis, bone metastasis number, total metastasis number, and TNFR1 expression in patients with known survival data. The Breslow estimate of the baseline hazard function was used to obtain the Cox partial likelihood. Variables with P ≤ 0.2 at the univariate analysis were used for the multivariate analysis.
Independent validation of protein analytes that correlated with survival after bone metastasis diagnosis
Analyte subsets correlating with survival were validated using 2 different methods. First, the highly correlated analytes in a training set were independently validated in a separate validation set (Supplementary Table S1). The training (n = 72) and validation sets (n = 71) were randomly selected among the 143 specimens with known survival data. The primary tumor categories were represented in similar proportions for the 2 sets (Supplementary Table S1). After unblinding the training set data, the analytes strongly correlated (P < 0.05) with survival were assessed in the validation set.
Second, the data were independently analyzed by a machine learning classification method (see below) yielding a ranking of the 88 end points based on their contribution to survival.
Machine learning classification methods for validation
Classification models were built using Support Vector Machine (SVM) algorithm in the Weka (Waikato environment for knowledge analysis) suite of machine learning tools (20). SVM uses a kernel function to nonlinearly map training set examples into a higher dimensional feature space, to build an optimal separating hyperplane, providing maximal separation margins between set examples from 2 differing classes and corresponds to a nonlinear decision boundary in the original space. Two cross-validation (CV) methods were used to examine the effectiveness of classifiers: 10-fold cross validation (10-CV) and jackknife test (21). A stratified 10-CV testing procedure was used to generate receiver operating characteristic (ROC) curves. The 10-CV approach randomized instances from the training set into 10 equally sized subsets and used each subset as a test set after a SVM classifier was trained with the remaining 9 subsets combined. Stratification ensured that the survival length class proportions in the original training set were maintained in each of the 10 subsets. A class prediction was obtained for every sample in the original training set using generic class labels positive (P) and negative (N). The overall accuracy (Acc) was computed as
Acc = (TP+TN)/(TP+TN+FP+FN),
where TP = true positive, TN = true negative, FP = false positive, FN = false negative.
Area under the ROC curve (AUC), a measure of the discriminatory capability of the classifier (22), was calculated. An AUC value near 0.5 is suggestive of a model that does not perform better than random guessing, whereas a value of 1.0 indicates a perfect classifier. Individual end points were ranked using the SVM classifier (23) as implemented in the SVMAttributeEval attribute evaluator program in Weka (20). In this approach, the attributes are ranked by the square of the weight assigned by the SVM.
Results
Microenvironment “soil” dominates proteomic signaling in bone metastases derived from carcinomas and sarcomas
We evaluated the relative contribution of the primary tumor site to the proteomic signaling in bone metastases to assess the influence of seed (primary tumor) versus soil (bone) on the signaling network. In the original set of 159 samples, the 5 most numerous primary carcinomas were: kidney, colorectal, breast, lung, and melanoma. The Kruskal–Wallis 1-way analysis of variance was used to compare the expression of 88 end points across the 5 groups. The expression of 16 proteins (BMP4, CREB S133, c-Src family Y416, HSP90, IGF1Rβ, IRS1 S612, LRP6, Prolactin Receptor, RANK, Rb S608, Rb S780, Shc Y317, STAT5 Y694, STAT6 Y641, TLR3, and Vimentin) was significantly different across the 5 primary tumor groups (P value, respectively, 0.001, 0.0350, 0.0122, 0.0226, 0.0176, 0.0260, <0.0001, 0.0496, 0.0484, 0.0350, 0.0369, 0.0137, 0.0131, 0.0061, 0.0088, and 0.0036; Supplementary Table S3). Seventy two end points showed no statistically significant differences in expression level between the 5 groups of metastases.
Sarcomas are derived from mesenchymal cells, whereas carcinomas arise from epithelial cells (24). As such, it would be expected that their bone metastasis cell signaling patterns may be somewhat different, but it is unknown whether their respective metastases would be more similar to the cellular origin (seed) or the metastasis microenvironment (soil). The Mann–Whitney U test was used to compare the protein signaling intensities of the metastases derived from primary sarcomas to the group of metastases derived from primary carcinomas. The expression of 38/88 proteins [Akt T308, androgen receptor (AR), β-arrestin, β-catenin S33/34/T41, BMP4, CaM kinase II, CREB, c-Src, DKK1, E-cadherin, Erk 1/2, ERα S118, Ezrin T567, Ezrin Y353, FAK Y576/577, HSP27, HSP70, IGF1R Y1135, IkBα S32/36, IL-10, insulin receptor, MMP14, mTOR, mTOR S2448, p53 S15, PDGFRβ, PDGFRβ Y716, PDGFRβ Y751, PGR S190, Prolactin Receptor, Rb S780, SAPK/JNK, STAT2 Y690, STAT5 Y694, TIMP2, TLR3, TLR9, and TNFR1] was different in these 2 classes (P–value, respectively, <0.0001, 0.0228, 0.0085, 0.0061, 0.0005, 0.0030, 0.0029, 0.0319, 0.0413, <0.0001, 0.0402, <0.0001, 0.0077, 0.0197, 0.0034, 0.0030, <0.0001, 0.0365, 0.0019, 0.0024, <0.0001, 0.0119, 0.0402, 0.0157, 0.0005, 0.0327, <0.0001, 0.0274, <0.0001, <0.0001, 0.0026, <0.0001, 0.0041, 0.0190, 0.0403, 0.0229, 0.0129, 0.0164; Supplementary Table S3). The expression level of 50 proteins was not statistically different between metastases arising from primary sarcomas compared with primary carcinomas. These data suggest that sarcomas and carcinomas influence their bone metastasis signaling pathways differently, with a generally higher expression of hormone and growth factor-driven pathways and a lower expression of mTOR pathway proteins in sarcomas.
Bone metastasis proteomic networks differ between long and short survivors
Spearman ρ rank correlation analyses for the 88 proteins analyzed revealed that 25 of the 2,556 protein pairs were strongly correlated (P > 0.75) across all samples (Supplementary Table S4). The data set was successively subdivided into quartiles by patients' survival. Top (survival ≥38 months) and bottom (survival ≤5 months) quartiles were compared with each other to identify differences in protein network interconnections between long and short survivors. Fifty-three of the 2,556 protein pairs were highly correlated in the long survivor group (Fig. 1, Supplementary Table S5), whereas 81 protein pairs showed a significant Spearman ρ in the short survivor group (Fig. 1, Supplementary Table S6). ERα S118, TNFα, HSP70, IL-6, SAPK/JNK, RANKL, Syndecan, MMP11, MMP14, serotonin, TIMP2, DKK1, Ezrin Y353, and EGFR were strongly linked to each other in the short survivor cohort (Fig. 1). The protein interaction network of metastatic lesions associated with a long survival was strikingly different from the short survival cohort network (Fig. 1). In the short survivor cases, a highly focused network among bone resorption, inflammation and neuroendocrine proteins radiated outward to interconnect with diverse pathways. In the long survivor, the diversity of coordinated strong interconnections was less focused resulting in an overall “weaker” network.
Bone metastasis molecular network in short and long survivors. Interacting proteins comprising the bone metastasis molecular network are strikingly different in short survivors (survival ≤ 5 months) compared with long survivors (survival ≥ 38 months), using Spearman ρ correlation analysis data for the sample set in this study. Correlations with ρ ≥ 0.75 and P ≤ 0.01 were used to build the network graph (Gephi 0.8.1 beta, The Gephi Consortium, Paris, France, www.gephi.org). Each node represents a molecule: the bigger the node, the more the significant correlations relative to the molecule. Each line connecting 2 nodes represents a significant correlation between the nodes: the thicker the line, the higher the Spearman ρ correlation. Proteins are grouped on the basis of Spearman ρ values and the number of connections among a group of nodes: strongly correlated nodes are represented closed to each other and with the same color. In short survivors, the molecular network is dominated by a group of interconnections between TNFα, RANKL, ERα S118, and serotonin. Interconnections between nodes are wider and less numerous in long-term survivors and lack a dominant network.
TNFR1 protein in bone metastasis samples is related to patients' survival after bone metastasis diagnosis
We explored which proteins within the TNFα and neuroendocrine pathways were associated with overall survival, in both training (n = 72) and validation sets (n = 71), using Kaplan–Meier estimator analysis. When the median protein analyte level was used as the cut-point, in the training set, TNFR1 was strongly correlated with survival (P = 0.0013; Fig. 2), with a median survival difference of 29 months between high-expressing (survival range 0–78 months) and low-expressing patients (survival range 2–83 months). The importance of TNFR1 was confirmed in the validation set (P = 0.0014) in which TNFR1 was associated with a median survival difference of 15.5 months between high-expressing (survival range 1–58 months) and low-expressing patients (survival range 1–98 months; Fig. 3). TNFR1 levels were compared over the 4 quartiles of the survival time course for all bone metastasis combined. TNFR1 expression in the 143 cases was highest in the shortest survival quartile and declined progressively for longer survivors. The Kruskal–Wallis analysis was significant across the 4 survival quartiles (P = 0.007). The Mann–Whitney t test between the top and bottom quartiles showed a significant difference in TNFR1 expression in the 2 groups (P = 0.0009; Fig. 4A). In all cases combined (n = 143), TNFR1 was correlated with patients' survival (Kaplan–Meier analysis; P ≤ 0.0001), with a median of survival difference of 22 months between high-expressing (survival range 0–78 months) and low-expressing patients (survival range 2–98 months; Fig. 4B). Kaplan–Meier analysis applied to the groups of patients in the 25th and 75th quartiles of TNFR1 expression showed a significant difference between the group survivals (P = 0.0001) and a median of survival difference of 29 months between high-expressing (survival range 1–39 months) and low-expressing patients (survival range 3–98 months; Fig. 4C). These data support a general functional relationship between TNFR1 levels and survival.
TNFR1 and bone regulating proteins correlate with patient's survival in the training set. Graphs representing Kaplan–Meier survival curves for TNFR1 and combination of TNFR1 with 1 or more proteins regulating bone metabolism in the training set. The median of the expressions of the end point (or combination of end points) was used as the cut point to create the 2 groups, high-expressing patients (high) and low-expressing patients (low). Combination of end points was achieved by summing together the end point expression values. Survival was defined as the number of months each patient survived following the first diagnosis of bone metastasis.
Correlation of TNFR1 and bone regulating proteins with patient's survival is validated in an independent set. Graphs represent Kaplan–Meier survival curves for TNFR1 and combination of TNFR1 with 1 or more proteins regulating bone metabolism in the validation set. The median of the expressions of the end point (or combination of end points) was used as the cut-point to create the 2 groups, high expressing patients (high) and low expressing patients (low).Combination of end points was achieved by adding together the end point expression values. Survival was defined as the number of months each patient survived following the first diagnosis of bone metastasis.
TNFR1 expression inversely correlates with survival over the full range of the survival curve over all 143 samples including all primary tumor groups. In all 4 graphs, survival is defined as the number of months each patient survived following the first diagnosis of bone metastasis. A, the total set of 143 samples with known survival data was subdivided into 4 quartiles based on survival time. The mean and median TNFR1 expression in the 4 quartiles descends over the shortest to the longest survivors. The Mann–Whitney test was applied to medians of the lowest (x ≤ 5) and highest (x ≥ 38) survival quartiles. The Kruskal–Wallis test was conducted on the medians of 4 groups of samples. B, this graph depicts the Kaplan–Meier survival curves in the total set of 143 samples, using the median of TNFR1 expressions as the cut point to create the high expressing (high) and the low expressing (low) patients' groups. C, Kaplan–Meier survival curves for the 25th and 75th quartile of TNFR1 expression in the 143 sample set. D, Kaplan–Meier curves showing the survival distribution in each primary tumor origin group of patients (color lines) and in the total set of 143 samples (black line).
TNFR1 levels predict survival independent of common prognostic factors
TNFR1 has not previously been recognized as a marker for bone metastasis survival, therefore we evaluated whether or not TNFR1 was associated with a true influence on patients' survival, independent of the contribution by other clinical variables. The possible correlation between patients' survival following the diagnosis of bone metastasis and the site of the metastatic tumor's primary origin was investigated. The data set was divided into 8 groups: head–neck tumors, urogenital tumors, gastrointestinal tumors, breast cancer, lung cancer, melanoma, unknown carcinomas, and sarcomas. The Kaplan–Meier analysis showed a significant correlation (P = 0.0118) between survival and primary tumor site (Fig. 4D). However, TNFR1 expression levels were unrelated to the primary tumor site (P > 0.05), overruling the hypothesis of a codependence between primary site and TNFR1. Nonetheless, the median TNFR1 level was higher in bone metastases from sarcomas when compared with those from all types of carcinomas. A possible explanation is that the survival range of patients in the primary sarcoma group (2–42 months) was much shorter than that of patients in the primary carcinoma group (0–98 months).
Age at diagnosis, sex, site of metastasis, number of bone metastases, and number of total metastases at follow-up were investigated for their potential correlation with the patients' survival. None of these parameters was related to survival by Kaplan–Meier analysis (P > 0.05). With univariate analysis, P values for the number of bone metastases and the number of total metastases in all organs were P ≤ 0.2, therefore these 2 variables were included in the multivariate Cox proportional hazard model with the level of TNFR1. TNFR1 expression level was the only contributing variable associated with survival time. With P = 0.0004 and HR = 1.007, the increase of TNFR1 expression by 1 unit reduces the survival by 0.7%.
We also investigated whether the sample specimen age was correlated with TNFR1 expression level. The Mann–Whitney t test applied to the median TNFR1 levels of the 25th and 75th quartiles was not significantly different providing additional support to the finding that TNFR1 was independent of common prognostic factors or specimen collection/specimen age.
TNFR1 and bone modulating proteins/hormones correlate with overall survival
To investigate potential cell signaling interactions between TNFR1 and neuroendocrine hormones that could modulate the bone microenvironment, we used the sum of the expression values of TNFR1 and known modulators of bone for further Kaplan–Meier survival analyses. TNFR1+Serotonin, TNFR1+RANK, TNFR1+TNFα, TNFR1+AR, TNFR1+Serotonin and AR, TNFR1+Serotonin+Ezrin Y353, and TNFR1+Serotonin+Ezrin Y353+ERα S118 correlated significantly with survival both in the training (P values, respectively, 0.0004, 0.0226, 0.0214, 0.0038, 0.0004, 0.0003, 0.0021) and validation sets (P values, respectively, 0.0087, 0.0228, 0.0244, 0.0192, 0.0252, 0.0204, 0.0174; Figs. 2 and 3). These results suggest that TNFR1 signaling couples with the neuroendocrine signaling to promote bone metastasis aggressiveness.
Second independent validation conducted using SVM
Machine learning can provide important insights into possible correlations between quantitative expression profiles and expected survival time. SVM has been shown to work well in a number of bioinformatics applications, including proteomics (25), protein function analysis (26, 27), cancer modeling (28), and others. We tested whether SVM classifiers could discriminate between long and short survivors and validate our results obtained by the Kaplan–Meier analysis. We used the expression level of all 88 end points measured in the 143 samples with known survival data. Each patient was represented by a feature vector, consisting of 89 attributes: 88 end point expression levels and an assigned class label, L for long survivors and S for short survivors. For the machine learning experiments, we created 10 partially overlapping data sets, BM15, BM20, BM25, BM30, BM35, BM40, BM45, BM50, BM55, and BM60, where the number corresponds to the number of samples in the data set from each of the 2 classes. The patients for each data set were selected from the respective ends of the original set, which was sorted in ascending order of survival time [i.e., BM15 set contained 15 patients with the shortest survival time (class S) and 15 patients with the longest survival time (class L)], so that each set was perfectly balanced, (Supplementary Table S7).
We used SVM classification with 2 validation algorithms, 10-fold-CV and jack-knife test [leave-one-out CV (LOOCV)]. LOOCV is the only method that is deterministic and displays low bias (29, 30). Hence, LOOCV has been increasingly used and widely recognized by investigators as an objective choice for examining the performance quality of various predictors (21). In 9 of 10 sets (BM15–BM55), SVM classifier was carried out with 62% to 76% accuracy, calculated by both 10-CV and LOOCV, which strongly suggests that short survivors and long survivors are well separated by protein expression values.
To test the baseline bias, we randomly reshuffled the class labels in several sets: in all the cases, the resulting accuracy closely matched the expected 50% from random guessing. Because all the sets were balanced, accuracy rates and AUC strongly correlated. AUCs calculated using LOOCV displayed more stable behavior over the range of set sizes, however, the method performance deteriorated on the BM60 set, where AUC dropped to 0.55 calculated by both validation techniques (Supplementary Fig. S1B). The decrease in performance was largely due to loss of sensitivity for class L, with more than half (33 of 60) instances misclassified. This was likely due to the fact that BM60 set included patients with a survival close to the median of the survivals. Diminished classification accuracy for the BM60 set suggests that patients in the 23- to 28-month survival range [respectively the L class boundaries for BM60 and BM55 (Supplementary Table S7)] display protein expression patterns characteristic of the class S and that the boundary between class S and L is located near the 28-month point on the survival time scale.
SVM ranking confirms the training set data
To determine the relative contribution of the 88 end point expression profiles to the SVM discriminating power, we ranked the proteins using the attribute evaluator program SVMAttributeEval. Results showed serotonin, TNFR1, AR, ERα S118, and CREB as the top 5 proteins in the ranking, confirming the Kaplan–Meier analysis results. Then we conducted a new classification on a single set (BM45) by SVM and 10-fold CV with the reduced feature vectors, using the top 5 proteins, top 10 (top 5 + Erk1/2, osteopontin, E-cadherin, β-catenin S33/34/T41, Rb S608), top 15 (top 10 + Insulin Receptor, β-catenin, STAT3 S727, c-Src Family Y416, SAPK/JNK T183/Y185), top 20 (top 15 + HSP70, NFkB p65 S536, FADD S194, IRS1 S612, EGFR Y1068), and top 25 (top 20 + RANK, HSP27, SAPK/JNK, PGR S190, vimentin; Table 2). Set BM45 was selected because it was the largest set with the highest AUC (Supplementary Fig. 1B). Although with 5 or less top ranking proteins, the classifier performance was worse than with the full feature vector, using high-ranking proteins in the range of 10 to 25 top selections improved the performance dramatically. The accuracy reached 83% for 20 and 25 top ranking proteins, implying that proteins in the top 25 ranking play a very important role in defining survival time.
Ten-fold cross validation results for the full BM45 set (see Materials and Methods) and 5 subsets with 5, 10, 15, 20, and 25 top-ranked attributes
Discussion
This study aimed to explore the protein signal pathway network of the bone metastasis microenvironment that correlates with patient survival following initial diagnosis of bone metastasis. We measured 88 proteins using RPMA applied to bone metastasis specimens. The analyte levels were evaluated for correlation with survival data within a unique study set of 159 bone metastasis from carcinomas, melanomas, and sarcomas (8 types of primary tumors) over a follow-up of 100 months following the initial diagnosis of bone metastasis. The data provided insights about prognostic markers and therapeutic targets that could be common to all patients with bone metastasis regardless of the primary tumor of origin.
Metastases originating from different primary tumors showed a similar pattern of proteomic cell signaling across tissue types for the majority of 88 inflammation, proliferation, invasion, and adhesion pathway proteins analyzed. Hormone receptors and mTOR pathway expression was similar among metastasis samples from different carcinoma types, whereas bone metastases from sarcomas showed a significantly higher expression of hormone receptors and lower expression of mTOR pathway proteins, compared with those originating from carcinomas (Supplementary Table S3). These data suggest that the bone microenvironment influences the metastatic tumor signaling profiles (seed vs. soil hypothesis), although proteins strongly related to the mesenchymal cell biology show a higher expression in bone metastases from primary sarcomas.
This study revealed several new signaling pathway classes associated with bone metastasis survival. In contrast to previous studies focusing only on the ligand TNFα, we included an evaluation of its receptor TNFR1 and explored the interconnections of this receptor with other pathways. TNFR1 has not been previously recognized as a determinant of bone metastasis prognosis. In the present study, TNFR1 was strongly associated with poor survival regardless of the primary tumor of origin (Figs. 2–4).
Spearman ρ rank comparison analyses were conducted to determine and clarify the network linkages within the bone metastasis microenvironment. Statistically significant correlations between signal protein pairs were different in short and long survivors suggesting the deregulation of cell signaling in the more aggressive metastases (Fig. 1 and Supplementary Tables S5 and S6). Such correlations may reflect direct or indirect interactions between protein pairs within the same cells or between host and tumor cells.
In the short-survivor cohort, the linkages were dominated by interconnections between TNFα and ERα S118, RANKL, serotonin, MMPs, HSP70, IL6, and Ezrin Y353 (Fig. 1). Kaplan–Meier analyses revealed a strong correlation of TNFR1, alone or in combination with RANK and/or serotonin and/or AR and/or TNFα and/or Ezrin Y353, with the patients' survival (Figs. 2 and 3). The strong interconnections between the ligands TNFα and RANKL in short-survivor group point to a concomitant role for their receptors, TNFR1 and RANK, to influence patients' survival. Moreover, a second new insight revealed in this study was the role of neuroendocrine pathways including serotonin and AR in bone metastasis survival. Serotonin and AR proved to be strongly linked to TNF and RANKL pathways.
The role of TNFR1 and serotonin was confirmed in 2 different validation methods. In addition to a traditional training and independent validation set comparison, we used machine learning classification models of survival span as a second independent validation method. The machine learning models confirmed the important role of TNFR1, serotonin, RANK, AR, and ERα S118 as the top ranking end points associated with survival. Including the next highest ranked 20 proteins augmented the discrimination power between short and long survivors, exemplifying the complex microenvironment interactions in bone metastasis (Table 2).
A weakness of our study is that we cannot know the subcellular localization of proteins we have measured in the bone metastasis microenvironment. We are using bone metastasis specimens rapidly frozen at the time of surgery. Consequently, the specimens contained tumor cells, osteoclasts and osteoblasts, inflammatory cells, stroma, blood, and bone matrix. The calcified bone content renders these specimens unsuitable for cryosectioning. For this reason, we pulverized and lysed the specimens directly from the frozen state. Normalization of RPMA proteomic data used cellular DNA content (19) as the best approximation of the total cellular content of the bone metastases and was less biased by blood and noncellular protein content. We can postulate that the elevation of TNFR1 in bone metastasis is taking place at the host bone cell level, but we cannot confirm this because bone decalcification and formalin fixation drastically alters the antigenicity and preservation of membrane antigens. Thus, it would not be meaningful to make quantitative conclusions about the cellular distribution of TNFR1 and serotonin in our frozen bone specimens based on immunohistochemistry of decalcified specimens. We have recently developed a novel fixation chemistry that will decalcify bone while stabilizing tissue phosphoproteins for RPMA and maintaining full diagnostic histomorphology and complete preservation of immunohistochemistry antigens (13). Examination of bone using this fixative revealed anti-serotonin immunohistochemical staining of bone cells consistent with bone remodeling cells (31, 32). In the future, this new preservation chemistry can be used for immunohistochemistry of bone metastasis to study the cellular localization and expression levels of analytes correlating with survival in this study.
These data support a previously unappreciated, but not unexpected, role for TNFR1 as a major driver of bone metastasis pathogenesis. TNFα has been linked to chronic inflammation and tumor progression. (33–35) TNFα upregulates RANK expression through the TNFR1 signal pathway. Therefore, elevation of TNFR1 would be expected to mobilize osteoclastogenesis and thereby mediate severe bone resorption (36). In keeping with this interpretation, a recent study reported that TNFα inhibition reduces cell motility and bone metastases in a metastatic breast cancer cell line, MDA-MB-231 (37). These data support a role for TNFα/TNFR1 in cancer progression and suggest a novel therapeutic strategy for bone metastasis. Selective TNFR1 inhibitors, such as Atrosab (Baliopharm), are currently being evaluated for the treatment of arthritis. (38–40) The effect of this drug in animal models of bone metastasis should be investigated.
These data also reveal an unexpected role for the neuroendocrine pathways in bone metastasis. Endocrine-regulated pathways may modulate, or link together, survival, motility and osteoclast differentiation and activation. Peripheral serotonin levels are known to regulate TNFα and other cytokines, and serotonin is known to regulate normal bone turnover (31, 41, 42), thus serotonin may also constitute a novel mechanism of seed-soil cross talk in the bone metastasis microenvironment.
Although it is difficult to directly investigate the molecular mechanisms of bone metastasis in archival frozen human bone samples, we believe that we have gained mechanistic insights in the present study in 2 ways: (i) correlation of key protein signal pathway proteins with survival post diagnosis of bone metastasis, that are confirmed by independent validation, and (ii) Spearman ρ network analysis of signaling proteins that indicates a marked difference in the crosstalk between signaling pathway nodes in the bone microenvironment associated with survival. These 2 types of analysis complement each other and provide evidence for a mechanistic interplay between neuroendocrine pathways and TNFR1. Preliminary data show that TNFR1 expression levels increase by 2-fold after 24 hours of treatment with serotonin or Tianeptine, a compound that increases serotonin uptake, compared with controls, in the U266 cell line which expresses TNFR1. These data are in keeping with a feedback loop between serotonin and TNFR1. Moreover, a recent study revealed that osteoclast precursors are capable of producing serotonin and that the availability of serotonin locally stimulates osteoclastogenesis (32). Drugs that modulate the level of peripheral serotonin or its action in the bone microenvironment may constitute a possible strategy to treat bone metastasis and improve survival (43).
We hypothesize that TNFR1 and serotonin participate in the cross talk between the tumor cells and the bone microenvironment in a manner that is independent of the primary tumor of origin. These data support a model of “seed versus soil” in which a common set of signaling pathways permit the metastasis to target and flourish in the bone, regardless of the primary tumor of origin.
We hypothesize that a combined therapy targeting TNFR1 and serotonin could improve the patients' survival through a double feedback mechanism: (i) low levels of peripheral serotonin downregulate TNFα expression, decreasing peripheral serotonin levels; (ii) the blockade of TNFR1 impedes the proinflammatory activity of TNFα, lowering peripheral serotonin, which in turn further downregulates TNFα expression.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: A. Chiechi, E.F. Petricoin, M.S. Benassi, P. Picci, V. Espina, L.A. Liotta
Development of methodology: E.F. Petricoin, I.I. Vaisman, L.A. Liotta
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Chiechi, C. Novello, G. Magagnoli, M.S. Benassi, P. Picci
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A. Chiechi, J. Deng, P. Picci, I.I. Vaisman, L.A. Liotta
Writing, review, and/or revision of the manuscript: A. Chiechi, I.I. Vaisman, V. Espina, L.A. Liotta
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Chiechi, E.F. Petricoin, P. Picci, L.A. Liotta
Study supervision: P. Picci, V. Espina, L.A. Liotta
Grant Support
A. Chiechi was supported by Istituto Ortopedico Rizzoli, George Mason University, and the Italian Istituto Superiore di Sanita' in the framework of the Italy/USA cooperative agreement between the U.S. Department of Health and Human Services, George Mason University, and the Italian Ministry of Public Health. C. Novello was supported by Fondazione Italiana per la Ricerca sul Cancro (FIRC), triennal fellowship “Mario e Valeria Rindi”. G. Magagnoli was supported by Istituto Orotpedico Rizzoli and the Italian Istituto Superiore di Sanita' in the framework of the Italy/USA cooperation agreement between the U.S. Department of Health and Human Services, George Mason University and the Italian Ministry of Public Health. M.S. Benassi and P. Picci were supported by Istituto Orotpedico Rizzoli. I. Vaisman, E.F. Petricoin III, V. Espina, and L.A. Liotta were supported by George Mason University.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
Acknowledgments
The authors thank the IOR pathology group for providing tissue specimens for this study. Sally Ross and Noor Lach-Hab conducted antibody validation. Daphne Fong assisted with RPMA detection and analysis.
Footnotes
-
Note: Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).
- Received November 2, 2012.
- Revision received February 27, 2013.
- Accepted February 27, 2013.
- ©2013 American Association for Cancer Research.
References
Volume 19, Issue 9
