Pan-Cancer Analyses Reveal Prognostic Value of Osteomimicry Across 20 Solid Cancer Types

Background Osteomimicry of cancer cells had been widely reported in prostate cancer and breast cancer. However, the prognostic value of osteomimicry in various cancer types remained unclear. We hypothesized that osteomimicry would result in remodeling of the tumor microenvironment and was eligible to predict patient prognosis. Methods A comprehensive transcriptomic analysis of the osteomimicry, which was characterized by mRNA expression of SPARC, SPP1, and BGLAP, across 20 solid tumors (7564 patients) using RNA-seq data from The Cancer Genome Atlas (TCGA) was conducted. Samples of each cancer type were classified into subgroups (high vs. low) based on median value of osteomimetic markers, the associations of these markers with clinical outcomes, immune cell infiltration and immune checkpoints expression were explored. Results Each osteomimetic marker harbored prognostic value in the pan-cancer analyses [SPARC: hazard ratio (HR) = 1.10, p = 0.028; SPP1: HR = 1.25, p < 0.001; BGLAP: HR = 1.13, p = 0.005]. Patients with high expression of all the three genes also had significantly unfavorable survival (HR = 1.61, p < 0.0001) compared with those of low expression. Correlation analyses demonstrated that osteomimicry was closely related to tumor purity, dendritic cells (DC) infiltration and expression of immune checkpoints. Conclusion Osteomimicry had prognostic value in various cancer types and the underlying mechanism might correlate to the trapping and dysfunction of DCs in the tumor microenvironment, revealing the potential of osteomimicry as a target of immunotherapy.


INTRODUCTION
Osteomimicry referred to the acquisition of genotypic and phenotypic properties of bone cells, predominantly osteoblast, in cancer cells (Rucci and Teti, 2018). It was reported that osteomimicry influenced all aspects of cancer development and metastasis including cell proliferation, adhesion, migration, invasion, and could promote cancer cell survival and angiogenesis (Kruger et al., 2014).
Specially, it was considered that osteomimicry was closely related to cancer metastasis to skeleton. (Rucci and Teti, 2018;Shupp et al., 2018) Once reaching the bone, cancer cells with osteomimicry feature could produce factors that directly [i.e., RANKL (receptor activator of nuclear factor kappa B ligand),  or indirectly (PTHrP, parathyroid hormone-related peptide) promote osteoclastogenesis and bone resorption (Rucci and Teti, 2018).
However, the phenomenon of osteomimicry was mainly discussed in prostate cancer (Koeneman et al., 1999) and breast cancer (Croset et al., 2018;Wang et al., 2019). It was scarcely reported in other cancer types and whether it had prognostic value remained unknown. Thus, we conducted a pan-cancer transcriptomic analysis of osteomimicry across a broad spectrum of solid tumors to define its impact on cancer prognosis using large-scale RNA-sequencing (RNA-seq) data of The Cancer Genome Atlas (TCGA) tumor samples. Since osteonectin (ON, encoded by SPARC), osteopontin (OPN, encoded by SPP1) and osteocalcin (OCN, encoded by BGLAP) were typical bone matrix proteins expressed in the advanced stage of osteoblast differentiation, Teti, 2010, 2018) which indicated extracellular matrix mineralization (Shupp et al., 2018), they were selected as the markers of osteomimicry as reported in other literatures (Koeneman et al., 1999;Rucci and Teti, 2010). Besides, previous studies showed that SPP1 also mediated inflammatory responses by functioning as a chemoattractant for immune cells (Lamort et al., 2019), and SPARC was able to suppress the migration of dendritic cells (DCs) after antigen stimulation (Sangaletti et al., 2005). Thus, the relationship between osteomimicry and immune cell infiltration in the tumor microenvironment was also investigated in this study.

Dataset and Tumor Types
The dataset used consisted of RNA-seq data from TCGA tumor samples (data accessed in Jan 2020). All samples were assayed by RNA-seq, as described by the TCGA research Network. Gene expression values were represented as RNA-Seq by Expectation Maximization (RSEM) data normalized within each sample to the upper quartile of total reads. Acute myeloid leukemia (LAML, n = 173) was excluded because it was hematologic tumor. Lower grade glioma (LGG, n = 530) and pheochromocytoma and paraganglioma (PCPG, n = 187) were excluded because they were not classically malignant tumors. Sarcoma (SARC, n = 265) was excluded because it was bone-derived tumor. Adrenocortical cancer (ACC, n = 79), bile duct cancer (CHOL, n = 45), large B-cell lymphoma (DLBC, n = 48), mesothelioma (MESO, n = 87), ocular melanomas (UVM, n = 80), uterine carcinosarcoma (UCS, n = 57) were excluded for the small sample size (n < 90).

Measurement of Osteomimicry
Osteomimicry was indicated by the gene expression of SPARC, SPP1, and BGLAP, using log 2-transformed values in RSEM. To further explore the effect of osteomimicry on patient's prognosis and tumor microenvironment, patients were divided into high versus low subgroups according to median value of the expression of the osteomimetic markers in the corresponding cancer type.

ESTIMATE (Estimation of STromal and Immune cells in
MAlignant Tumor tissues using Expression data) (Yoshihara et al., 2013) was calculated to predict tumor purity.
The correlation between osteomimicry and immune cell infiltration in the tumor microenvironment was investigated by two independent tools, TIMER (tumor immune estimation resource) 1 (Li et al., 2017) and ImmuCellAI 2 (Miao et al., 2020).

Gene Set Enrichment Analysis
To understand the differences in biological functions and pathways between different subgroups of osteomimicry, gene set enrichment analysis (GSEA 3 , accessed at January, 2020) was performed. We employed the molecular signatures Database (MSigDB) H (hallmark gene sets) collection of chemical and genetic perturbations (n = 20250 gene sets). Calculations were repeated 1000 times for each analysis according to the default weighted enrichment statistical method. GSEA results were shown using normalized enrichment scores (NES), accounting for the size and degree to which a gene set is overrepresented at the top or bottom of the ranked list of genes (nominal p-value < 0.05 and FDR ≤ 0.25).

Statistical Analysis
Associations between subgroups and categorical variables (e.g., sex and disease stage) were analyzed using the chi-square test (Fisher's exact test or Pearson's chi-square test where appropriate), and the Mann-Whitney U test for continuous variables (e.g., age). Correlations between gene expression were evaluated using the spearman correlation test. The correlation was considered weak if the spearman coefficient was less than 0.2, moderate if less than 0.4, relatively strong if less than 0.6, strong if less than 0.8 and very strong if not less than 0.8. The prognostic significance was estimated using Kaplan-Meier survival curves. Cox proportional hazards model was used to calculate the hazard ration (HR) and corresponding 95% confidence interval (CI), incorporating age, sex and disease stage for adjustment. All statistical analyses were performed with R version 3.6.1. Statistical significance was set at p < 0.05 (two-sided).
Overall, the median log 2 -transformed expression value of SPARC, SPP1 and BGLAP were 14.33 (interquartile range, [IQR], 13.17 to 15.38), 11.33 (IQR, 9.02 to 13.59) and 5.08 (IQR, 4.37 to 5.76), respectively. We then divided the samples into two subgroups (high vs. low) by the median value in each cancer type. The demographic and clinical features of the TCGA patients were summarized in Table 2. Generally, both age and gender distribution were similar between subgroups. Patients with tumor stage I/II were allocated into the early stage group and those with tumor stage III/IV were allocated into the advanced stage group.
Osteomimicry was significantly associated with advanced tumor stage in BLCA, BRCA, CESC, COAD, HNSC, KIRC, and STAD, and was marginally significantly associated with advanced tumor stage in ESCA, LIHC, SKCM, and UCEC. Multivariable analysis was performed to explore whether osteomimicry was a prognostic factor for survival outcomes, incorporating clinically relevant covariates (including tumor stage) for adjustment (Figure 3). The results of multivariable modeling largely supported the findings in the univariable analysis. Higher expression of SPARC was an unfavorable prognosticator for patients with COAD (HR = 1.       (Figure 3).

High Expression of All the Three Genes Indicated Unfavorable Prognosis
We further compared the survival outcome of patients with high expression of all the three genes with those with low expression (Figure 4). High expression subgroups had significant unfavorable survival in COAD, HNSC, KIRC, LIHC, and had marginally significantly unfavorable survival in GBM, KIRP, LUAD, LUSC, READ, and UCEC.
Similarly, multivariable analysis was also performed to explore whether patients with high expression of all the three genes had a different survival outcome. We found that patients with high expression of all the three genes had significantly unfavorable prognosis in COAD (HR As expected, when all the 19 cancer types (PRAD was excluded because there was only one death) were pooled, we found that patients with high expression of all the three genes had significantly unfavorable prognosis (HR = 1.61, 95% CI = [1.31, 1.93], p < 0.001) compared with those of low expression (Figure 5).

Correlation Between Osteomimicry and Immune Cell Infiltration
It was reported that the osteomimicry may contribute to an immunosuppressive microenvironment in cancer (Reinstein et al., 2017). Therefore, the patterns of immune cell infiltration were evaluated in the cancer types in which osteomimicry indicated unfavorable survival.
First of all, we evaluated the correlation between osteomimetic markers and tumor purity. The result showed that SPARC and SPP1 expression were positively correlated to ESTIMATE score, but BGLAP expression was weakly and negatively correlated to ESTIMATE score (Figure 6).
Second, we investigated the correlation of osteomimicry with tumor-infiltrating immune cells via two independent methods. Results from TIMER demonstrated that SPARC and FIGURE 5 | Cox proportional hazards analyses of all the three genes across 20 cancer types. Note: age (young vs. old, the median age as the threshold), sex (male vs. female) and disease stage (stage III-IV vs. stage I-II) were included in the multivariable model for adjustment. p < 0.05 represents significant difference in survival outcomes.
SPP1 expression were positively associated with infiltration of DCs and neutrophil cells in all the study models (Figure 7). However, no prominent correlation was found between BGLAP expression and immune cell infiltration (R 2 < 0.40 for all cancer types and all the six immune cell types) (Figure 7). While using ImmuCellAI, we compared 24 infiltrative immune cell types between subgroups, and found that the high expression subgroups had significant more DCs infiltration in the majority of study models (Figure 8).

Correlation Analysis Between Osteomimicry and Immune Checkpoint Expression
We analyzed the correlation between osteomimicry and immune checkpoint expression in the eight cancer types in which osteomimicry harbored prognostic value. All the three markers were found to be correlated to the expression of a series of immune checkpoints. Particularly, SPARC and SPP1 expression were positively correlated to DC-related markers CD86, NRP1, and inhibitory checkpoint markers including LAIR1, HAVCR2, and PDCDLG2 in most cancer types. SPARC expression was positively and significantly correlated to CD276 (B7-H3) in all the eight cancer types. BGLAP expression was positively and significantly correlated to TNFRSF25 in all these cancer types except CESC (Figure 9).

Gene Set Enrichment Analyses of Osteomimetic Expression Subgroups
GSEA was performed to explore the potential biological pathways of osteomimicry. We chose COAD (SPARC), LUSC (SPARC), CESC (SPP1), GBM (SPP1), LIHC (SPP1), LUAD (SPP1), HNSC FIGURE 6 | Correlation between osteomimicry and ESTIMATE score in the cancer types in which osteomimicry indicated unfavorable survival. SPARC (A) and SPP1 (B) expression were positively correlated to ESTIMATE score, but BGLAP (C) expression was weakly and negatively correlated to ESTIMATE score. p < 0.05 represents that difference is statistically significant.
(BGLAP), and LIHC (BGLAP) as the study models, in which the osteomimicry influenced both prognosis and immunological features in these tumors. The common enriched gene sets in the high expression subgroups included unfolded protein response and glycolysis. Under the conditions of FDR < 25% and p-value < 5%, for SPARC, there were 35 and 27 enriched gene sets in the high expression subgroup of COAD and LUSC, respectively, for SPP1, the gene set P53 pathway was enriched in the high expression subgroup of CESC, GBM, LIHC, for BGLAP, the gene set of DNA repair was enriched in the high expression subgroup of HNSC and LIHC.

DISCUSSION
Here, we presented several key aspects of the tumor osteomimicry, identified from RNA-seq data across large-scale TCGA solid tumor samples. Although the concept of osteomimicry was scarcely reported in cancer types other than prostate cancer and breast cancer, we observed the expression of bone markers in all of the twenty cancer types included in our study, indicating osteomimicry might be a common phenomenon in tumor pathophysiology. To our knowledge, this was the first comprehensive transcriptomic investigation of the osteomimetic expression across a large spectrum of solid tumors.
We divided the samples into subgroups based on the median value of each osteomimetic marker and investigated their roles in predicting survival outcomes. Generally, subgroups with high expression of osteomimetic markers had more advanced tumor stage in most of the cancer types. Accordingly, survival analysis showed that the osteomimicry harbored prognostic value in many cancer types, including CESC (SPP1), COAD (SPARC, SPP1, BGLAP, SPARC + SPP1 + BGLAP), GBM (SPP1), HNSC FIGURE 7 | Correlation between osteomimicry and immune cell infiltration in the cancer types in which osteomimicry indicated unfavorable survival. Results from TIMER demonstrated that SPARC (A) and SPP1 (B) expression were positively associated with infiltration of DCs and neutrophil cells in all study models. BGLAP (C) expression showed no prominent correlation with immune cell infiltration. p < 0.05 represents that difference is statistically significant.
There were some reports in regard to the relation between SPARC/SPP1/BGLAP and the prognosis of solid tumors. Ma et al. (2019) performed a systematic review recently and concluded that higher SPARC was associated with poor prognosis in gastrointestinal tumors (gastric cancer, esophageal squamous cell carcinoma, colorectal cancer, and biliary tract cancer), pancreatic cancer and respiratory tract tumors (nonsmall cell lung cancer and nasopharyngeal carcinoma). Higher SPP1 expression had been reported to be correlated to poor prognosis in liver cancer, urothelial carcinoma, breast cancer, colorectal cancer (Zhao et al., 2018). Lewenhaupt et al. found serum BGLAP level was associated with poor prognosis in prostate cancer Lewenhaupt et al., 1988Lewenhaupt et al., , 1990. However, in previous studies, SPARC/SPP1/BGLAP was deemed independent functional molecule rather than marker of osteomimicry, and was investigated separately in different cancer types and no one had combined these molecules and explored their prognostic value in pan-cancer analyses. In our study, for the first time to our knowledge, we revealed that high expression of these osteomimetic markers was associated with unfavorable survival across a large spectrum of solid tumors. Furthermore, to better distinguish osteomimicry, patients were divided into subgroups with high versus low expression of all these three genes. Similarly, subgroups of high expression of all the three genes had worse survival, and showed higher HR (1.61) than that calculated with only one marker (1.10 for SPARC, 1.25 for SPP1, 1.13 for BGLAP). The above observations further revealed the importance of osteomimicry in the prediction of prognosis of various cancer types.
Tumor microenvironment has long been shown playing an important role in cancer progression, response to clinical intervention including immunotherapy, and clinical outcome of cancer patients (Schelker et al., 2017;Binnewies et al., 2018;Deng et al., 2020). However, whether osteomimicry had any impact on the tumor microenvironment remained obscure. The correlation between osteomimicry and tumor purity was then investigated, and result showed that SPARC and SPP1 expression had a moderate-to-strong and positive correlation with tumor purity, while BGLAP expression had a poor-to-moderate and negative correlation with tumor purity. This was consistent with the FIGURE 8 | Violin plots of comparison of immune cell infiltration between high and low subgroups in the cancer types in which osteomimicry indicated unfavorable survival. Results from ImmuCellAI indicated that the high expression subgroups of SPARC (A), SPP1 (B), and BGLAP (C) had significant more DCs infiltration in the majority of study models. p < 0.05 represents that difference is statistically significant.
process of osteoblast differentiation, which could be divided into three stages: proliferation, extracellular matrix maturation, and extracellular matrix mineralization. SPARC, SPP1, and BGLAP genes were sequentially expressed in the three stages of osteoblast differentiation (Shupp et al., 2018). High SPARC expression was observed in newly mineralizing bone but not in osteoid, while increase in BGLAP synthesis was observed during the transition from osteoid to mature mineralized matrix (Nakase et al., 1994;Licini et al., 2019). Therefore, higher BGLAP expression might indicate complete remodeling of the tumor microenvironment and cause a relatively low tumor purity.
Further, we investigated the immune cell infiltration in the tumor microenvironment using two different databases. Interestingly, the results from both TIMER and ImmuCellAI demonstrated that patients with a high osteomimetic expression harbored a high intratumoral DCs infiltration. DCs had been known as the antigen presenting cells and were able to prime, activate and direct the T cells to target tumor cells (Saxena and Bhardwaj, 2018). Previous reports had demonstrated that SPARC and SPP1 were able to suppress the migration of DCs to lymph nodes, which would result in the accumulation of DCs at primary tumor site (Sangaletti et al., 2005;Begum et al., 2007). Recent studies also showed that local factors within tumor microenvironment, including hypoxia, adenosine and lactate accumulations and decreased pH, could lead to dysfunction of DCs which were trapped at primary tumor site (Veglia and Gabrilovich, 2017). IL-10 produced by tumor-associatedmacrophage could also inhibit IL-12 production by DCs and alter their ability to mount antigen specific T cell responses (Ruffell et al., 2014). The higher expression of inhibitory molecules (such as PD-L1) (Salmon et al., 2016) and accumulating lipid (Herber et al., 2010) in the tumor microenvironment could also lead to the dysfunction of DCs. It was reported that DCs with abnormal function in the tumor microenvironment can produce tumor promoting IL-6 and immunosuppressive galectin-1, which can impair the local anti-tumor immunity (Tesone et al., 2016). Therefore, it was reasonable to speculate that osteomimicry might trap the DCs at primary tumor site, leading to dysfunction and even immunosuppression of DCs and thus a poor prognosis of patients.
Certainly, the bone matrix protein SPARC/SPP1/BGLAP were functional molecules and played an important role in various malignancies according to previous literatures. Besides regulating bone mineralization, SPP1 has also been shown playing a role in tumor cell proliferation, adhesion, invasion and metastasis (Zhao et al., 2018;Pang et al., 2019). It was reported that SPP1 could regulate MMP-2 and MMP-9 expression and promote extracellular matrix degradation via activating αvβ-NF-κB pathway, and thus accelerate the growth, migration and invasion of cancer cells (Fong et al., 2009;Liu et al., 2010;Qin et al., 2018). Some studies also showed that SPP1 contributed to tumor cell migration, invasion and survival through PI3K-Akt (Fong et al., 2009;Zhang et al., 2014;Wei et al., 2018) and STAT3 (Behera et al., 2010) signaling, and could also promote tumor growth and metastasis by inducing angiogenesis (Dai et al., 2009;Wang et al., 2011). SPARC was a matricellular protein modulating cell-matrix interactions and was found up-regulated in tumor stroma and associated with poor prognosis in many cancer types, including colorectal cancer (Kurtul et al., 2017;Drev et al., 2019), pancreatic cancer (Vaz et al., 2015), cervical carcinoma (Shi et al., 2016), and non-small cell lung cancer (Koukourakis et al., 2003). It's worth noting that the influence of SPARC on tumor progression depended on the initiating cell type, tumor stage and context of the microenvironment (Tai and Tang, 2008;Arnold and Brekken, 2009;Rossi et al., 2016). Notably, Munasinghe et al. reported that depletion of fibronectin FIGURE 9 | Correlation between osteomimetic markers (A) SPARC, (B) SPP1 and (C) BGLAP and immune checkpoints in the cancer types in which osteomimicry indicated unfavorable survival. SPARC (A) and SPP1 (B) expression were positively correlated to DC-related markers CD86, NRP1 and inhibitory checkpoint markers including LAIR1, HAVCR2 and PDCDLG2 in most cancer types. SPARC (A) expression was positively and significantly correlated to CD276 (B7-H3) in all the eight cancer types. BGLAP (C) expression was positively and significantly correlated to TNFRSF25 in all these cancer types except CESC. *indicated p < 0.05. **indicated p < 0.01. ***indicated p < 0.001. p < 0.05 represents that difference is statistically significant.
switched the function of SPARC from promoting cancer cell proliferation to growth inhibition and induction of apoptosis (Munasinghe et al., 2020). Although BGLAP was known crucial during bone remodeling, the function remained largely unknown and it arose great interest recently as a bone-derived hormone in regulation of energy metabolism (Li et al., 2016;Zoch et al., 2016;Mera et al., 2018;Moser and van der Eerden, 2018).
Among the solid tumors investigated in this study, osteomimicry was found unable to predict prognosis in some cancer types. Although osteomimicry was used to explain the predilection of skeleton metastasis in prostate cancer and breast cancer, it showed no correlation with advanced tumor stage or poor survival in PRAD in this study. BGLAP was correlated to advanced tumor stage in BRCA, but showed no prognostic value. The failure to detect association might due to lack of adequate follow-up time, limited event rates, biased population distribution or medical intervention.
The main limitation of this study is lack of ability of osteomimicry to predict the response to immunotherapy, and it requires further validation in cancer patients treated with immune checkpoint inhibitors. Nevertheless, our findings are important and provide new insights into the prognostic and immunological features of osteomimicry in various solid tumors. Future molecular studies are warranted to shed light on how osteomimicry affects different cancers and survival outcomes.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
Ethical review and approval was not required for the study on Human participants in accordance with the Local Legislation and Institutional Requirements. Written informed consent from the participants' legal guardian/next of kin was not required to participate in this study in accordance with the National Legislation and the Institutional Requirements.

AUTHOR CONTRIBUTIONS
CY designed this study and contributed to data collection, data analysis, data interpretation, and writing the manuscript. LS helped to design this study and contributed to data collection, data analysis and data interpretation, and revising the manuscript. HP helped to design this study and contributed to data interpretation and revising the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
The results published here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga.