Identification of candidate genes for myeloma-induced osteocyte death based on microarray data
© The Author(s). 2016
Received: 26 May 2016
Accepted: 30 June 2016
Published: 12 July 2016
The study was aimed to investigate the molecular mechanisms of osteocyte death in multiple myeloma (MM) patients.
GSE27372 was downloaded from Gene Expression Omnibus, including three HOB-01 (osteocyte cell line) control samples and three HOB-01 samples co-cultured with JJN3 (human MM cell line). After the differentially expressed genes (DEGs) were identified by Student’s t test method, enrichment analyses were performed for them using DAVID software. Using TRANSFAC, TSGene, and tumor-associated gene (TAG) databases, functional annotation was conducted for the DEGs. Additionally, protein-protein interaction (PPI) network and sub-network analyses were performed using STRING database and Cytoscape software.
Total 393 DEGs were identified, including 22 transcription factors (e.g., KLF4 and IRF8) and 37 TAGs. Enrichment analysis suggested that EGF, S1PR1, and NPY1R were enriched in the function of circulatory system development. EGF (degree = 31) and EGR1 (degree = 19) had high degrees and interactions in the PPI network. In the sub-network, S1PR1, C3AR1, and NPY1R could interact with each other.
These DEGs might participate in the osteocyte apoptosis induced by myeloma cells. These findings might provide a theoretical basis for a better understanding of the osteolysis in MM patients.
KeywordsMultiple myeloma Osteocyte Differentially expressed genes Enrichment analysis Protein-protein interaction network
Multiple myeloma (MM) originates in neoplastic plasma cell disorder, and it is characterized by the clonal proliferation of malignant plasma cells in the bone marrow. As the second most general hematological cancer, the incidence of MM worldwide is about 1.5/100,000 new case . It is also found that the incidence of MM is higher in men than in women, as well as in black than white in the USA . MM is always associated with organ dysfunction, especially osteolysis . Although the survival rate of MM patients is increased due to the advanced medicines , MM continues to be considered as an incurable disease, and further study is required to fully understand the potential molecular mechanisms of osteolysis in MM patients.
To adopt a volume appropriate for the local environment, the bone continuously remodels to keep the balance between bone formation and resorption mediated by osteocytes, osteoblasts, and osteoclasts . Osteocyte represents the most abundant cell type in the bone and is actively involved in bone turnover. Dickkopf 1(DKK1) levels increase in the peripheral blood and bone marrow plasma of MM patients and relate to the development of osteolytic lesions [6, 7]. Tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) produced by myeloma cells is positively correlated with osteolytic markers (such as urinary deoxypyridinoline and serum calcium), indicating that TRAIL may function in osteolysis of MM patients . Via promoting the expression of receptor activator of NF-kB (RANK) in osteoclast precursors, c-Akt (AKT) plays a role in the osteoclast formation and bone osteolysis induced by MM . Giuliani et al. found that MM cells promoted osteocyte death and altered the transcriptional profile in osteocyte . However, they did not further perform comprehensive bioinformatics analysis to investigate the internal causes. In spite of the above researches, the molecular mechanism of the increased osteocyte death in MM patients is still unclear.
In our study, in order to investigate the molecular mechanism of osteocyte death in MM patients, we reanalyzed the gene expression profile in Giuliani et al. study and identified the differentially expressed genes (DEGs) between normal osteocytes and osteocytes affected by MM cells. The enriched functions and pathways of DEGs were further identified. In addition, we constructed the protein-protein interaction (PPI) network of DEGs. Then, the most significant sub-network was screened out.
Affymetrix microarray data
Gene transcriptional profile of GSE27372  was downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), and GSE27372 was based on the platform of GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 array. There were six specimens in this dataset, including three HOB-01 (osteocyte cell line) control samples and three HOB-01 samples co-cultured with JJN3 (human myeloma cell line).
AFFY package (version 1.28.0, http://www.bioconductor.org/packages/release/bioc/html/affy.html) of Bioconductor  was used to pre-process the Affymetrix microarray data. All raw data were scaled according to the robust multi-array average (RMA) method  with default settings. After background correction, quantile normalization, and probe summarization, the gene expression matrix was obtained. The Student’s t test method  was adopted to analyze the expression differences between control and co-culture group. For each significant DEG, both P value <0.05 and |log2 fold change (FC)| >0.58 need to be met.
Function and pathway enrichment analyses of DEGs
Functional annotation for DEGs was performed by Gene Ontology . Kyoto Encyclopedia of Genes and Genomes (KEGG http://www.genome.jp/kegg/pathway.html) pathway enrichment analysis was used to identify main functional and metabolic pathways involving DEGs. We used P value <0.01 as the cut-off criterion for the enrichment analysis which was conducted by the Database for Annotation, Visualization and Integrated Discovery (DAVID; version 2.1b, http://david.abcc.ncifcrf.gov/) online software .
Transcription factors (TFs) and tumor-associated genes (TAGs) in DEGs
DEGs with the function of transcriptional regulation, namely, differentially expressed transcription factors were selected based on the TRANSFAC database (http://www.gene-regulation.com/pub/databases.html) . According to the TSGene database (http://bioinfo.mc.vanderbilt.edu/TSGene/search.cgi)  and TAG database (http://www.binfo.ncku.edu.tw/TAG/) , we extracted the TAGs from DEGs, including tumor suppressor genes and oncogenes.
Construction of PPI network and sub-network
The interaction pairs of DEGs were analyzed via online tool Search Tool for the Retrieval of Interacting Genes (STRING; version 9.0, http://string-db.org) , and interaction data were downloaded on June 27, 2014. Only the gene pairs which were recorded in database, and experimental validated, text-mined, or co-expressed were used and combined score >0.4 was set as the criterion of PPI. Then, the PPI network was constructed using Cytoscape (version 2.8, http://cytoscape.org) .
The sub-network analysis of PPI network was performed using the ClusterONE plug of Cytoscape . The pathway enrichment analysis of the genes in the most significant sub-network was performed by using the DAVID online software .
MM-induced DEGs in osteocytes
After analyzing the microarray data of control and co-culture groups, a total of 393 DEGs were screened out, including 167 down-regulated genes and 226 up-regulated genes in co-culture group.
Functions and pathways enriched by DEGs
Function enrichment analysis of DEGs
Cardiovascular system development
Circulatory system development
Anatomical structure morphogenesis
Extracellular matrix part
Extracellular matrix structural constituent
Structural molecule activity
Regulation of response to stimulus
Regulation of signaling
Regulation of cell communication
Extracellular region part
Platelet dense tubular network membrane
Sequence-specific DNA binding
RNA polymerase II core promoter proximal region sequence-specific DNA binding transcription factor activity
Sequence-specific DNA binding transcription factor activity
Pathway enrichment analysis of DEGs
Protein digestion and absorption
Regulation of actin cytoskeleton
Hypertrophic cardiomyopathy (HCM)
JAK-STAT signaling pathway
Nicotinate and nicotinamide metabolism
Calcium signaling pathway
VEGF signaling pathway
Phosphatidylinositol signaling system
Cytokine-cytokine receptor interaction
Differentially expressed TFs and TAGs
Differentially expressed TFs and TAGs
Tumor suppressor genes
TGFB1I1, NKX2-2, KLF5, ID2
RUNX1T1, PDGFRB, CD24
TPM1, SPTBN1,PRICKLE1, KLF5, IGFBP3, FAT4,EFNA1, DNAJB4, DAB2, CADM1
TGFB2, LRRC17, BAMBI
SOX9, SOX7, NR3C2, MSX1, MAFF, KLF4, JUNB, IRF8, HOXB6, HMGA2, HES1, OXF1, FOXA2, ETV5, ETV4, ELK3, EGR3, EGR1
VEGFA, SPHK1, PDGFRA, JUNB, HMGA2, FOSL1, FGF5, ETV1, CCND1
SPRY2, SOX7, IRF8, IL24, HOPX, FOXA2, EPHB2, ENC1, EGR1, DUSP6
Furthermore, 16 down-regulated TAGs and 21 up-regulated TAGs were screened out (Table 3). Among the down-regulated TAGs, there were 10 tumor suppressor genes and 3 oncogenes. Meanwhile, among the up-regulated TAGs, there were 10 tumor suppressor genes (e.g. IRF8) and 9 oncogenes.
Construction of PPI network and sub-network
Pathways enriched by the genes in sub-network
Neuroactive ligand-receptor interaction
C3AR1, S1PR1, NPY1R
To gain insight into the molecular mechanisms of the myeloma-induced osteocyte death, gene expression profiles in osteocytes co-cultured with or without myeloma cells were systematically analyzed here. A total of 226 up-regulated genes and 167 down-regulated genes were identified. Then, PPI network was constructed, and sub-network was identified.
TFs like KLF4 and IRF8 were up-regulated in osteocytes co-cultured with myeloma cells in comparison with osteocytes cultured alone. Previous study reported that KLF4 expression is essential for blocking cell cycle and increasing the resistance of MM cells to alkylating agents . As one of the Yamanaka reprogramming factors and TFs, KLF4 can promote the expression of autophagy-related genes . In addition, silencing of the interferon consensus sequence-binding protein (ICSBP/IRF8) gene may be induced by DNA methylation or other mechanisms and correlates with the malignant phenotype of MM . In our study, KLF4 and IRF8 were up-regulated TFs, and IRF8 was a tumor suppressor gene. These suggested that KLF4 and IRF8 might play a role in osteolysis in MM patients.
In this study, EGF and EGR1 had high connectivity degrees in PPI network. A heparin-binding factor in EGF family, amphiregulin (AREG) is overexpressed in primary myeloma cells and can promote growth of bone marrow stromal cell . EGR-1 induces apoptosis in MM via interacting with JUN, and decreased JUN/EGR-1 can enhance resistance of MM cells to bortezomib . Functional enrichment indicated that EGF was involved in circulatory system development. It has been reported that the substances of the circulatory system can induce the apoptosis of tumor cells . EGF could interact with EGR1 in the PPI network, indicating that EGF and EGR1 might be involved in osteocytes apoptosis induced by MM cells through interacting with each other.
Enrichment analysis suggested that S1PR1, C3AR1, and NPY1R in sub-network were enriched in the pathway of neuroactive ligand-receptor interaction, meanwhile S1PR1 and NPY1R were involved in the function of circulatory system development. S1PR1 is a G-protein-coupled receptor which binds to the bioactive signaling molecule sphingosine 1-phosphate (S1P) . It is reported that S1PR1 and S1PR2 regulate osteoclast precursor migration between the bone marrow cavities and the circulation . It is shown that component 3 (C3A) and its receptor C3AR play a role in osteoclast formation [30, 31], implying a potential role of C3AR in osteocyte death. Thus, S1PR1 and C3AR1 may promote the death of osteocytes, and this is consistent with the finding of Giuliani et al. . NPY1R is a receptor of neuropeptide Y (NPY), which is a neurotransmitter. There is evidence that NPY regulates bone homeostasis via actions in peripheral tissues . NPY1R is the only Y receptor which is robustly expressed in bone marrow stromal cells  and osteoblasts , implying a direct role of NPY1R in bone remodeling. In addition, the absence of peripheral Y1 receptor will lead to pronounced anabolic effects on the bone . Here, NPY1R was significantly down-regulated in osteocytes co-cultured with myeloma cells, inhibiting bone remodeling. In the sub-network, S1PR1, C3AR1, and NPY1R had interactions with each other, suggesting that S1PR1, C3AR1, and NPY1R might function in osteocytes apoptosis induced by MM cells via interactions.
In conclusion, a total of 393 DEGs were identified between osteocytes co-cultured with and without myeloma cells, and KLF4, IRF8, EGF, EGR1, S1PR1, C3AR1, and NPY1R might be involved in osteocyte cell apoptosis induced by MM cells. Here, we studied the potential molecular mechanism of the increased osteocyte death, which may provide a theoretical basis for understanding and treatment of osteolysis in MM patients. However, the present study analyzed gene expression data generated from not primary cells isolated from myeloma patients and healthy controls but cell lines, and the sample size was small. Thus, confirmation of results by quantitative real-time polymerase chain reaction (qRT-PCR) in cell lines followed by the investigation of the roles and functions of candidate genes in cells isolated from MM patients and healthy controls is still needed.
DEGs, differentially expressed genes; DKK1, Dickkopf 1; ECM, extracellular matrix; EGF, epidermal growth factor; KEGG, Kyoto Encyclopedia of Genes and Genomes; MM, multiple myeloma; NPY1R, neuropeptide Y1 receptor; PPI, protein-protein interaction; RMA, robust multi-array average; S1P, sphingosine 1-phosphate; S1PR1, sphingosine-1-phosphate receptor 1; STRING, Search Tool for the Retrieval of Interacting Genes; TAG, tumor-associated gene; TAGs, tumor-associated genes; TFs, transcription factors; TRAIL, tumor necrosis factor-related apoptosis-inducing ligand
This work was supported by special fund for medical service of Jilin finance department (No.: SCZSY201507).
Availability of data and materials
The dataset “Gene transcriptional profile of GSE27372” supporting the conclusions of this article can be downloaded from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). GSE27372 was based on the platform of GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array.
HT was responsible to the conception and design of the research, acquisition of data. obtaining funding, drafting the manuscript and revision of the manuscript for important intellectual content.
The author declares that they have no competing interests.
Ethics approval and consent to participate
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Campa D, Martino A, Varkonyi J, Lesueur F, Jamroziak K, Landi S, Jurczyszyn A, Marques H, Andersen V, Jurado M. Risk of multiple myeloma is associated with polymorphisms within telomerase genes and telomere length. Int J Cancer. 2015;136:E351–E8.View ArticlePubMedGoogle Scholar
- Howlader N, Noone A, Krapcho M, Neyman N, Aminou R, Altekruse S, Kosary C, Ruhl J, Tatalovich Z, Cho H. SEER cancer statistics review, 1975–2009 (vintage 2009 populations). Bethesda: National Cancer Institute; 2012.Google Scholar
- Rajkumar SV. Multiple myeloma: 2011 update on diagnosis, risk‐stratification, and management. Am J Hematol. 2011;86:57–65.View ArticlePubMedGoogle Scholar
- Bergsagel PL, Mateos M-V, Gutierrez NC, Rajkumar SV, San Miguel JF. Improving overall survival and overcoming adverse prognosis in the treatment of cytogenetically high-risk multiple myeloma. Blood. 2013;121:884–92.View ArticlePubMedPubMed CentralGoogle Scholar
- Crockett JC, Rogers MJ, Coxon FP, Hocking LJ, Helfrich MH. Bone remodelling at a glance. J Cell Sci. 2011;124:991–8.View ArticlePubMedGoogle Scholar
- Tian E, Zhan F, Walker R, Rasmussen E, Ma Y, Barlogie B, Shaughnessy Jr JD. The role of the Wnt-signaling antagonist DKK1 in the development of osteolytic lesions in multiple myeloma. N England J Med. 2003;349:2483–94.View ArticleGoogle Scholar
- Kocemba KA, Groen R, Van Andel H, Kersten MJ, Mahtouk K, Spaargaren M, Pals ST. Transcriptional silencing of the Wnt-antagonist DKK1 by promoter methylation is associated with enhanced Wnt signaling in advanced multiple myeloma. PLoS One. 2012;7:e30359.View ArticlePubMedPubMed CentralGoogle Scholar
- Kawano Y, Ueno S, Abe M, Kikukawa Y, Yuki H, Iyama K, Okuno Y, Mitsuya H, Hata H. TRAIL produced from multiple myeloma cells is associated with osteolytic markers. Oncol Rep. 2012;27:39–44.PubMedGoogle Scholar
- Cao H, Zhu K, Qiu L, Li S, Niu H, Hao M, Yang S, Zhao Z, Lai Y, Anderson JL. Critical role of AKT protein in myeloma-induced osteoclast formation and osteolysis. J Biol Chem. 2013;288:30399–410.View ArticlePubMedPubMed CentralGoogle Scholar
- Giuliani N, Ferretti M, Bolzoni M, Storti P, Lazzaretti M, Dalla Palma B, Bonomini S, Martella E, Agnelli L, Neri A. Increased osteocyte death in multiple myeloma patients: role in myeloma-induced osteoclast formation. Leukemia. 2012;26:1391–401.View ArticlePubMedGoogle Scholar
- Gautier L, Cope L, Bolstad BM, Irizarry RA. affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004;20:307–15.View ArticlePubMedGoogle Scholar
- Irizarry RA, Hobbs B, Collin F, Beazer‐Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4:249–64.View ArticlePubMedGoogle Scholar
- Haynes W. Student’s t-test. Encyclopedia of systems biology. New York; Springer; 2013; p. 2023–5.Google Scholar
- Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protocols. 2008;4:44–57.View ArticleGoogle Scholar
- Matys V, Fricke E, Geffers R, Gößling E, Haubrock M, Hehl R, Hornischer K, Karas D, Kel AE, Kel-Margoulis OV. TRANSFAC®: transcriptional regulation, from patterns to profiles. Nucleic Acids Res. 2003;31:374–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhao M, Sun J, Zhao Z. TSGene: a web resource for tumor suppressor genes. Nucleic Acids Res. 2013;41:D970–D6.View ArticlePubMedGoogle Scholar
- Chen J-S, Hung W-S, Chan H-H, Tsai S-J, Sun HS. In silico identification of oncogenic potential of fyn-related kinase in hepatocellular carcinoma. Bioinformatics. 2013;29:420–7.View ArticlePubMedGoogle Scholar
- Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, Lin J, Minguez P, Bork P, von Mering C. STRING v9. 1: protein-protein interaction networks, with increased coverage and integration. Nucleic Acids Res. 2013;41:D808–D15.View ArticlePubMedGoogle Scholar
- Smoot ME, Ono K, Ruscheinski J, Wang P-L, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics. 2011;27:431–2.View ArticlePubMedGoogle Scholar
- Demchak B, Hull T, Reich M, Liefeld T, Smoot M, Ideker T, Mesirov JP. Cytoscape: the network visualization tool for GenomeSpace workflows. F1000Research. 2014;3:151.PubMedPubMed CentralGoogle Scholar
- Schoenhals M, Kassambara A, Veyrune J-L, Moreaux J, Goldschmidt H, Hose D, Klein B. Krüppel-like factor 4 blocks tumor cell proliferation and promotes drug resistance in multiple myeloma. Haematologica. 2013;98:1442–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Riz I, Hawley TS, Hawley RG. KLF4-SQSTM1/p62-associated prosurvival autophagy contributes to carfilzomib resistance in multiple myeloma models. Oncotarget. 2015;6:14814.View ArticlePubMedPubMed CentralGoogle Scholar
- Tshuikina M, Jernberg-Wiklund H, Nilsson K, Oberg F. Epigenetic silencing of the interferon regulatory factor ICSBP/IRF8 in human multiple myeloma. Exp Hematol. 2008;36:1673–81. e1.View ArticlePubMedGoogle Scholar
- Mahtouk K, Hose D, Rème T, De Vos J, Jourdan M, Moreaux J, Fiol G, Raab M, Jourdan E, Grau V. Expression of EGF-family receptors and amphiregulin in multiple myeloma. Amphiregulin is a growth factor for myeloma cells. Oncogene. 2005;24:3512–24.View ArticlePubMedPubMed CentralGoogle Scholar
- Chen L, Wang S, Zhou Y, Wu X, Entin I, Epstein J, Yaccoby S, Xiong W, Barlogie B, Shaughnessy Jr JD. Identification of early growth response protein 1 (EGR-1) as a novel target for JUN-induced apoptosis in multiple myeloma. Blood. 2010;115:61–70.View ArticlePubMedPubMed CentralGoogle Scholar
- Kulcsár G. Apoptosis of tumor cells induced by substances of the circulatory system. Cancer Biother Radiopharm. 1997;12:19–26.View ArticlePubMedGoogle Scholar
- Maceyka M, Harikumar KB, Milstien S, Spiegel S. Sphingosine-1-phosphate signaling and its role in disease. Trends Cell Biol. 2012;22:50–60.View ArticlePubMedGoogle Scholar
- Ishii M, Kikuta J. Sphingosine-1-phosphate signaling controlling osteoclasts and bone homeostasis. Biochimica et Biophysica Acta (BBA)-Molecular and Cell Biology of Lipids. 2013;1831:223–7.View ArticleGoogle Scholar
- Abe M, Hiura K, Wilde J, Shioyasono A, Moriyama K, Hashimoto T, Kido S, Oshima T, Shibata H, Ozaki S. Osteoclasts enhance myeloma cell growth and survival via cell-cell contact: a vicious cycle between bone destruction and myeloma expansion. Blood. 2004;104:2484–91.View ArticlePubMedGoogle Scholar
- Ignatius A, Schoengraf P, Kreja L, Liedert A, Recknagel S, Kandert S, Brenner RE, Schneider M, Lambris JD, Huber‐Lang M. Complement C3a and C5a modulate osteoclast formation and inflammatory response of osteoblasts in synergism with IL‐1β. J Cell Biochemistry. 2011;112:2594–605.View ArticleGoogle Scholar
- Sousa D, Herzog H, Lamghari M. NPY signalling pathway in bone homeostasis: Y1 receptor as a potential drug target. Current Drug Targets. 2009;10:9–19.View ArticlePubMedGoogle Scholar
- Lundberg P, Allison SJ, Lee NJ, Baldock PA, Brouard N, Rost S, Enriquez RF, Sainsbury A, Lamghari M, Simmons P. Greater bone formation of Y2 knockout mice is associated with increased osteoprogenitor numbers and altered Y1 receptor expression. Journal of Biological Chemistry. 2007;282:19082–91.View ArticlePubMedGoogle Scholar
- Teixeira L, Sousa DM, Nunes AF, Sousa MM, Herzog H, Lamghari M. NPY revealed as a critical modulator of osteoblast function in vitro: new insights into the role of Y1 and Y2 receptors. Journal of cellular biochemistry. 2009;107:908–16.View ArticlePubMedGoogle Scholar
- Lee NJ, Nguyen AD, Enriquez RF, Doyle KL, Sainsbury A, Baldock PA, Herzog H. Osteoblast specific Y1 receptor deletion enhances bone mass. Bone. 2011;48:461–7.View ArticlePubMedGoogle Scholar