Edited by: Tao Liu, University of New South Wales, Australia
Reviewed by: Prashant Trikha, Nationwide Children's Hospital, United States; Sabarish Ramachandran, Texas Tech University Health Sciences Center, United States
This article was submitted to Molecular and Cellular Oncology, a section of the journal Frontiers in Oncology
†These authors have contributed equally to this work
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Breast cancer, a complex and heterogeneous disease, is the leading cause of cancer-associated death in the women worldwide (
It is reported that the pathology of bone metastasis is characterized by a vicious cycle where the balance between bone forming cells “osteoblasts” and the bone resorbing cells “osteoclasts” is imbalanced (
In this study, we aimed to find differentially expressed genes (DEGs) in breast cancer bone metastasis by integrated analysis. Then, functional enrichment analysis including Gene Ontology (GO) and the Kyoto Encyclopedia of Genes Genomes (KEGG) was used to investigate the biological function of DEGs. Protein-protein interaction (PPI) and transcriptional factors (TFs)-target genes network construction was performed to further investigate the function of identified DEGs. Quantitative reverse transcription-polymerase chain reaction (QRT-PCR) was applied to validate the expression of candidate DEGs. Finally, receiver operating characteristic (ROC) analyses was applied to analyze diagnosis ability of identified DEGs. Our study may be helpful in understanding the pathogenic mechanism of bone metastasis from breast cancer, and most difference from previous studies is that current study focuses on bone metastasis of breast cancer and aim to find the differential expressed genes and pathways from other site metastasis of breast cancer.
In this study, we searched datasets of breast cancer bone metastasis from the Gene Expression Omnibus (GEO) database (
Three datasets of selected breast cancer with bone metastasis.
GSE14020 | GPL96 [HG-U133A] Affymetrix Human Genome U133A Array | 28 | 8 | 2009 | USA | 19573813; 25670079 |
GPL570[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array | 19 | 10 | ||||
GSE54323 | GPL570[HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array | 15 | 14 | 2015 | Sweden | 25888067 |
GSE46141 | GPL10379 Rosetta/Merck Human RSTA Custom Affymetrix 2.0 microarray [HuRSTA-2a520709] | 86 | 5 | 2013 | Sweden | 24287398 |
Raw expression data in this study was downloaded. Limma and Meta-analysis for MicroArrays (metaMA) packages were used to identify the DEGs, and the inverse normal method was used to combine the
In order to obtain the biological function and signaling pathways of DEGs, the GeneCoDis3 (
It is useful to understand the molecular mechanism of bone metastasis from breast cancer by studying the interactions between proteins. In order to gain insights into the interaction between proteins encoded by DEGs and other proteins, the database of BioGRID (
TFs play crucial roles in regulating gene expression. We downloaded the TFs in the human genome and the motifs of genomic binding sites from the TRANSFAC. Moreover, the 2 KB sequence in the upstream promoter region of DEGs was downloaded from UCSC (
Twenty-five patients with bone metastasis of breast cancer, 25 patients with lung metastasis of breast cancer, and 25 patients with liver metastasis of breast cancer were incorporated in our study. Meanwhile, 25 primary breast cancer patients without metastasis was enrolled, and their adjacent non-cancer breast tissues was considered as normal tissue to be used as control in the present study. Their tumor tissues from metastatic sites were used for RNA isolation. Ethical approval was obtained from the ethics committee of the cancer hospital and informed written consent was obtained from all of subjects.
Total RNA of the tumor tissues from metastatic sites was extracted using the RNA liquid Reagent (TIANGEN) according to the manufacturer's protocols. One microgram RNA was applied to synthesize cDNA by Fast Quant Reverse Transcriptase (TIANGEN). Then real-time PCR was performed in an ABI 7300 Real-time PCR system with SYBR® Green PCR Master Mix (TIANGEN). All reactions were carried out in triplicate. GAPDH was used for the internal reference. Relative gene expression was analyzed by fold change method.
In order to assess the diagnostic value of IBSP, MMP9, TNFAIP6, DHRS3, RIPK4, CSPG4, and CD200 in bone metastasis of breast cancer, receiver operating characteristic (ROC) analyses were performed using pROC package in R language. The area under the curve (AUC) under binomial exact confidence interval was calculated to generate the ROC curve.
Human breast cancer cell line MCF-7 was purchased from ATCC (Manassas, VA, USA) and was cultured in RPMI 1640 medium (Invitrogen, CA, USA) supplemented with 10% fetal bovine serum (Invitrogen, Carlsbad, CA, USA) and incubated at 37°C and 5% CO2. CD200 and TNFAIP6 were selected from up-regulated gene groups and DHRS3, ASS1, and RIPK4 were selected from down-regulated gene groups. Commercial RNA interfering sets targeting CD200 and TNFAIP6 were purchased from Ribobio technology (Shenzhen, Guangdong, China), and negative control of siRNA was also transfected into MCF-7 cells to avoid the influence of transfection itself. Exogenous human full-length DHRS3, ASS1, and RIPK4 plasmid vectors were purchased from Sino biological Inc. (Beijing, China); and empty vector was also transfected into MCF7 cells as negative control. Lipofectamine 2000 (Invitrogen) was used for all transfection studies by following the manufacturer's protocol. Transfection efficiency was examined using qRT-PCR. Briefly, total RNA was extracted from cells using Trizol (Invitrogen, CA, USA) and further purified by RNAeasy mini kit plus DNase I treatment (Qiagen, Germany). The relative mRNA level for each gene was quantified by real-time RT-PCR with SYBR Green (Applied Biosystems, CA, USA), using GAPDH as a control.
The migration ability of cells was measured by wound-healing assay. Briefly, 2 × 104 cells were inoculated in 6-cm tissue culture dishes, cultured overnight, scratch was performed when cell growth fusion reached to 90%. Migration images were captured 24 and 48 h after scratching, the migration rate (width) of cells were calculated.
Corning™ BioCoat™ Matrigel™ Invasion Chamber with Corning™ Matrigel Matrix (Thermo Fisher Scientific, Waltham, MA, USA) was used for cell invasion assay. One hundred microliters of cells with concentration of 5 × 105/mL were added in the upper chamber with (8 μm pore), 600 μL medium containing 10% serum was added in the lower chamber. After cultured for 6 h, medium in the chamber were removed and unmigrated cells were swabbed. Cells was fixed by 4% polyformaldehyde for 10 min, stained with crystal violet for another 10 min. The filter membrane was photographed under the inverted microscope (200×) after sealed with the neutral gum. Cells were counted by Image Pro Plus Version 6 software, 3 wells in each group and 5 vision fields of each well were randomly selected to calculate the average number of cells.
A total of 749 DEGs were identified with the threshold of
Top 10 up- and down-regulated DEGs.
3381 | IBSP | 0.091052 | −1.13987 | 8.66E-15 | 1.06E-10 | Up |
4318 | MMP9 | 2.196633 | 0.711473 | 4.50E-12 | 2.74E-08 | Up |
4322 | MMP13 | 0.044185 | −0.65971 | 3.82E-10 | 1.55E-06 | Up |
4958 | OMD | 0.009911 | −0.95916 | 1.19E-09 | 3.64E-06 | Up |
1513 | CTSK | 1.700623 | 0.648007 | 3.23E-09 | 7.88E-06 | Up |
7130 | TNFAIP6 | 0.425532 | −0.36804 | 1.81E-08 | 3.67E-05 | Up |
23314 | SATB2 | 0.126353 | −0.60566 | 3.32E-08 | 5.78E-05 | Up |
25903 | OLFML2B | 0.748614 | 0.212792 | 4.36E-08 | 6.65E-05 | Up |
10085 | EDIL3 | −0.93194 | −1.45248 | 6.10E-08 | 8.27E-05 | Up |
5118 | PCOLCE | 1.030222 | 0.374014 | 1.25E-07 | 0.000151 | Up |
9249 | DHRS3 | 0.488833 | 0.944744 | 5.24E-07 | 0.000288 | Down |
445 | ASS1 | 0.439835 | 1.305644 | 2.88E-06 | 0.001086 | Down |
54101 | RIPK4 | 0.426273 | 0.68844 | 3.83E-06 | 0.001262 | Down |
7431 | VIM | 0.928514 | 1.352016 | 1.20E-05 | 0.002822 | Down |
1464 | CSPG4 | −0.49526 | −0.39281 | 1.57E-05 | 0.003477 | Down |
55905 | RNF114 | 0.739501 | 1.17908 | 1.67E-05 | 0.00364 | Down |
57447 | NDRG2 | 0.552197 | 1.037042 | 2.20E-05 | 0.004259 | Down |
3696 | ITGB8 | −1.35342 | −0.49282 | 4.02E-05 | 0.006208 | Down |
6533 | SLC6A6 | −0.32443 | −0.23874 | 4.86E-05 | 0.007051 | Down |
2357 | FPR1 | −0.56296 | −0.17572 | 8.85E-05 | 0.01113 | Down |
The heat map of top 200 DEGs. Diagram presents the result of a two-way hierarchical clustering of top 200 DEGs and samples. The clustering is constructed using the complete-linkage method together with the Euclidean distance. Each row represents a DEG and each column, a sample. The DEG clustering tree is shown on the right. The color scale illustrates the relative level of DEG expression: red, below the reference channel; green, higher than the reference.
To investigate the biological function of the identified DEGs in breast cancer bone metastasis, GO term and KEGG pathway enrichment analyses was performed. In GO term and KEGG pathway enrichment analyses, cell adhesion, angiogenesis, and signal transduction were the most significant enrichment in biological process; protein binding, calcium ion binding and collagen binding were the most remark enrichment in molecular function; extracellular matrix, extracellular region, and plasma membrane were the most significant enrichment in cellular component. ECM-receptor interaction, focal adhesion, and osteoclast differentiation were the most significant enrichment in KEGG signaling pathways. Top 25 GO terms of DEGs were presented in
The functional enrichment of GO terms and KEGG pathways for differential mRNA.
To obtain the interaction between the proteins encoded by DEGs and other proteins, PPI network was explored and visualized by Cytoscape. PPI networks of proteins encoded by DEGs (FDR < 0.05) and interacting proteins encoded by DEGs (
The PPI networks. The node denoted protein, and the edge denoted the interactions. The rose red and green ellipse presented the up- and down-regulated proteins encoded by genes, respectively. The rectangle presented high degree (degree > 5) proteins encoded by DEGs. The solid line and dotted line presented direct interaction and colocalization, respectively.
In order to investigate the TFs-target genes regulatory network for bone metastasis of breast cancer, we utilized TRANSFAC to obtain 53 TFs (from proteins encoded by all DEGs) regulating other DEGs (
The TFs-target genes networks. In the network, the hexagon and rectangle presented the TFs and DEGs, respectively. The rose red and green ellipse presented the up- and down-regulation, respectively.
In this study, eight candidate genes were selected from the top 10 up- or down-regulated DEGs, which were IBSP, MMP9, MMP13, TNFAIP6, CD200, DHRS3, ASS1, RIPK4, VIM, and PROM1 for validation of integrated analysis results (
Identification of differential mRNA by QRT-PCR using metastatic bone, liver, and lung tissues.
In order to access the discriminatory ability of the above 10 genes used for qRT-PCR validation in breast cancer bone metastasis, ROC curve analyses were conducted and AUC were calculated. As shown in
The ROC curve analyses of IBSP, MMP9, MMP13, TNFAIP6, CD200, DHRS3, ASS1, RIPK4, VIM, and PROM1 in breast cancer bone metastasis. The x axis and y axis presented specificity and sensitivity, respectively.
To evaluate the impact of these predicted genes on cell migration and invasion, the wound healing assay and matrigel invasion assay were employed. As shown in
Wound healing assay and transwell invasion assay.
Among numerous malignancies, metastatic breast cancer is the second leading cause of woman death (
In this study, a total of 749 DEGs were identified, consisting of 469 up-regulated genes and 280 down-regulated genes and eight DEGs were selected from the top 10 up- or down-regulated DEGs, such as IBSP, MMP9, MMP13, TNFAIP6, DHRS3, ASS1, RIPK4, and VIM for the following qRT-PCR validation. CD200 and PROM1 was also selected for qRT-PCR validation as cluster of differentiation markers. Although there is no statistical difference except CD200 and DHRS3, the tendency of expression was generally consist with the bioinformatics analysis results. In addition, their diagnostic ability was also evaluated, and it appeared that IBSP, MMP9, MMP13, TNFAIP6, CD200, DHRS3, ASS1, RIPK4, and VIM may be considered as specific prognostic markers for bone metastasis to identify primary breast cancers that have potential to metastasize to bone tissue.
There are some preclinical evidences of factors associated with breast cancer cells homing to bone tissue. Awolaran et al. pointed out that 15 proteins expressed by breast cancer cells mediated breast cancer bone metastasis: ICAM-1, CDH11, osteoactivin, bone sialoprotein, CCN3, IL-11, etc. (
In our study, we detected some published potential biomarker of breast cancer bone metastasis, or their function has been reported in breast cancer bone metastasis such as SPARC, IBSP, MMP9, MMP13. SPARC was found as a DEG in breast cancer bone metastasis, while wasn't in the list of the up- or down-regulated DEGs. SPARC, also known as osteonectin, is involved in the chemoattraction of both breast and prostate cancer, and has been detected as a biomarker for bone metastasis (
IBSP, an osteoblast differentiation marker, is associated with osteoblastogenesis, bone formation and tumor-induced osteolysis (
Matrix metallopeptidase 9 (MMP9) and matrix metallopeptidase 13 (MMP13), members of the peptidase M10 family of matrix metalloproteinases (MMPs), is normally expressed by resident bone cells primarily by osteoclasts and osteoblasts (
Besides some published potential biomarkers of breast cancer bone metastasis were identified in our study, we also detected some novel biomarkers or their relationship with breast cancer bone metastasis has never been reported, such as TNFAIP6, DHRS3, ASS1, RIPK4, VIM, CD200. However, their function in breast cancer has been investigated except RIPK4.
Cancer cells can create an inflammatory microenvironment to enhance the tumor metastatic process (
Dehydrogenase/reductase 3 (DHRS3) has been found to be involved in the tumor suppressive pathway and constitutively expressed in breast cancer cell lines (
For RIPK4, their relationship to breast cancer still remains unreported. Receptor interacting serine/threonine kinase 4 (RIPK4) is involved in a number of signaling pathways and plays an important role in regulating cells proliferation, cell differentiation and cell apoptosis (
According to the KEGG pathway enrichment analysis, we found two important signaling pathways including osteoclast differentiation and rheumatoid arthritis in the bone metastasis of breast cancer. Bone metastasis of breast cancer induces severe osteolysis with increased bone resorption (
To further confirm the results of bioinformatics prediction, the wound healing assay and matrigel invasion assay were performed. Total five genes were selected, two are up-regulated in breast cancer and three are down-regulated in breast cancer. As shown in
In summary, we found several DEGs including IBSP, MMP9, MMP13, TNFAIP6, DHRS3, ASS1, RIPK4, VIM, CD200, and PROM1 may play important roles in the process of bone metastasis from breast cancer. Among which, IBSP, MMP9, TNFAIP6, DHRS3, RIPK4, and CD200 had a diagnose value for patients with breast cancer bone metastasis. Additionally, osteoclast differentiation and rheumatoid arthritis signaling pathways and related DEGs including FOS and CTSK may be involved in the development of bone metastasis from breast cancer. However, there are limitations to our study. Firstly, sample size in the QRT-PCR was small and a larger-scale study is further needed. Secondly, the biological function of identified DEGs wasn't investigated. Experiments in animal model and cell culture are further needed to explore underlying functional mechanism of breast cancer bone metastasis.
YZ is in charge of data collection and data mining. WH is in charge of molecular biological experiments and revising manuscript. SZ is in charge of data analysis, molecular biological experiments, and manuscript draft.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.