Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis,Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer

Abstract

Background:

Bone metastasis is a leading cause of the high mortality rate of prostate cancer (PCa), but curative strategies remain lacking. Recent studies suggest long non-coding RNAs (lncRNAs) may be potential targets to develop drugs. However, PCa bone metastasis-specifically-related lncRNAs were rarely reported. This study aimed to identify crucial lncRNAs and reveal their function mechanisms.

Methods:

GSE32269 and GSE26964 microarray datasets, downloaded from the Gene Expression Omnibus database, were used to analyze differentially expressed genes (DEGs)/lncRNAs (DELs) and miRNAs (DEMs), respectively. Weighted gene co-expression network analysis was performed to screen PCa bone metastasis-associated modules. The co-expression and competing endogenous RNAs (ceRNAs) networks were constructed to identify hub lncRNAs. Univariate Cox regression analysis was conducted to determine their prognostic values. The correlation of lncRNAs with immune infiltrating cells was analyzed by using Tumor IMmune Estimation Resource. Therapeutic drugs were predicted by querying the Connectivity Map (CMap) and the Comparative Toxicogenomics Database (CTD).

Results:

A total of 18 DELs, 2,614 DEGs and 86 DEMs were screened between bone metastatic and primary PCa samples. Four modules enriched by DEGs were shown to be bone metastasis-associated. LncRNA HCG18 and MCM3AP-AS1 were identified to be important because they existed in both of the co-expression and ceRNA networks (forming the relationship pairs: HCG18/MCM3AP-AS1-KNTC1, MCM3AP-AS1-hsa-miR-508-3p-DTL and HCG18/MCM3AP-AS1-hsa-miR-127-3p-CDKN3). All the genes in these interaction pairs were significantly associated with overall survival of PCa patients. Also, HCG18, MCM3AP-AS1 and their target mRNAs were positively correlated with various tumor-infiltrated immune cells, especially increased M2 macrophages. Valproic acid and trichostatin A may be effective to treat PCa bone metastasis by targeting HCG18 and MCM3AP-AS1.

Conclusion:

HCG18 and MCM3AP-AS1 that regulate M2 macrophage infiltration may be important targets to treat PCa bone metastasis and improve prognosis.

Keywords

prostate cancer bone metastasis lncRNA prognosis tumor-infiltrating

Introduction

Prostate cancer (PCa) is a leading cause of cancer-related death among men, with 174,650 new cases and 31,620 deaths estimated to occur in 2019 in the United States.¹ Metastatic growth in distant organs, especially bone (the most common metastatic site)^2,3 is the major cause of the high mortality rate. Although there have several strategies (including surgery, radiotherapy, chemotherapy, androgen deprivation therapy, radionuclide, diphosphate, zoledronic acid, et al.) recommended for the treatment of PCa bone metastasis,^4,5 most of them are palliative, not curative. Therefore, it is urgent to explore new therapeutic approaches for this disease.

With the developments in sequencing technology, recent scholars devote themselves to explain the mechanisms of PCa bone metastasis from the molecular levels and propose some concepts for targeted therapy. Proteins are complex macromolecules that perform all the specific functions in cells. Thus, messenger RNAs (mRNAs) that encode proteins were previously considered as the main study targets.⁶ The expression of Myc-associated zinc-finger protein (MAZ) was found to be upregulated in PCa tissues with bone metastasis compared with those without bone metastasis.⁷ High expression level of MAZ was positively correlated with poor overall and bone metastasis-free survival in PCa patients.⁷ These findings indicated that silencing of MAZ may be a potential therapeutic method to repress bone metastasis of PCa cells, which was proved by subsequent in vitro and in vivo experiments.⁷ Kaplan-Meier analysis of The Cancer Genome Atlas (TCGA) dataset revealed that overall survival (OS) was significantly shorter in patients who had the high expression level of G Protein-Coupled Receptor Class C Group 5 Member A (GPRC5A).⁸ Hence, GPRC5A may also be a possible therapeutic target, which was demonstrated in xenograft mice models (that is, prostate cancer cells with GPRC5A knockout failed to establish bone metastasis in mice).⁸

Recently, several studies suggest that long non-coding RNAs (lncRNAs, > 200 bp in length) contribute to the progression of PCa by influencing the expression of protein-encoding mRNAs via a model of co-expression or competing endogenous RNAs (ceRNAs) for competing microRNA (miRNAs).^9,10 LncAMPC was identified to be highly expressed in tumor tissues of metastatic PCa patients compared with localized PCa.¹⁰ Suppression of lncAMPC by small interfering RNAs (siRNAs) reduced proliferation, migration and invasion capacities of PCa cells.¹⁰ Mechanism studies showed that lncAMPC upregulated the expression of leukemia inhibitory factor (LIF) by sponging miR-637 and inhibiting its activity in the cytoplasm; while enhanced the transcription of LIF receptor (LIFR) by decoying histone H1.2 in the nucleus. Increased expressions of LIF/LIFR were beneficial to maintain the stability of programmed death-ligand 1 (PD-L1) protein which promoted CD8 T-cell exhaustion and led to immunosuppression.¹⁰ NAP1L6 expression was observed to be upregulated in PCa tissues and cell lines compared with controls.¹¹ Compared with PCa patients having a low level of NAP1L6 expression, those with high NAP1L6 expression had significantly increased distant metastasis rate and reduced OS rate.¹¹ siRNA-mediated knockdown of NAP1L6 significantly inhibited the migration and invasion of PCa cells by decreasing the expression of inhibin-β A.¹¹ Overexpressed SNHG7 in PCa was revealed to act as a sponge for miR-324-3p to result in the elevated expression of miR-324-3p target gene WNT2B, which facilitated the migration and invasion of PCa cells by inducing epithelial-mesenchymal transition (EMT).¹² These findings indicate lncRNAs may represent potential targets to develop therapeutic drugs for PCa bone metastasis. However, PCa bone metastasis-specifically-related lncRNAs were rarely reported, except for a recently reported PCAT7-miR-324-5p-transforming growth factor beta receptors I (TGFBR1) ceRNA axis.¹³

In this study, we aimed to: 1) screen PCa bone metastasis-associated protein-encoding mRNAs by weighted gene co-expression network analysis (WGCNA)¹⁴; 2) identify crucial lncRNAs based on the co-expression and ceRNA network analyses; 3) extract the hub protein-encoding mRNAs, lncRNAs and miRNAs by survival analysis; 4) predict the potential drugs that may treat PCa bone metastasis by reversing the expression levels of target genes; and 5) further reveal the mechanisms of lncRNAs by correlating with immune infiltration.^15,16

Materials and Methods

Data Collection

Microarray expression profiles of GSE32269 and GSE26964 were downloaded from the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) database (http://www.ncbi.nlm.nih.gov/geo/). GSE32269 contained 22 localized and 29 bone metastatic PCa samples, which analyzed the mRNA expression profiles based on the platform of Affymetrix Human Genome U133A Array (GPL96).¹⁷ GSE26964 included 6 primary and 7 bone metastatic PCa samples, which analyzed the miRNA expression profiles based on the platform of Capitalbio mammal microRNA V3.0 (GPL8469). In addition, the mRNA and miRNA expression profiles, as well as clinical survival information of 490 prostate adenocarcinoma (PRAD) patients were collected from the TCGA database (https://portal.gdc.cancer.gov/).

Screening of Differentially Expressed RNAs

The gene symbols of lncRNAs and mRNAs in the GSE32269 dataset were re-annotated using the data management system BioMart (http://asia.ensembl.org/biomart/martview/05d7d26eb91dd919eed0287709c9acbb). The differentially expressed genes (DEGs) and lncRNAs (DELs) in the GSE32269 dataset and differentially expressed miRNAs (DEMs) in the GSE26964 dataset were identified between bone metastatic and primary PCa samples using the linear models for microarray data (LIMMA) package in R (version 3.34.7; https://bioconductor.org/packages/release/bioc/html/limma.html).¹⁸ The statistical threshold of DEGs, DELs and DEMs was set as false discovery rate (FDR) < 0.05 and |log₂FC(fold change)| > 1. Hierarchical clustering analysis was performed to reveal the gene expression differences of 2 sample groups using package pheatmap in R (version: 1.0.8; https://cran.r-project.org/web/packages/pheatmap).

Identification of Modules by WGCNA Analysis

WGCNA¹⁴ is a popular systems biology algorithm that not only could effectively identify highly-connected and function-shared genes (in the same module) from thousands of genes in the microarray or sequencing data based on the correlation between genes due to their similar expression profile, but also could associate co-expression modules with clinical characteristics. Thus, WGCNA is a widely applied method to screen specific phenotype-associated genes in many cancers, including prostate cancer (castration-resistant-related,¹⁹ Gleason score-related²⁰ and lymph node metastasis-related²¹). Similar to these studies, we also used the WGCNA package in R (version 1.61; https://cran.r-project.org/web/packages/WGCNA/index.html)¹⁴ to identify modules highly correlated with bone metastasis of PCa in the GSE32269 dataset. Briefly, the soft threshold power (β) was chosen based on the scale-free topology criterion. The hierarchical clustering dendrogram was created by calculating the topological overlap matrix-based dissimilarity measure between genes. The modules were obtained with the DynamicTreeCut algorithm²² by defining the minimum module size of 100 and the minimum cut height of 0.99. The DEGs identified in the GSE32269 dataset were mapped to the above modules screened by WGCNA analysis. The enrichment of DEGs in each module was assessed by using a hypergeometric-based test [f(k, N, M, n) = C(k, M)*C(n-k, N-M)/C(n, N)].²³ Those with p < 0.05 and fold enrichment > 1 were considered as significant modules, the genes in which may be hub genes. Finally, the association between module eigengenes (the expression profiles of module genes) and clinical phenotypes was investigated according to Pearson’s correlation tests to reveal the importance of modules.

Construction of an lncRNA-mRNA Co-Expression Network

The DEGs (enriched in the significant modules identified by WGCNA) were selected to construct a co-expression network with all DELs to identify specific lncRNAs that regulated bone metastasis-related mRNAs. The Pearson correlation coefficients (PCC) between the expression values of each DEL-DEG pair across samples in the GSE32269 dataset were calculated using the cor.test function (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.test.html) in R. The DEL-DEG co-expression pairs with an absolute value of PCC > 0.4 were used to construct a DEL-DEG co-expression network, which was visualized by Cytoscape software (version 3.6.1; www.cytoscape.org/).

Construction of an lncRNA-miRNA-mRNA ceRNA Network

It has been reported that lncRNAs can serve as ceRNAs to competitively bind with miRNAs via miRNA response elements (MREs) and inhibit their expressions. The sequestered miRNAs could not interact with mRNAs that harbor common MREs with lncRNAs and thereby prevent the degradation or translational inhibition of mRNAs.^24,25 Theoretically, the expression trend of lncRNAs and mRNAs were consistent, but opposite to miRNAs. Based on this regulatory theory among lncRNAs, miRNAs and mRNAs, we constructed a ceRNA network by the following steps: 1) lncRNA-miRNA interaction prediction: the DELs (identified in GSE32269 dataset) and DEMs (identified in the GSE26964 dataset) were mapped into the lncRNAs-miRNAs interactions deposited in the miRwalk database (version 2.0; http: //www.zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2).²⁶ Only the relationship pairs with the opposite expression change trend of DELs and DEMs were retained for further analysis; 2) miRNA-mRNA interaction prediction: the target genes of DEMs that were interacted with DELs were further predicted with the miRwalk database. Likewise, only the target genes with the opposite expression change trend to DEMs and enriched in the significant modules were selected; and 3) construction of the ceRNA network: the ceRNA crosstalk network was constructed with the previously mentioned interactions between lncRNAs-miRNAs and between miRNAs-mRNAs and visualized through the Cytoscape software.

Function Enrichment Analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (version 6.8; http://david.abcc.ncifcrf.gov)²⁷ was used to predict the biological functions of DEGs in the co-expression and ceRNA networks. The Gene Ontology (GO) biological process terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched, with p-value < 0.05 defined as the statistical threshold.

Prognosis Analysis

Using the TCGA data, DEGs, DELs and DEMs in the co-expression and ceRNA networks were subjected to univariate Cox regression analysis with “survival” package (version 2.41 -1; http://bioconductor.org/packages/survivalr/) to identify candidate genes that were associated with OS. A p-value < 0.05 was considered significant. Subsequently, the samples were divided into a high-expression group and a low-expression group according to the median expression values of each gene and then the Kaplan-Meier survival curve was plotted.

Small Molecule Drug Analysis

The module DEGs were uploaded into the Connectivity Map (CMap) database (https://portals.broadinstitute.org/cmap/) to identify small molecules (showing negative enrichment scores) that may be potential drugs to treat PCa patients with bone metastasis. Furthermore, the Comparative Toxicogenomics Database (CTD) (http://ctdbase.org) was also searched using the name of selected hub genes to obtain the chemical-gene interaction pairs. The shared small molecular drugs identified in the CMap database and the CTD database, were considered to be especially important. The corresponding interaction pairs were used to construct the gene-drug network.

Association Between Gene Expression and Tumor-Infiltrating Immune Cells

The associations between the expression levels of selected hub genes and the abundance of 6 immune infiltrates of PRAD (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) were explored using the online tool TIMER (Tumor IMmune Estimation Resource; version 1.0, https://cistrome.shinyapps.io/timer/). In addition, correlations between hub genes and gene markers of tumor-infiltrating immune cells²⁸ were analyzed using Spearman’s correlation. CIBERSORT, TIDE, EPIC, quanTIseq and xCell algorithms in TIMER (version 2.0; http://timer.cistrome.org) were used to investigate the association with macrophage subtype. A p-value < 0.01 was defined as a threshold value.

Results

Identification of DEGs, DELs and DEMs

A total of 10,217 (including 95 lncRNAs and 10,122 protein-encoding mRNAs) were re-annotated in the GSE32269 dataset. Of them, 2,632 (including 18 lncRNAs and 2,614 mRNAs) were identified to be differentially expressed between bone metastatic and primary PCa samples using the LIMMA method (Table S1; Figure 1A). In the analysis of GSE26964 dataset, 86 of 546 miRNAs were found to meet the statistical criteria of FDR < 0.05 and |log2FC| > 1 (Table S1; Figure 1B). Hierarchical clustering analysis based on the expression levels of these DEGs, DELs and DEMs revealed that the PCa samples in the GSE32269 (Figure 1C) and GSE26964 (Figure 1D) datasets could be clustered into 2 different groups.

Figure 1.

Identification of differentially expressed RNAs in bone metastatic vs primary PCa samples. A, C: volcano plot (A) and heatmap (C) of differentially expressed lncRNAs and mRNAs identified in the GSE32269 dataset; B, D: volcano plot (B) and heatmap (D) of differentially expressed miRNAs identified in the GSE26964 dataset. Red, high expression; blue, low expression. FC, fold change; FDR, false discovery rates.

Screening of PCa Bone Metastasis- and Prognosis-Associated DEGs, DELs and DEMs

WGCNA was firstly performed on all the mRNAs in the GSE32269 dataset. As shown in Figure 2A and B, power = 6 was chosen as the soft-thresholding to fit a scale-free network distribution (scale-free R² = 0.9; mean connectivity = 1). Twelve modules were clustered with thresholds of gene number ≥ 100 and CutHeight = 0.99 (Figure 2C). Then, the DEGs were enriched into each module using a hypergeometric-based test. As a result, 4 modules (blue, brown, green and yellow) were predicted to be PCa bone metastasis-associated because of fold enrichment > 1 and p < 0.05 (Table 1). The DEGs included in these 4 modules are displayed in Table S2. Furthermore, the module-trait relationship analysis showed that the blue (r = -0.73), yellow (r = -0.72), green (r = 0.54) and brown (r = 0.72) modules were strongly correlated with the phenotype of patients (Figure 2D), indicating the DEGs in these modules may be important genes for bone metastasis of PCa.

Figure 2.

Weighted gene co-expression network analysis. A: selecting the soft-threshold power β (= 6) based on the scale-free topology threshold of 0.9; B: displaying the mean connectivity (= 1) when β = 6; C: showing the clustering dendrogram of genes; D, showing the correlations between gene modules and clinical traits.

Table 1.

Crucial Modules Screened by WGCNA.

ID	Color	Module size	#DEmRNAs	Enrichment infor
ID	Color	Module size	#DEmRNAs	Enrichment fold[95%CI]	P_hyper
Module 1	Black	151	26	0.484[0.305-0.739]	3.44E-04
Module 2	Blue	292	213	2.048[1.695-2.473]	8.76E-14
Module 3	Brown	261	172	1.851[1.507-2.269]	3.75E-09
Module 4	Green	228	134	1.651[1.315-2.066]	1.41E-05
Module 5	Greenyellow	111	20	0.506[0.297-0.822]	3.49E-03
Module 6	Gray	1287	293	0.639[0.555-0.734]	5.55E-11
Module 7	Magenta	132	53	1.128[0.800 -1.570]	4.99E-01
Module 8	Pink	147	45	0.860[0.599 -1.213]	4.06E-01
Module 9	Purple	116	36	0.871[0.580 -1.282]	5.15E-01
Module 10	Red	178	29	0.458[0.297-0.683]	3.57E-05
Module 11	Turquoise	2008	639	0.894[0.804-0.992]	3.44E-02
Module 12	Yellow	252	178	1.984[1.616-2.431]	4.36E-11

To reveal the functions of DELs, we constructed co-expression and ceRNA networks based on their regulatory relationships with DEGs in the above 4 modules. By calculation of the PCC, 18 DELs were found to interact with 444 DEGs, which constituted 903 interaction pairs (Figure 3; Table S3). Furthermore, 17 DELs were predicted to interact with 32 DEMs, while 21 of 32 DEMs were predicted to regulate 155 DEGs. Based on 21 overlapped DEMs, an lncRNA-miRNA-mRNA ceRNA network was constructed (Figure 4; Table S4).

Figure 3.

The co-expression relationships between differentially expressed lncRNAs and differentially expressed mRNAs identified in blue, brown, green and yellow modules.

Figure 4.

Construction of a ceRNA network using the differentially expressed lncRNAs, miRNAs and mRNAs in crucial modules.

Univariate Cox regression analysis of the TCGA data showed that 7 DELs (INE1, SNHG20, LINC00342, PVT1, HCG18, MCM3AP-AS1 and LINC01399), 7 DEMs (hsa-miR-301b, hsa-miR-508, hsa-miR-326, hsa-miR-127, hsa-miR-330, hsa-miR-532 and hsa-miR-331) and 73 DEGs in the co-expression and ceRNA networks were significantly associated with OS (Table S5; Figure 5).

Figure 5.

The target relationships between differentially expressed RNAs and small molecular drugs.

Comparison of the mRNAs regulated by 7 OS-related DELs in the co-expression network with OS-related DEGs showed that 17 were common (ACYP1, ACPP, PNN, SRSF7, CDKN3, HSPE1, CDC6, KLK3, KPNA2, SPDL1, KIF2C, SHMT2, CYTH2, VWA1, CKS2, KNTC1 and CPE). Of them, 12 (except ACPP, CPE, VWA1, KLK3 and CYTH2) had the same expression trend with 6 DELs (HCG18, SNHG20, LINC00342, PVT1, LINC01399 and MCM3AP-AS1) (all upregulated), indicating the corresponding co-expression interaction axes (HCG18-KNTC1, LINC00342-SPDL1/SHMT2, LINC01399-KPNA2/ACYP1/CDKN3/PNN/CKS2/SRSF7/HSPE1, MCM3AP-AS1-KIF2C/KNTC1, PVT1-KIF2C, SNHG20-KIF2C/CDC6) were correlative to PCa bone metastasis and prognosis. In the ceRNA network, only 6 OS-related DELs (INE1, SNHG20, LINC00342, PVT1, HCG18 and MCM3AP-AS1; HR > 1, upregulated) could interact with 3 OS-related DEMs (hsa-miR-326, hsa-miR-508-3p, hsa-miR-127-3p; HR < 1, downregulated). Seven target genes of 2 miRNAs (hsa-miR-508-3p, hsa-miR-127-3p) were also risk factors for poor OS (CDKN3, CDC6, DTL, ZWINT, SMC2, RAD51AP1, CDC25A; HR > 1, upregulated). These findings suggested the corresponding ceRNA interaction axes (MCM3AP-AS1-hsa-miR-508-3p-RAD51AP1/CDC6/CDC25A/DTL/SMC2, HCG18/INE1/LINC00342/PVT1/SNHG20/MCM3AP-AS1-hsa-miR-127-3p-CDKN/CDC25A/ZWINT) were also important candidate biomarkers for PCa bone metastasis and prognosis.

Function Enrichment Analysis

Function enrichment analysis showed that DEGs in the co-expression network were enriched into 98 GO biological process terms and 16 KEGG pathways (Table 2), in which only KPNA2, ACYP1 and HSPE1 genes were not enriched. This result indicated that the other 9 genes in the crucial co-expression axes may be especially pivotal for PCa bone metastasis and prognosis. The related GO biological process terms and KEGG pathways included GO:0051301∼cell division (KNTC1, KIF2C, CDC6, SPDL1, CKS2), GO:0000082∼G1/S transition of mitotic cell cycle (CDC6, CDKN3), GO:0007155∼cell adhesion (PNN), GO:0008284∼positive regulation of cell proliferation (SHMT2), hsa05200: Pathways in cancer (CKS2) and hsa05168: Herpes simplex infection (SRSF7).

Table 2.

Function Enrichment for the Genes in the Co-Expression Network.

Category	Term	P-value	Genes
GO BP	GO:0051301∼cell division	1.21E-12	KIFC1, NEK2, USP9X, KNTC1, NEDD9, AURKA, PTTG1, KIF2C, NUMA1, SAC3D1, NCAPG2, NCAPG, BUB1, FBXO5, STAG2, CDCA3, ERCC6 L, KIF14, CDK1, CDC6, KIF11, TPX2, SPDL1, CENPE, CDC20, CDC25C, NCAPD3, SMC4, CCNB1, CCNB2, MAD2L1, SPAG5, CKS2, BUB1B, FBXL7, SMC1A
GO BP	GO:0007067∼mitotic nuclear division	7.17E-11	NEK2, USP9X, KNTC1, NEDD9, AURKA, PTTG1, CEP55, KIF2C, NUMA1, SAC3D1, NCAPG2, BUB1, FBXO5, ASPM, STAG2, CDCA3, ERCC6 L, CDC6, CDK1, KIF11, KIF15, TPX2, CDC20, PBK, CDC25C, CCNB2, BUB1B, FBXL7
GO BP	GO:0000070∼mitotic sister chromatid segregation	8.38E-08	KIFC1, MAD2L1, CENPA, NEK2, SPAG5, NUSAP1, ESPL1, SMC1A, SMC4
GO BP	GO:0030198∼extracellular matrix organization	3.14E-07	OLFML2B, COL3A1, ITGA10, OLFML2A, SPINT1, SPARC, SOX9, COL5A2, COL5A1, CSGALNACT1, LAMB3, FBLN1, LAMA4, BGN, LAMB2, COL1A2, VCAN, VWA1, LAMB1, SPP1
GO BP	GO:0007062∼sister chromatid cohesion	1.22E-06	CENPM, KNTC1, SPDL1, CDC20, CENPE, KIF2C, MAD2L1, CENPA, BUB1, BUB1B, CENPU, SMC1A, STAG2, ERCC6L
GO BP	GO:0007080∼mitotic metaphase plate congression	2.32E-06	KIF14, CCNB1, KIFC1, KIF2C, CHMP4A, PIBF1, CENPE, SPDL1, CEP55
GO BP	GO:0008283∼cell proliferation	3.44E-06	TSPAN1, POLA1, ENPEP, KIF2C, TYMS, MCM7, COL4A3BP, PEMT, BUB1, CSF1 R, CDK1, DLGAP5, KIF15, PAK1IP1, TPX2, SKP2, CDC25C, STAT3, CSGALNACT1, CHRM3, DLX5, CKS2, BUB1B, FBXL7, TXNRD1, MELK
GO BP	GO:0007018∼microtubule-based movement	2.54E-05	KIF14, KIFC1, KIF2C, KIF4A, KIF11, KIF15, KIF5C, CENPE, RACGAP1, KIF20A, KIF13B
GO BP	GO:0007051∼spindle organization	2.86E-05	KIZ, KIF11, SPAG5, TTK, AURKA, ASPM
GO BP	GO:0000086∼G2/M transition of mitotic cell cycle	2.90E-05	CDK1, NEK2, FOXM1, SKP2, TPX2, AURKA, CDC25C, CCNB1, PRKAR2B, CCNB2, CSNK1E, FGFR1OP, FBXL7, MELK
GO BP	GO:0045926∼negative regulation of growth	7.17E-05	ING4, MT1 M, MT1E, MT1 H, MT1G, MT1F
GO BP	GO:0000082∼G1/S transition of mitotic cell cycle	1.83E-04	CDK1, CDC6, TYMS, MCM7, RRM2, POLA1, SKP2, FBXO5, ORC6, CDKN3, MCM4
GO BP	GO:0000083∼regulation of transcription involved in G1/S transition of mitotic cell cycle	1.91E-04	CDK1, CDC6, TYMS, RRM2, POLA1, FBXO5
GO BP	GO:0006936∼muscle contraction	2.72E-04	ACTG2, DES, ACTA2, COL4A3BP, LMOD1, TPM2, SSPN, VIPR1, TPM1, MYLK, MYL9
GO BP	GO:0019882∼antigen processing and presentation	3.47E-04	AZGP1, RAB3B, MICA, HLA-C, HLA-E, CTSH, PSMB8, HLA-G
GO BP	GO:0000281∼mitotic cytokinesis	6.00E-04	KIF4A, CENPA, NUSAP1, CEP55, RACGAP1, KIF20A
GO BP	GO:0031145∼anaphase-promoting complex-dependent catabolic process	6.46E-04	CCNB1, CDK1, MAD2L1, SKP2, BUB1B, AURKA, CDC20, PTTG1, PSMB8
GO BP	GO:0007093∼mitotic cell cycle checkpoint	9.58E-04	MAD2L1, KNTC1, BUB1, TTK, BUB1B, SMC1A
GO BP	GO:0071294∼cellular response to zinc ion	9.87E-04	MT1 M, MT1E, MT1 H, MT1G, MT1F
GO BP	GO:0045638∼negative regulation of myeloid cell differentiation	9.87E-04	ZFP36, MEIS2, WDR61, ITPKB, PAF1
GO BP	GO:0002480∼antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-independent	0.001065	HLA-C, HLA-E, HLA-G, HLA-F
GO BP	GO:0034501∼protein localization to kinetochore	0.001494	CDK1, TTK, BUB1B, SPDL1
GO BP	GO:0001501∼skeletal system development	0.001871	BGLAP, DLX5, MMP9, COL3A1, COL1A2, VCAN, NPR3, PAPSS1, SOX9, COL5A2, SERP1
GO BP	GO:0007010∼cytoskeleton organization	0.001918	SYNE3, DES, APOE, PPL, KRT15, NEDD9, KRT13, OBSL1, DPYSL2, MSN, SOX9, TPM1
GO BP	GO:0007155∼cell adhesion	0.001926	BGLAP, MYBPC1, FERMT1, ITGA10, CD99, NEDD9, ACKR3, CDH2, SSPN, CD151, COL5A1, PNN, AZGP1, LAMB3, LAMA4, OMD, LAMB2, SORBS2, VCAN, LAMB1, THBS2, SPP1, ADAM9
GO BP	GO:0042493∼response to drug	0.00316	CDK1, BGLAP, ADCY1, HMGB2, DPYSL2, FOSB, PNP, STAT3, CCNB1, FOS, TYMS, ACSL1, MCM7, JUN, PEMT, ABCC4, SRD5A1
GO BP	GO:0007566∼embryo implantation	0.00332	FBLN1, MMP9, TRIM28, PCSK5, MMP2, SPP1
GO BP	GO:0006890∼retrograde vesicle-mediated transport, Golgi to ER	0.003726	KIF2C, KIF4A, KIF11, KIF15, ARCN1, TMED10, CENPE, RACGAP1
GO BP	GO:0042476∼odontogenesis	0.003837	FGFR2, BGLAP, OSR2, COL1A2, LAMB1
GO BP	GO:0030574∼collagen catabolic process	0.004481	CTSK, MMP9, COL3A1, COL1A2, MMP2, COL5A2, COL5A1
GO BP	GO:0006260∼DNA replication	0.004567	CDK1, GINS2, ING4, CDC6, MCM7, DTL, RRM2, POLA1, ORC6, CDC25C, MCM4
GO BP	GO:0051384∼response to glucocorticoid	0.004838	TYMS, BGLAP, PAPPA, SPARC, ANXA3, GHR, ADAM9
GO BP	GO:0007076∼mitotic chromosome condensation	0.005177	NCAPG, NUSAP1, NCAPD3, SMC4
GO BP	GO:0030335∼positive regulation of cell migration	0.005333	WNT5B, SYNE2, F3, FGFR1OP, ELP5, HBEGF, HSPA5, LRRC15, LAMB1, CTSH, MYLK, CSF1R
GO BP	GO:0002474∼antigen processing and presentation of peptide antigen via MHC class I	0.005659	PDIA3, HLA-C, HLA-E, HLA-G, HLA-F
GO BP	GO:0007052∼mitotic spindle organization	0.005659	CCNB1, KIF11, TTK, AURKA, SMC1A
GO BP	GO:0051603∼proteolysis involved in cellular protein catabolic process	0.005938	CTSK, CTSO, CLPP, HSPA5, CTSH, PSMB8
GO BP	GO:0007059∼chromosome segregation	0.006035	KIF11, HJURP, NEK2, SPAG5, USP9X, CENPE, TOP2A
GO BP	GO:0007179∼transforming growth factor beta receptor signaling pathway	0.006996	FOS, LTBP2, USP9X, JUN, COL3A1, COL1A2, GDF15, ADAM9
GO BP	GO:0006270∼DNA replication initiation	0.00715	CDC6, MCM7, POLA1, ORC6, MCM4
GO BP	GO:0051436∼negative regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle	0.007434	CCNB1, CDK1, MAD2L1, FBXO5, BUB1B, CDC20, PSMB8
GO BP	GO:0071276∼cellular response to cadmium ion	0.007464	MT1E, MT1 H, MT1G, MT1F
GO BP	GO:0032873∼negative regulation of stress-activated MAPK cascade	0.008241	FOXM1, FOXO1, PBK
GO BP	GO:0030705∼cytoskeleton-dependent intracellular transport	0.008798	KIF14, MYO10, KIF5C, KIF13B
GO BP	GO:0034976∼response to endoplasmic reticulum stress	0.009646	UFL1, HYOU1, HERPUD1, HSP90B1, PDIA3, COL4A3BP, MBTPS1
GO BP	GO:0022617∼extracellular matrix disassembly	0.010265	CTSK, LAMB3, KLK2, MMP9, MMP2, ADAMTS5, SPP1
GO BP	GO:0051437∼positive regulation of ubiquitin-protein ligase activity involved in regulation of mitotic cell cycle transition	0.010265	CCNB1, CDK1, MAD2L1, FBXO5, BUB1B, CDC20, PSMB8
GO BP	GO:0051726∼regulation of cell cycle	0.01068	CCNB1, ICK, CCNB2, DTL, JUN, FOXM1, SKP2, CDC25C, STAT3
GO BP	GO:0090307∼mitotic spindle assembly	0.010858	KIFC1, KIF11, NEK2, TPX2, PIBF1
GO BP	GO:0051549∼positive regulation of keratinocyte migration	0.011353	MMP9, HBEGF, ADAM9
GO BP	GO:0051988∼regulation of attachment of spindle microtubules to kinetochore	0.011353	NEK2, SPAG5, RACGAP1
GO BP	GO:1902895∼positive regulation of pri-miRNA transcription from RNA polymerase II promoter	0.01186	FOS, JUN, STAT3, SMARCA4
GO BP	GO:0007094∼mitotic spindle assembly checkpoint	0.01186	MAD2L1, BUB1, TTK, BUB1B
GO BP	GO:0001568∼blood vessel development	0.013103	LAMA4, COL1A2, FOXO1, PLPP3, COL5A1
GO BP	GO:0009612∼response to mechanical stimulus	0.013999	CCNB1, BGLAP, MEIS2, JUN, COL3A1, FOSB
GO BP	GO:0030199∼collagen fibril organization	0.014328	FMOD, COL3A1, COL1A2, COL5A2, COL5A1
GO BP	GO:1903071∼positive regulation of ER-associated ubiquitin-dependent protein catabolic process	0.014896	HERPUD1, UBQLN2, SGTA
GO BP	GO:0009629∼response to gravity	0.014896	FOS, BGLAP, SPARC
GO BP	GO:1901215∼negative regulation of neuron death	0.015623	HTRA2, APOE, CHMP4A, STAT3, GHR
GO BP	GO:0009636∼response to toxic substance	0.017155	CDK1, TYMS, FOS, SLC30A4, EPHX2, SLC18A2, CNP
GO BP	GO:0051439∼regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle	0.017467	CCNB1, CDK1, FBXO5, CDC20
GO BP	GO:0035338∼long-chain fatty-acyl-CoA biosynthetic process	0.018428	ACSL1, ELOVL5, ACOT1, HACD1, TECR
GO BP	GO:0036500∼ATF6-mediated unfolded protein response	0.018848	HSP90B1, HSPA5, MBTPS1
GO BP	GO:0060337∼type I interferon signaling pathway	0.019342	IRF9, HLA-C, HLA-E, PSMB8, HLA-G, HLA-F
GO BP	GO:0007088∼regulation of mitotic nuclear division	0.01961	NEK2, FBXO5, OBSL1, CDC25C
GO BP	GO:0006268∼DNA unwinding involved in DNA replication	0.023186	MCM7, MCM4, TOP2A
GO BP	GO:0090501∼RNA phosphodiester bond hydrolysis	0.023186	AZGP1, ANG, RNASE6
GO BP	GO:2000009∼negative regulation of protein localization to cell surface	0.023186	NEDD4 L, ASTN2, TAX1BP3
GO BP	GO:0032870∼cellular response to hormone stimulus	0.023189	FOS, JUN, FOSB, STAT3, GHR
GO BP	GO:0006997∼nucleus organization	0.024314	NUMA1, CHMP4A, CEP55, GOLM1
GO BP	GO:0040014∼regulation of multicellular organism growth	0.024314	FGFR2, HTRA2, STAT3, GHR
GO BP	GO:0035987∼endodermal cell differentiation	0.026873	LAMB3, MMP9, LAMB1, MMP2
GO BP	GO:0009314∼response to radiation	0.029571	JUN, COL3A1, KRT13, SMC1A
GO BP	GO:0016032∼viral process	0.030998	SATB1, FBLN1, KRT18, CENPA, POLA1, BUB1, HLA-C, ACKR3, CENPU, CDC25C, SGTA, PSMB8, STAT3, NCKAP1
GO BP	GO:0045444∼fat cell differentiation	0.032008	STEAP4, WNT5B, NR4A2, FOXO1, NR4A1, PSMB8
GO BP	GO:0071872∼cellular response to epinephrine stimulus	0.03294	PDE4B, SRD5A1, PDE4D
GO BP	GO:0036499∼PERK-mediated unfolded protein response	0.03294	HERPUD1, ATF3, HSPA5
GO BP	GO:0043206∼extracellular fibril organization	0.03294	LTBP2, COL3A1, COL5A1
GO BP	GO:0071277∼cellular response to calcium ion	0.034771	FOS, ADCY1, JUN, ITPKB, FOSB
GO BP	GO:0008284∼positive regulation of cell proliferation	0.035662	KLF5, FGFR2, KIF14, SHMT2, CCPG1, FOXM1, TTK, CDC20, VIPR1, SOX9, STAT3, UFL1, OSR2, ATF3, FGFR1OP, HBEGF, CTSH, DPP4, CSF1R
GO BP	GO:0006508∼proteolysis	0.036419	DPP3, TMPRSS2, HPN, KLK2, PDIA3, KLK3, MMP9, AGBL5, ANPEP, ESPL1, MMP2, PRSS8, CTSK, HTRA2, PAPPA, CTSO, CTSH, MBTPS1, ADAMTS5, DPP4
GO BP	GO:0034097∼response to cytokine	0.036974	TYMS, FOS, JUN, COL3A1, SPARC
GO BP	GO:0060021∼palate development	0.037156	TCF21, DHRS3, OSR2, DLX5, PAK1IP1, TIPARP
GO BP	GO:0006629∼lipid metabolic process	0.037647	NAAA, APOF, PTPRN2, PEMT, PLA1A, ACAD8, TECR, PLPP3, ASAH1
GO BP	GO:0042761∼very long-chain fatty acid biosynthetic process	0.038316	ELOVL5, HACD1, TECR
GO BP	GO:0036109∼alpha-linolenic acid metabolic process	0.038316	ACSL1, ELOVL5, SCP2
GO BP	GO:0034975∼protein folding in endoplasmic reticulum	0.038316	HSP90B1, PDIA3, HSPA5
GO BP	GO:0007049∼cell cycle	0.03914	EID1, DIXDC1, KRT18, CCPG1, HJURP, RB1CC1, FOXM1, AURKA, CDC20, CDKN3, ING1
GO BP	GO:0043923∼positive regulation by host of viral transcription	0.043999	HPN, JUN, SMARCA4
GO BP	GO:0001503∼ossification	0.044769	CHRDL1, MMP9, SPARC, SOX9, COL5A2, SPP1
GO BP	GO:0021987∼cerebral cortex development	0.046573	KIF14, COL3A1, SRD5A1, CDH2, ASPM
GO BP	GO:0044236∼multicellular organism metabolic process	0.047889	TIPARP, GHR
GO BP	GO:1903225∼negative regulation of endodermal cell differentiation	0.047889	COL5A2, COL5A1
GO BP	GO:0007057∼spindle assembly involved in female meiosis I	0.047889	FBXO5, AURKA
GO BP	GO:1901898∼negative regulation of relaxation of cardiac muscle	0.047889	PDE4B, PDE4D
GO BP	GO:0001837∼epithelial to mesenchymal transition	0.048607	FGFR2, TRIM28, SOX9, HNRNPAB
GO BP	GO:0043588∼skin development	0.048607	COL3A1, COL5A2, COL5A1, DHCR24
GO BP	GO:0006271∼DNA strand elongation involved in DNA replication	0.049972	GINS1, GINS2, POLA1
KEGG	hsa04110: Cell cycle	5.91E-08	CDC6, CDK1, SKP2, TTK, CDC20, ESPL1, PTTG1, CDC25C, MCM4, CCNB1, MAD2L1, CCNB2, MCM7, BUB1, BUB1B, ORC6, SMC1A, STAG2
KEGG	hsa04114: Oocyte meiosis	1.34E-05	CDK1, ADCY1, CDC20, ESPL1, AURKA, PTTG1, CDC25C, CCNB1, CCNB2, MAD2L1, BUB1, FBXO5, PPP3CC, SMC1A
KEGG	hsa04512: ECM-receptor interaction	1.59E-04	LAMA4, LAMB3, LAMB2, COL3A1, COL1A2, ITGA10, LAMB1, THBS2, COL5A2, COL5A1, SPP1
KEGG	hsa05146: Amoebiasis	0.003011	ADCY1, LAMA4, IL1R1, LAMB3, LAMB2, COL3A1, COL1A2, LAMB1, COL5A2, COL5A1
KEGG	hsa05166: HTLV-I infection	0.005513	ZFP36, IL1R1, ADCY1, WNT5B, HLA-C, CDC20, PTTG1, HLA-E, HLA-G, HLA-F, FOS, ATF3, MAD2L1, JUN, PPP3CC, BUB1B
KEGG	hsa04510: Focal adhesion	0.005561	COL3A1, ITGA10, COL5A2, COL5A1, MYL9, LAMB3, LAMA4, LAMB2, JUN, COL1A2, LAMB1, THBS2, MYLK, SPP1
KEGG	hsa04940: Type I diabetes mellitus	0.006311	CPE, PTPRN2, HLA-C, HLA-E, HLA-G, HLA-F
KEGG	hsa04978: Mineral absorption	0.007698	MT1 M, MT1E, TRPV6, MT1 H, MT1G, MT1F
KEGG	hsa01040: Biosynthesis of unsaturated fatty acids	0.026169	ELOVL5, ACOT1, HACD1, TECR
KEGG	hsa00230: Purine metabolism	0.027605	ADCY1, RRM2, PDE4B, ENPP4, POLA1, PDE4D, ENTPD4, PAPSS1, GMPR, PNP, CANT1
KEGG	hsa05200: Pathways in cancer	0.029381	FGFR2, ADCY1, WNT5B, KLK3, MMP9, CYCS, SKP2, FOXO1, MMP2, STAT3, FOS, LAMB3, HSP90B1, LAMA4, LAMB2, JUN, CKS2, LAMB1, CSF1R
KEGG	hsa00062: Fatty acid elongation	0.032616	ELOVL5, ACOT1, HACD1, TECR
KEGG	hsa05222: Small cell lung cancer	0.033052	LAMA4, LAMB3, LAMB2, CYCS, CKS2, SKP2, LAMB1
KEGG	hsa05168: Herpes simplex infection	0.034803	IRF9, CDK1, FOS, SRSF7, JUN, CYCS, SKP2, HLA-C, HLA-E, HLA-G, HLA-F
KEGG	hsa04914: Progesterone-mediated oocyte maturation	0.036454	CCNB1, CDK1, ADCY1, MAD2L1, CCNB2, BUB1, CDC25C
KEGG	hsa04115: p53 signaling pathway	0.040761	CCNB1, CDK1, CCNB2, RRM2, CYCS, RPRM

Function enrichment analysis of DEGs in the ceRNA network showed that 63 GO biological process terms and 7 KEGG pathways were enriched (Table 3), in which only hub gene RAD51AP1 was not enriched, but HSPE1 was included. Hereby, the 6 genes in the crucial ceRNA axes may be especially key for PCa bone metastasis and prognosis. The related GO biological process terms and KEGG pathways included GO:0051301∼cell division (CDC6, SMC2, CDC25A, ZWINT), GO:0000082∼G1/S transition of mitotic cell cycle (CDC6, CDKN3, CDC25A), GO:0006260∼DNA replication (CDC6, DTL, CDC25A), GO:0007093∼mitotic cell cycle checkpoint (ZWINT) and hsa04110: Cell cycle (CDC6, CDC25A).

Table 3.

Function Enrichment for the Genes in the ceRNA Network.

Category	Term	P-value	Genes
GO BP	GO:0051301∼cell division	4.83E-16	KIF14, CDK1, CDC6, KIF11, USP9X, TPX2, CENPE, AURKA, CDC25C, TACC3, SMC2, CDC37, CDC25A, SMC4, CCNB1, KIF2C, CCNB2, MAD2L1, NCAPG2, ZWINT, CKS2, BUB1, FBXO5, FBXL7, SMC1A, NUP43
GO BP	GO:0007067∼mitotic nuclear division	6.64E-11	CDC6, CDK1, KIF11, USP9X, KIF15, TPX2, AURKA, PBK, CEP55, CDC25C, CDC25A, KIF2C, CCNB2, NCAPG2, BUB1, FBXO5, FBXL7, NUP43
GO BP	GO:0000086∼G2/M transition of mitotic cell cycle	2.38E-10	CDK1, SKP2, TPX2, AURKA, CDC25C, CDC25A, HMMR, CCNB1, PRKAR2B, CCNB2, FGFR1OP, FBXL7, MELK, DNM2
GO BP	GO:0008283∼cell proliferation	5.51E-10	CDK1, KIF15, POLA1, SKP2, TPX2, PRKDC, CDC25C, TACC3, CDC25A, KIF2C, EPS8, DLX5, BUB1, CKS2, PCNA, FBXL7, TXNRD1, IGFBP4, MELK, EMP1
GO BP	GO:0000082∼G1/S transition of mitotic cell cycle	2.87E-07	CDK1, CDC6, PCNA, POLA1, SKP2, FBXO5, MCM2, CDKN3, CDC25A, MCM6
GO BP	GO:0006260∼DNA replication	1.15E-06	CDK1, GINS2, CDC6, DTL, PCNA, POLA1, MCM2, RBM14, CDC25C, CDC25A, MCM6
GO BP	GO:0000070∼mitotic sister chromatid segregation	2.29E-06	MAD2L1, CENPA, ZWINT, NUSAP1, SMC1A, SMC4
GO BP	GO:0007062∼sister chromatid cohesion	3.50E-05	KIF2C, MAD2L1, CENPA, ZWINT, BUB1, CENPE, SMC1A, NUP43
GO BP	GO:0000083∼regulation of transcription involved in G1/S transition of mitotic cell cycle	4.50E-05	CDK1, CDC6, PCNA, POLA1, FBXO5
GO BP	GO:0007018∼microtubule-based movement	7.98E-05	KIF14, KIF2C, KIF11, KIF15, CENPE, RACGAP1, KIF20A
GO BP	GO:0000281∼mitotic cytokinesis	1.16E-04	CENPA, NUSAP1, CEP55, RACGAP1, KIF20A
GO BP	GO:0007052∼mitotic spindle organization	1.33E-04	CCNB1, KIF11, TTK, AURKA, SMC1A
GO BP	GO:0007093∼mitotic cell cycle checkpoint	1.72E-04	MAD2L1, ZWINT, BUB1, TTK, SMC1A
GO BP	GO:0006271∼DNA strand elongation involved in DNA replication	2.82E-04	GINS1, GINS2, PCNA, POLA1
GO BP	GO:0007080∼mitotic metaphase plate congression	3.05E-04	KIF14, CCNB1, KIF2C, CENPE, CEP55
GO BP	GO:0000079∼regulation of cyclin-dependent protein serine/threonine kinase activity	3.75E-04	CDC6, CDC25C, CDKN3, CDC37, CDC25A
GO BP	GO:0051726∼regulation of cell cycle	8.06E-04	CCNB1, CCNB2, E2F5, DTL, SKP2, CDC25C, CDC25A
GO BP	GO:0019886∼antigen processing and presentation of exogenous peptide antigen via MHC class II	0.001321	KIF2C, KIF11, KIF15, CENPE, RACGAP1, DNM2
GO BP	GO:0021987∼cerebral cortex development	0.001491	KIF14, COL3A1, SRD5A1, CDH2, TACC3
GO BP	GO:0006977∼DNA damage response, signal transduction by p53 class mediator resulting in cell cycle arrest	0.002174	CCNB1, CDK1, PCNA, AURKA, CDC25C
GO BP	GO:0001649∼osteoblast differentiation	0.002271	DLX5, TNC, GJA1, VCAN, HSPE1, GPNMB
GO BP	GO:0030574∼collagen catabolic process	0.002443	COL4A2, COL3A1, COL15A1, MMP2, COL5A1
GO BP	GO:0006270∼DNA replication initiation	0.002759	CDC6, POLA1, MCM2, MCM6
GO BP	GO:0006268∼DNA unwinding involved in DNA replication	0.003315	MCM2, TOP2A, MCM6
GO BP	GO:0016925∼protein sumoylation	0.003778	PCNA, TDG, SMC6, SMC1A, TOP2A, NUP43
GO BP	GO:0032467∼positive regulation of cytokinesis	0.004186	KIF14, CDC6, RACGAP1, ECT2
GO BP	GO:0043206∼extracellular fibril organization	0.004806	COL3A1, MFAP5, COL5A1
GO BP	GO:0031145∼anaphase-promoting complex-dependent catabolic process	0.00522	CCNB1, CDK1, MAD2L1, SKP2, AURKA
GO BP	GO:0006890∼retrograde vesicle-mediated transport, Golgi to ER	0.005955	KIF2C, KIF11, KIF15, CENPE, RACGAP1
GO BP	GO:0007077∼mitotic nuclear envelope disassembly	0.006822	CCNB1, CDK1, CCNB2, NUP43
GO BP	GO:0007076∼mitotic chromosome condensation	0.007515	NUSAP1, SMC2, SMC4
GO BP	GO:1904874∼positive regulation of telomerase RNA localization to Cajal body	0.007515	CCT5, CCT6A, SHQ1
GO BP	GO:0030198∼extracellular matrix organization	0.007863	COL4A2, LUM, TNC, COL3A1, VCAN, MFAP5, COL5A1
GO BP	GO:0007051∼spindle organization	0.008539	KIF11, TTK, AURKA
GO BP	GO:0042752∼regulation of circadian rhythm	0.009192	ROCK2, EZH2, PRKDC, TOP2A
GO BP	GO:0048511∼rhythmic process	0.011987	ROCK2, EZH2, PRKDC, TOP2A
GO BP	GO:0007094∼mitotic spindle assembly checkpoint	0.013212	MAD2L1, BUB1, TTK
GO BP	GO:0070301∼cellular response to hydrogen peroxide	0.013874	CDK1, EZH2, PCNA, ECT2
GO BP	GO:0051439∼regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle	0.017293	CCNB1, CDK1, FBXO5
GO BP	GO:0007057∼spindle assembly involved in female meiosis I	0.01755	FBXO5, AURKA
GO BP	GO:0007155∼cell adhesion	0.020552	SRPX, TNC, COL15A1, NINJ1, VCAN, CDH2, ADAM12, GPNMB, COL5A1, PNN
GO BP	GO:0007059∼chromosome segregation	0.022162	KIF11, USP9X, CENPE, TOP2A
GO BP	GO:0007507∼heart development	0.023077	KCNJ8, RB1CC1, COL3A1, PCNA, GJA1, PRKDC
GO BP	GO:0051436∼negative regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle	0.024799	CCNB1, CDK1, MAD2L1, FBXO5
GO BP	GO:0009314∼response to radiation	0.025108	COL3A1, SMC1A, CDC25A
GO BP	GO:0051383∼kinetochore organization	0.02621	SMC2, SMC4
GO BP	GO:0006272∼leading strand elongation	0.02621	PCNA, POLA1
GO BP	GO:0097421∼liver regeneration	0.026814	EZH2, PCNA, AURKA
GO BP	GO:0021510∼spinal cord development	0.028566	EED, SRD5A1, DPYSL2
GO BP	GO:0051437∼positive regulation of ubiquitin-protein ligase activity involved in regulation of mitotic cell cycle transition	0.029556	CCNB1, CDK1, MAD2L1, FBXO5
GO BP	GO:0000722∼telomere maintenance via recombination	0.032203	PCNA, POLA1, SMC6
GO BP	GO:0051984∼positive regulation of chromosome segregation	0.034795	CDC6, SMC6
GO BP	GO:0045132∼meiotic chromosome segregation	0.034795	SMC2, SMC4
GO BP	GO:0007100∼mitotic centrosome separation	0.043304	KIF11, AURKA
GO BP	GO:0032354∼response to follicle-stimulating hormone	0.043304	PAPPA, SRD5A1
GO BP	GO:0071394∼cellular response to testosterone stimulus	0.043304	ROCK2, SRD5A1
GO BP	GO:0010032∼meiotic chromosome condensation	0.043304	SMC2, SMC4
GO BP	GO:0030199∼collagen fibril organization	0.046235	LUM, COL3A1, COL5A1
GO BP	GO:0000302∼response to reactive oxygen species	0.046235	PRDX6, CYCS, TXNRD1
GO BP	GO:0042787∼protein ubiquitination involved in ubiquitin-dependent protein catabolic process	0.046514	CCNB1, CDK1, MAD2L1, RMND5A, AURKA
GO BP	GO:0018108∼peptidyl-tyrosine phosphorylation	0.046514	BAZ1B, FGFR1OP, TTK, CDC37, MELK
GO BP	GO:0045184∼establishment of protein localization	0.048393	CORO1C, KIF14, CEP55
GO BP	GO:0043065∼positive regulation of apoptotic process	0.049935	HTRA2, PRKDC, ECT2, TOP2A, MELK, DNM2, SHQ1
KEGG	hsa04110: Cell cycle	1.63E-13	CDC6, CDK1, YWHAZ, E2F5, SKP2, TTK, PRKDC, MCM2, CDC25C, CDC25A, MCM6, CCNB1, MAD2L1, CCNB2, BUB1, PCNA, SMC1A
KEGG	hsa04114: Oocyte meiosis	2.76E-06	CCNB1, CDK1, YWHAZ, MAD2L1, CCNB2, BUB1, FBXO5, AURKA, CDC25C, SMC1A
KEGG	hsa04914: Progesterone-mediated oocyte maturation	4.01E-05	CCNB1, CDK1, MAD2L1, CCNB2, BUB1, CDC25C, PIK3R3, CDC25A
KEGG	hsa03030: DNA replication	0.006866	PCNA, POLA1, MCM2, MCM6
KEGG	hsa05222: Small cell lung cancer	0.013573	COL4A2, CYCS, CKS2, SKP2, PIK3R3
KEGG	hsa04512: ECM-receptor interaction	0.014681	COL4A2, TNC, COL3A1, COL5A1, HMMR
KEGG	hsa04115: p53 signaling pathway	0.036265	CCNB1, CDK1, CCNB2, CYCS

Identification of Small-Molecule Drugs

A total of 585 small molecules were predicted to have the potential to treat PCa patients with bone metastasis. Of them, 32 were found to target crucial DELs, DEMs and DEGs (Table 4; Figure 6). Ionomycin (DTL, CDKN3, KNTC1, SPDL1), irinotecan (CDKN3), resveratrol (CDKN3, KNTC1, SPDL1), etoposide (CDKN3), vorinostat (CDKN3) may be especially effective for the treatment of PCa patients with bone metastasis because their enrichment score approximated to -1. Furthermore, valproic acid (HCG18, MCM3AP-AS1) and trichostatin A (CDKN3, MCM3AP-AS1) may also be pivotal treatment drugs because they could target DELs.

Table 4.

Crucial Small Molecular Drugs Predicted by CMap.

Cmap name	Enrichment	P-value
ionomycin	-0.979	0.00004
irinotecan	-0.965	0.0001
colforsin	-0.898	0.00004
resveratrol	-0.794	0
etoposide	-0.741	0.00891
daunorubicin	-0.733	0.01026
vorinostat	-0.707	0
simvastatin	-0.689	0.02009
thioridazine	-0.686	0
miconazole	-0.685	0.00715
clofibrate	-0.653	0.24073
thapsigargin	-0.638	0.09651
niclosamide	-0.635	0.01676
clobetasol	-0.61	0.12827
methotrexate	-0.588	0.00379
ivermectin	-0.514	0.09244
zidovudine	-0.484	0.21383
glafenine	-0.432	0.33792
disulfiram	-0.389	0.33958
trichostatin A	-0.357	0
lovastatin	-0.341	0.63524
wortmannin	-0.299	0.06655
dexamethasone	-0.266	0.54038
Sirolimus	-0.242	—
Valproic Acid	-0.104	—
cantharidin	-0.837	—
doxorubicin	-0.364	—
sulindac sulfide	-0.783	—
tretinoin	-0.327	—
exemestane	-0.501	—
quercetin	-0.369	—
carbamazepine	-0.226	—

Figure 6.

Kaplan-Meier survival curves comparing the survival difference of the high-expression group and the low-expression group. The groups were classified according to the median expression values of each crucial lncRNA, miRNA and mRNA in prostate adenocarcinoma samples from the TCGA database. A: lncRNA HCG18; B: lncRNA MCM3AP-AS1; C: hsa-miR-508-3p; D: hsa-miR-127-3p: E: DTL; F: CDKN3; G: KNTC1. HR, hazard ratio.

Association of Hub Genes’ Expression With Immune Infiltration

To further explore the possible mechanisms of hub genes, we analyzed the association of their expressions with the abundance of tumor-infiltrated immune cells which were previously demonstrated to be contributors for PCa progression.^29,30 TIMER analysis showed that DTL and KNTC1 were significantly associated with all 6 immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) (Figure 7A and B). Among the 3 target genes of hsa-miR-127-3p, the FDR of CDKN3 was the minimum. The expression level of CDKN3 was strongly associated with the abundance of B cells, neutrophils, and dendritic cells (p < 0.001) (Figure 7D). Also, the expression levels of HCG18 and MCM3AP-AS1 were positively associated with the abundance of B cells, CD8+ T cells, neutrophils, macrophages and dendritic cells. Furthermore, we also found the expression levels of DTL (31), KNTC1 (35), CDKN3 (29), HCG18 (31) and MCM3AP-AS1 (23) were correlated with most of the markers of various immune cells (Table 5), especially CD163, VSIG4 and MS4A4A of M2 macrophages (all associated and these biomarkers were DEGs). Their associations with M2 macrophages were also confirmed using the estimations from CIBERSORT, quanTIseq, TIDE or XCELL (Table 6). Thus, HCG18/MCM3AP-AS1-KNTC1, HCG18/MCM3AP-AS1-hsa-miR-127-3p-CDKN3 and MCM3AP-AS1-hsa-miR-508-3p-DTL interaction axes may be involved in PCa bone metastasis by regulating immune microenvironment, especially M2 macrophages.

Figure 7.

The correlation between the expression levels of lncRNAs/mRNAs and immune infiltration in prostate adenocarcinoma samples from the TCGA database. A: DTL; B: KNTC1; C: CDKN3; D: HCG18; E: MCM3AP-AS1.

Table 5.

Correlation Analysis Between Hub Genes and Markers of Immune Cells.

		DTL		KNTC1		CDKN3		HCG18		MCM3APAS
		Cor	P	Cor	P	Cor	P	Cor	P	Cor	P
CD8+ T cell	CD8A	0.153	*	0.306	***	0.0797	0.105	0.206	***	0.167	**
	CD8B	0.0681	0.165	0.162	**	-0.0834	0.0892	0.214	***	0.0498	0.311
B cell	CD19	0.0360	0.464	0.113	0.0216	0.0611	0.214	-0.0127	0.797	-0.0196	0.691
	CD79A	0.112	0.0218	0.160	*	0.130	*	0.0518	0.292	-0.00387	0.937
TAM	CCL2	0.0169	0.731	0.047	0.339	0.0771	0.117	0.0634	0.197	-0.0436	0.375
	CD68	0.433	***	0.391	***	0.423	***	0.144	*	0.165	**
	IL10	0.236	***	0.285	***	0.224	***	0.210	***	0.176	**
M1 Macrophage	NOS2	0.0419	0.394	0.114	0.0205	-0.0817	0.0962	0.250	***	0.126	0.0101
	IRF5	0.451	***	0.517	***	0.429	***	0.248	***	0.249	***
	PTGS2	0.110	0.0255	0.165	**	-0.0638	0.194	0.276	***	0.0843	0.0860
M2 Macrophage	CD163^a	0.365	***	0.363	***	0.304	***	0.324	***	0.306	***
	VSIG4^a	0.371	***	0.360	***	0.326	***	0.258	***	0.230	***
	MS4A4A^a	0.369	***	0.363	***	0.378	***	0.186	**	0.178	**
Dendritic cells	HLA-DPB1	0.144	*	0.193	***	0.140	*	-0.00241	0.961	-0.0672	0.171
	HLA-DQB1	0.199	***	0.188	**	0.168	**	0.0246	0.616	-0.0058	0.906
	HLA-DRA	0.334	***	0.340	***	0.291	***	0.192	***	0.125	0.0108
	HLA-DPA1	0.296	***	0.317	***	0.225	***	0.166	**	0.112	0.0221
	CD1C	0.235	***	0.311	***	0.141	*	0.204	***	0.148	*
	NRP1	0.395	***	0.191	***	0.318	***	0.404	***	0.274	***
	ITGAX	0.348	***	0.381	***	0.493	***	0.189	**	0.255	***
Th1	TBX21	0.197	***	0.192	***	0.358	***	0.150	*	0.135	*
	STAT4^a	0.156	*	0.122	0.0125	0.401	***	0.212	***	0.220	***
	IFNG^a	0.203	***	0.182	**	0.273	***	0.125	0.0106	0.0623	0.205
	CXCR3	0.208	***	0.140	*	0.311	***	0.147	*	0.108	0.0272
	CCR5	0.307	***	0.233	***	0.429	***	0.273	***	0.255	***
Th2	GATA3	-0.00591	0.904	0.145	*	-0.173	**	0.107	0.0292	0.0361	0.462
	STAT6^b	0.00169	0.973	0.281	***	-0.171	**	0.395	***	0.325	***
	CCR4	0.380	***	0.469	***	0.240	***	0.483	***	0.396	***
Th17	STAT3^b	0.299	***	0.315	***	0.030	0.545	0.559	***	0.428	***
	IL6	-0.0615	0.211	-0.0429	0.383	-0.0231	0.639	0.0575	0.242	-0.0883	0.0722
	CCR6	0.243	***	0.313	***	0.193	***	0.300	***	0.207	***
	RORA	0.341	***	0.487	***	0.114	0.0205	0.589	***	0.432	***
tTreg	FOXP3	0.364	***	0.426	***	0.264	***	0.308	***	0.258	***
	TGFB1	0.235	***	0.302	***	0.235	***	0.136	*	0.0911	0.0635
	IL10	0.236	***	0.285	***	0.224	***	0.210	***	0.176	**
T-cell exhaustion	PDCD1^a	0.131	*	0.224	***	0.111	0.0239	0.104	0.0335	0.0329	0.5046
	CTLA4	0.150	*	0.311	***	0.200	***	0.0827	0.0921	0.0615	0.210
	LAG3	-0.0721	0.142	0.119	0.0147	-0.0685	0.1630	-0.167	**	-0.142	*
	HAVCR2	0.398	***	0.418	***	0.405	***	0.164	**	0.152	*
	GZMB	0.0615	0.210	0.103	0.0355	0.0704	0.152	0.0203	0.680	-0.0604	0.219
Neutrophils	CD66b(CEACAM8)	0.0598	0.223	0.0902	0.0660	0.0173	0.725	0.0517	0.293	0.0255	0.604
	CD11b(ITGAM)^a	0.2981	***	0.434	***	0.231	***	0.245	***	0.237	***
	CCR7	0.194	***	0.257	***	0.147	0.00267	0.166	**	0.0979	0.0460

* P < 0.01; **P < 0.001; ***P < 0.0001. Bold indicated the genes were differentially expressed between bone metastatic and primary PCa samples in the GSE32269 dataset. ^a, upregulated; ^b, downregulated.

Table 6.

Association With Macrophage Subtype.

	Infiltrates	rho	adj.p
CDKN3	Macrophage M0_CIBERSORT	0.0493	0.513
	Macrophage M0_CIBERSORT-ABS	0.114	0.0707
	Macrophage M1_CIBERSORT	0.0626	0.391
	Macrophage M1_CIBERSORT-ABS	0.180	**
	Macrophage M1_QUANTISEQ	-0.0109	0.921
	Macrophage M1_XCELL	0.309	***
	Macrophage M2_CIBERSORT	0.304	***
	Macrophage M2_CIBERSORT-ABS	0.356	***
	Macrophage M2_QUANTISEQ	0.00532	0.959
	Macrophage M2_TIDE	-0.0553	0.457
	Macrophage M2_XCELL	0.145	0.0140
	Macrophage_EPIC	0.345	***
	Macrophage_TIMER	0.00279	0.978
	Macrophage_XCELL	0.268	***
	Macrophage/Monocyte_MCPCOUNTER	0.0323	0.708
DTL	Macrophage M0_CIBERSORT	0.00469	0.967
	Macrophage M0_CIBERSORT-ABS	0.0742	0.236
	Macrophage M1_CIBERSORT	0.109	0.0727
	Macrophage M1_CIBERSORT-ABS	0.255	***
	Macrophage M1_QUANTISEQ	0.0105	0.897
	Macrophage M1_XCELL	0.105	0.0841
	Macrophage M2_CIBERSORT	0.265	***
	Macrophage M2_CIBERSORT-ABS	0.400	***
	Macrophage M2_QUANTISEQ	0.108	0.0735
	Macrophage M2_TIDE	-0.121	0.0425
	Macrophage M2_XCELL	0.0400	0.554
	Macrophage_EPIC	0.272	***
	Macrophage_TIMER	0.166	**
	Macrophage_XCELL	0.138	0.0173
	Macrophage/Monocyte_MCPCOUNTER	0.188	***
KNTC1	Macrophage M0_CIBERSORT	0.0319	0.674
	Macrophage M0_CIBERSORT-ABS	0.115	0.0523
	Macrophage M1_CIBERSORT	0.0516	0.446
	Macrophage M1_CIBERSORT-ABS	0.258	***
	Macrophage M1_QUANTISEQ	0.0215	0.808
	Macrophage M1_XCELL	-0.0131	0.874
	Macrophage M2_CIBERSORT	0.117	0.0480
	Macrophage M2_CIBERSORT-ABS	0.422	***
	Macrophage M2_QUANTISEQ	0.195	***
	Macrophage M2_TIDE	-0.205	***
	Macrophage M2_XCELL	-0.123	0.0377
	Macrophage_EPIC	0.264	***
	Macrophage_TIMER	0.130	0.0262
	Macrophage_XCELL	0.0250	0.763
	Macrophage/Monocyte_MCPCOUNTER	0.201	**
HCG18	Macrophage M0_CIBERSORT	-0.00481	0.944
	Macrophage M0_CIBERSORT-ABS	0.0401	0.577
	Macrophage M1_CIBERSORT	0.236	***
	Macrophage M1_CIBERSORT-ABS	0.311	***
	Macrophage M1_QUANTISEQ	0.0699	0.286
	Macrophage M1_XCELL	-0.311	***
	Macrophage M2_CIBERSORT	-0.125	0.0380
	Macrophage M2_CIBERSORT-ABS	0.146	0.0123
	Macrophage M2_QUANTISEQ	0.297	***
	Macrophage M2_TIDE	-0.221	***
	Macrophage M2_XCELL	-0.241	***
	Macrophage_EPIC	-0.0487	0.483
	Macrophage_TIMER	0.102	0.106
	Macrophage_XCELL	-0.175	*
	Macrophage/Monocyte_MCPCOUNTER	0.367	***
MCM3AP-AS1	Macrophage M0_CIBERSORT	0.0140	0.876
	Macrophage M0_CIBERSORT-ABS	0.0647	0.365
	Macrophage M1_CIBERSORT	0.165	*
	Macrophage M1_CIBERSORT-ABS	0.222	***
	Macrophage M1_QUANTISEQ	-0.0689	0.325
	Macrophage M1_XCELL	-0.239	***
	Macrophage M2_CIBERSORT	0.0230	0.786
	Macrophage M2_CIBERSORT-ABS	0.173	*
	Macrophage M2_QUANTISEQ	0.288	***
	Macrophage M2_TIDE	-0.0552	0.452
	Macrophage M2_XCELL	-0.152	0.0115
	Macrophage_EPIC	-0.00908	0.918
	Macrophage_TIMER	0.165	*
	Macrophage_XCELL	-0.0910	0.164
	Macrophage/Monocyte_MCPCOUNTER	0.227	***

* P < 0.01; **P < 0.001; ***P < 0.0001.

Discussion

In the present study, we identified HCG18 and MCM3AP-AS1 as PCa bone metastasis- and poor prognosis-associated lncRNAs. HCG18 and MCM3AP-AS1 may function by directly co-expressing with KNTC1 to regulate its transcription or indirectly modulating the expression of CDKN3 by sponging hsa-miR-127-3p. Furthermore, MCM3AP-AS1 also acted as a molecular sponge for hsa-miR-508-3p to increase the expression of DTL.

There had several studies to show that HCG18 was upregulated in cancer tissues and cells, including nasopharyngeal carcinoma,³¹ hepatocellular carcinoma,³² colorectal cancer,³³ gastric cancer³⁴ and lung adenocarcinoma.³⁵ A higher level of HCG18 was associated with positive lymph node metastasis,³¹ tumor node metastasis stage, invasion depth³⁶ and poor prognosis^31,34 of cancer patients. Knockdown of HCG18 inhibited the proliferation, migration, invasion, metastasis, while induced the apoptosis of cancer cells.^31
-33,36,37 Also, silencing of HCG18 suppressed the growth of xenografts in mice.³⁶ Mechanistically, HCG18 was reported to competitively bind to miR-140,³¹ miR-214-3p,³² miR-1271,³³ miR-34a-5p³⁵ and miR-152-3p³⁴ to respectively enhance the expressions of the target genes of these miRNAs, including cyclin D1, CENPM, MTDH, HMMR and DNAJB12. Until now, no studies reported the association of HCG18 with PCa. Also, the direct co-expression mechanisms with mRNAs were not investigated. In this study, we, for the first time, predicted HCG18 may be involved in PCa bone metastasis and poor prognosis by influencing the expression of KNTC1 and hsa-miR-127-3p-CDKN3. KNTC1 (kinetochore associated 1) encodes a protein that ensures proper chromosome segregation during cell division. Thus, it may mainly influence tumor cell growth. This hypothesis was demonstrated by a study in which knockdown of KNTC1 effectively inhibited cell viability and increased apoptosis.³⁸ Although the functions of KNTC1 for metastasis were not explored, the roles of associated genes [such as spindle and kinetochore-associated protein (SKA)] may indirectly explain its possible mechanisms. SKA1 was reported to be significantly correlated with clinical stage, extrathyroid invasion³⁹ and poor prognosis.⁴⁰ Knockdown of SKA1 repressed the ability of cell proliferation, migration and invasion of PCa⁴¹ and glioma.⁴⁰ SKA2 expression was also increased in breast cancer tissues and cells, and associated with lymph node metastasis. SKA2 knockdown significantly reduced the migration and invasion by reversing EMT.⁴² CDKN3 had been a known gene highly expressed in PCa.⁴³ Knockdown of CDKN3 significantly attenuated cell invasion, induced G1 phase arrest and increased apoptosis rates in PCa cells⁴³; while cells with increased expression of CDKN3 exhibited elevated migratory and invasive potentials.⁴⁴ miR-127-3p was directly proved to be reduced in bone metastasis-positive PCa tissues relative to that in bone metastasis-negative PCa tissues.⁴⁵ Overexpression of miR-127-3p inhibited the invasion and migration of PCa cells.⁴⁵ In line with these studies, we also found HCG18, KNTC1 and CDKN3 were upregulated, while miR-127-3p was downregulated in bone metastasis samples compared with those without bone metastasis. Patients with high levels of HCG18, KNTC1 and CDKN3, but low expression of miR-127-3p had a shorter OS. Hence, HCG18-KNTC1 and HCG18-hsa-miR-127-3p-CDKN3 may be crucial targets for developing therapeutic approaches.

MCM3AP-AS1 was previously identified as a PCa-associated lncRNA. It was abnormally upregulated in prostate cancer tissues compared with controls.^46,47 Lentivirus-mediated overexpression of MCM3AP-AS1 promoted the proliferation, invasion, and migration, while suppressing the apoptosis of PCa cells. Mechanism study revealed that MCM3AP-AS1 may promote PCa progression by inducing methylation of NPY1 R promoter and then downregulating the expression of NPY1 R, which subsequently activated the MAPK pathway.⁴⁷ MCM3AP-AS1-miR-455-EGFR,⁴⁸ MCM3AP-AS1-miR-340-5p-KPNA4,⁴⁹ MCM3AP-AS1-miR-211-5p-SPARC,⁵⁰ MCM3AP-AS1-miR-138-5p-FOXK1⁵¹ and MCM3AP-AS1-miR-28-5p-CENPF⁵² ceRNA axes were also identified to participate in cancer metastasis, but whether this mechanism was related to PCa bone metastasis had not been confirmed. In this study, we, for the first time, predicted MCM3AP-AS1 may be involved in PCa bone metastasis by directly interacting with KNTC1 or acting as a sponge of hsa-miR-508-3p or hsa-miR-127-3p to regulate DTL and CDKN3. The roles of KNTC1 and hsa-miR-127-3p-CDKN3 were described as above. DTL (denticleless E3 ubiquitin protein ligase homolog) is an ubiquitin-related protein that was reported to contribute to the migration, invasion and poor prognosis of cancer⁵³ by inducing the degradation of programmed cell death 4.⁵⁴ Low expressed miR-508-3p was demonstrated to be significantly associated with distant metastasis.⁵⁵ Ectopic expression of miR-508-3p significantly suppressed the invasion ability of cancer cells.^55,56 However, no studies validated the roles of hsa-miR-508-3p and DTL in PCa, indicating MCM3AP-AS1-hsa-miR-508-3p-DTL axis may be a novel target for explaining its pathogenesis and developing targeted drugs.

A previous study indicated that HCG18 was an immune-related lncRNA to predict the poor prognosis for patients with anaplastic gliomas.⁵⁷ To further explore the biological functions of our identified lncRNAs and their mRNA targets, we also performed the association analysis with infiltrating immune cells in PCa tissue samples. Our results showed that all of HCG18, MCM3AP-AS1 and their target mRNAs were positively correlated with an increased count of M2 macrophage as well as its corresponding marker genes (CD163, VSIG4 and MS4A4A). There is evidence to support that men with high numbers of CD163-positive M2 macrophages had an increased risk to die of PCa compared with those with a low number.⁵⁸ Higher expression of macrophage markers (CD163 and VSIG4) was associated with higher rates of biochemical relapse in aged prostate cancer patients.⁵⁹ Treatment of tumor cells with M2 type macrophages decreased the susceptibility of PCa cells to natural killer (NK) cell cytotoxicity and may lead to tumor progression.⁶⁰ Pretreatment to reduce M2-like macrophages suppressed skeletal metastatic tumor growth.⁶¹ Accordingly, we concluded that our identified lncRNAs and their mRNA targets may promote PCa bone metastasis by increasing M2 macrophages.

Based on the CMap and CTD database analysis, we predicted histone deacetylases (HDAC) inhibitors valproic acid and trichostatin A may be potentially effective for the treatment of PCa bone metastasis. Our prediction results can be validated according to published studies. Previous evidence had shown that valproic acid and trichostatin A inhibited the proliferation and induced the apoptosis of PCa cells,^62,63 with the mechanisms of inactivating EGFR-STAT3 pathway and downregulating cell cycle-related genes (cyclin D1 and CDK6, CDK1, cyclin B).^62,64 Valproic acid and trichostatin A were also reported to induce the reversal process of EMT (downregulation of Slug, SMAD4 and up-regulation of E-cadherin), resulting in attenuated invasion and migration abilities in PC3 cells.^65
-67 However, all these studies about PCa focused on the mRNAs changed by HDAC inhibitors and did not explore lncRNAs. Our studies, for the first time, indicated valproic acid and trichostatin A may reverse the high expression of HCG18 and MCM3AP-AS1 in PCa. More importantly, HDAC inhibitor CG-745 was demonstrated to enhance the anti-cancer effect of anti-PD-1 antibody by changing the immune microenvironment, including increasing cytotoxic T cells and NK cells, while decreasing regulatory T cells and M2 macrophages.⁶⁸ Trichostatin A treatment made the HepG2 cells more susceptible to NK cell-mediated killing,⁶⁹ which may be accompanied by reduced M2 type macrophages.⁶⁰ In this study, HCG18 and MCM3AP-AS1 were identified to be associated with M2 type macrophages. Thus, valproic acid and trichostatin A may be effective to treat PCa bone metastasis by reducing the number of HCG18- and MCM3AP-AS1-mediated M2 type macrophages.

There are some limitations in this study. First, the sample size in the GSE32269 and GSE26964 datasets were relatively small. Second, the TCGA data did not provide information on bone metastasis. The bone metastasis-free survival could not be performed and only OS was investigated. Thus, the expression and survival-associations of DELs, DEMs and DEGs in the co-expression or ceRNA axes needed to be validated in a larger clinical patient cohort with PCa bone metastasis. Third, the co-expression or ceRNA relationship pairs were only predicted by software and preliminarily constructed according to their expression levels. In vitro and in vivo experiments (co-precipitation, knockdown or overexpression) were required to confirm their interactions. Fourth, the treatment mechanisms of valproic acid and trichostatin A to regulate lncRNAs and then influence the infiltration of immune cells should be further explored by some wet experiments (drug treatment, co-incubation of PCa with M1 or M2 macrophages, flow cytometry, transwell assay).^68,70

Conclusion

In summary, our study suggests lncRNAs HCG18 and MCM3AP-AS1 may be important targets to treat PCa metastasis and improve poor prognosis. They functioned by a co-expression (HCG18/MCM3AP-AS1-KNTC1) or ceRNA (MCM3AP-AS1-hsa-miR-508-3p-DTL and HCG18/MCM3AP-AS1-hsa-miR-127-3p-CDKN3) fashion. They also regulated M2 macrophage infiltration. Valproic acid and trichostatin A may be potential therapeutic drugs by reversing the expression levels of HCG18 and MCM3AP-AS1.

Supplemental Material

Supplemental Material, sj-xlsx-1-tct-10.1177_1533033821990064 - Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer

Supplemental Material, sj-xlsx-1-tct-10.1177_1533033821990064 for Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer by Yanfang Chen, Zheng Chen, Jian Mo, Mao Pang, Zihao Chen, Feng Feng, Peigen Xie and Bu Yang in Technology in Cancer Research & Treatment

Supplemental Material

Supplemental Material, sj-xlsx-2-tct-10.1177_1533033821990064 - Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer

Supplemental Material, sj-xlsx-2-tct-10.1177_1533033821990064 for Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer by Yanfang Chen, Zheng Chen, Jian Mo, Mao Pang, Zihao Chen, Feng Feng, Peigen Xie and Bu Yang in Technology in Cancer Research & Treatment

Supplemental Material

Supplemental Material, sj-xlsx-3-tct-10.1177_1533033821990064 - Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer

Supplemental Material, sj-xlsx-3-tct-10.1177_1533033821990064 for Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer by Yanfang Chen, Zheng Chen, Jian Mo, Mao Pang, Zihao Chen, Feng Feng, Peigen Xie and Bu Yang in Technology in Cancer Research & Treatment

Supplemental Material

Supplemental Material, sj-xlsx-4-tct-10.1177_1533033821990064 - Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer

Supplemental Material, sj-xlsx-4-tct-10.1177_1533033821990064 for Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer by Yanfang Chen, Zheng Chen, Jian Mo, Mao Pang, Zihao Chen, Feng Feng, Peigen Xie and Bu Yang in Technology in Cancer Research & Treatment

Supplemental Material

Supplemental Material, sj-xlsx-5-tct-10.1177_1533033821990064 - Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer

Supplemental Material, sj-xlsx-5-tct-10.1177_1533033821990064 for Identification of HCG18 and MCM3AP-AS1 That Associate With Bone Metastasis, Poor Prognosis and Increased Abundance of M2 Macrophage Infiltration in Prostate Cancer by Yanfang Chen, Zheng Chen, Jian Mo, Mao Pang, Zihao Chen, Feng Feng, Peigen Xie and Bu Yang in Technology in Cancer Research & Treatment

Footnotes

Authors’ Note

Yanfang Chen and Zheng Chen contributed to this work equally. The datasets generated for this study can be found in the GEO (http://www.ncbi.nlm.nih.gov/geo/; GSE32269 and GSE26964) and the TCGA-PRAD () databases. No experiments on humans and animals were used in this study. Therefore, no ethical approval is required.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Bu Yang, PhD

Supplemental Material

Supplemental material for this article is available online.

References

Siegel

Miller

Jemal

. Cancer statistics, 2019. CA Cancer J Clin. 2019;69(1):7–34.

Zhuo

Cheng

Pan

, et al. Prostate cancer with bone metastasis in Beijing: an observational study of prevalence, hospital visits and treatment costs using data from an administrative claims database. BMJ Open. 2019;9(6):e028214.

Liu

Kuai

Zhu

, et al. Prognosis of prostate cancer and bone metastasis pattern of patients: a SEER-based study and a local hospital based study from China. Sci Rep. 2020;10(1):9104.

Orcajo-Rincon

Caresia-Aróztegui

Del Puig Cózar-Santiago

, et al. Radium-223 in the treatment of bone metastasis in patients with castration-resistant prostate cancer. Review and procedure. Rev Esp Med Nucl Imagen Mol. 2018; 37(5):330–337.

Cheng

Zhuo

Pan

, et al. Treatment patterns of prostate cancer with bone metastasis in Beijing: a real-world study using data from an administrative claims database. Pharmacoepidemiol Drug Saf. 2019; 28(11):1501–1509.

Zhu

Wen

Xuan

, et al. Identifying the key genes and microRNAs in prostate cancer bone metastasis by bioinformatics analysis. FEBS Open Bio. 2020;10(4):674–688.

Yang

Lang

, et al. MAZ promotes prostate cancer bone metastasis through transcriptionally activating the KRas-dependent RalGEFs pathway. J Exp Clin Cancer Res. 2019;38(1):391.

Sawada

Kikugawa

Iio

, et al. GPRC5A facilitates cell proliferation through cell cycle regulation and correlates with bone metastasis in prostate cancer. Int J Cancer. 2020;146(5):1369–1382.

Pang

Zhao

Wang

, et al. Competing endogenous RNA and coexpression network analysis for identification of potential biomarkers and therapeutics in association with metastasis risk and progression of prostate cancer. Oxid Med Cell Longev. 2019;2019:8265958.

10.

Zhang

Shi

Chen

, et al. Novel long non-coding RNA lncAMPC promotes metastasis and immunosuppression in prostate cancer by stimulating LIF/LIFR expression. Mol Ther. 2020;S1525-S0016(20)30303–30308.

11.

Zheng

Gao

, et al. Long non-coding RNA NAP1L6 promotes tumor progression and predicts poor prognosis in prostate cancer by targeting Inhibin-β A. Onco Targets Ther 2018;11:4965–4977.

12.

Han

Zhou

. Knockdown of LncRNA SNHG7 inhibited epithelial-mesenchymal transition in prostate cancer though miR-324-3p/WNT2B axis in vitro. Pathol Res Pract. 2019;215(10):152537.

13.

Lang

Dai

Yang

Zhang

. SMAD3/SP1 complex-mediated constitutive active loop between lncRNA PCAT7 and TGF-β signaling promotes prostate cancer bone metastasis. Mol Oncol. 2020; 14(4):808–828.

14.

Langfelder

Horvath

. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9(1):559.

15.

Song

Chao

Zhuo

Chen

. Identification of hub genes in prostate cancer using robust rank aggregation and weighted gene co-expression network analysis. Aging (Albany NY). 2019;11(13):4736–4756.

16.

Wang

Peng

. UBASH3B Is a novel prognostic biomarker and correlated with immune infiltrates in prostate cancer. Front Oncol. 2019;9:1517.

17.

Cai

Wang

, et al. ERG induces androgen receptor-mediated regulation of SOX9 in prostate cancer. J Clin Invest. 2013;123(3):1109–1122.

18.

Ritchie

Phipson

, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res .2015;43(7):e47.

19.

Cheng

Qin

. Identification of castration-resistant prostate cancer-related hub genes using weighted gene co-expression network analysis. J Cell Mol Med. 2020; 24(14):8006–8017.

20.

Wang

Yang

. A Gleason score-related outcome model for human prostate cancer: a comprehensive study based on weighted gene co-expression network analysis. Cancer Cell Int. 2020;20:159.

21.

Chen

Lin

, et al. Development and validation of hub genes for lymph node metastasis in patients with prostate cancer. J Cell Mol Med. 2020;24(8):4402–4414.

22.

Langfelder

Zhang

Horvath

. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut Package for R. Bioinformatics. 2008;24(5):719–720.

23.

Cao

Zhang

. A Bayesian extension of the hypergeometric test for functional enrichment analysis. Biometrics. 2014;70(1):84–94.

24.

Chan

Tay

Noncoding

RNA

: RNA regulatory networks in cancer. Int J Mol Sci. 2018; 19(5):1310.

25.

Huo

Lin

. Computational identification of cross-talking ceRNAs. Adv Exp Med Biol. 2018;1094:97–108.

26.

Dweep

Gretz

. miRWalk2.0: a comprehensive atlas of microRNA-target interactions. Nat Methods. 2015;12(8):697–697.

27.

Dennis

Sherman

Hosack

, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol. 2003;4(5):P3.

28.

Shang

Sun

Liu

Integrated analysis of lymphocyte infiltration-associated lncRNA for ovarian cancer via TCGA, GTEx and GEO datasets. PeerJ 2020;8:e8961.

29.

Xiang

Jin

Yang

, et al. Infiltrating CD4+ T cells attenuate chemotherapy sensitivity in prostate cancer via CCL5 signaling. Prostate. 2019;79(9):1018–1031.

30.

Maolake

Izumi

Shigehara

, et al. Tumor-associated macrophages promote prostate cancer migration through activation of the CCL22-CCR4 axis. Oncotarget. 2017;8(6):9739–9751.

31.

Jiang

. LncRNA HCG18 contributes to nasopharyngeal carcinoma development by modulating miR-140/CCND1 and Hedgehog signaling pathway. Eur Rev Med Pharmacol Sci. 2019;23(23):10387–10399.

32.

Zou

Sun

. LncRNA HCG18 contributes to the progression of hepatocellular carcinoma via miR-214-3p/CENPM axis. J Biochem 2020;168(5):535–546. doi:10.1093/jb/mvaa073

33.

Zhang

Sun

Zhang

. The long non-coding RNA HCG18 promotes the growth and invasion of colorectal cancer cells through sponging miR-1271 and upregulating MTDH/Wnt/β-catenin. Clin Exp Pharmacol Physiol. 2020;47(4):703–712.

34.

Liu

Zhao

. HNF1A-induced lncRNA HCG18 facilitates gastric cancer progression by upregulating DNAJB12 via miR-152-3p. Onco Targets Ther. 2020;13:7641–7652.

35.

Pan

Jiang

Zhao

. HCG18/miR-34a-5p/HMMR axis accelerates the progression of lung adenocarcinoma. Biomed Pharmacother. 2020;129:110217.

36.

. Silencing of long non-coding RNA-HCG18 inhibits the tumorigenesis of gastric cancer through blocking PI3K/Akt pathway. Onco Targets Ther. 2020;13:2225–2234.

37.

Liu

Lin

Dong

, et al. Long noncoding RNA HCG18 up-regulates the expression of WIPF1 and YAP/TAZ by inhibiting miR-141-3p in gastric cancer. Cancer Med. 2020;9(18):6752–6765. doi:10.1002/cam4.3288

38.

Liu

Min

Wang

Zhang

. shRNA-mediated knockdown of KNTC1 suppresses cell viability and induces apoptosis in esophageal squamous cell carcinoma. Int J Oncol. 2019;54(3):1053–1060.

39.

Dong

Wang

. Expression of spindle and kinetochore-associated protein 1 is associated with poor prognosis in papillary thyroid carcinoma. Dis Markers. 2015;2015:616541.

40.

Wang

Zeng

Zhou

, et al. SKA1 promotes malignant phenotype and progression of glioma via multiple signaling pathways. Cancer Cell Int. 2019;19:324.

41.

Wang

Sun

Teng

Zhong

Fan

. Overexpression of spindle and kinetochore-associated protein 1 contributes to the progression of prostate cancer. Tumour Biol. 2017;39(6):1010428317701918.

42.

Ren

Yang

Zhang

Liu

Wang

. SKA2 mediates invasion and metastasis in human breast cancer via EMT. Mol Med Rep. 2019;19(1):515–523.

43.

Cao

, et al. Cyclin-dependent kinase inhibitor 3 (CDKN3) plays a critical role in prostate cancer via regulating cell cycle and DNA replication signaling. Biomed Pharmacother. 2017;96:1109–1118.

44.

Yaqinuddin

Qureshi

Qazi

Farooq

Abbas

. DNMT1 silencing affects locus specific DNA methylation and increases prostate cancer derived PC3 cell invasiveness. J Urol. 2009;182(2):756–761.

45.

Fan

Zhang

, et al. Transcriptional downregulation of miR-127-3p by CTCF promotes prostate cancer bone metastasis by targeting PSMB5. FEBS Lett. 2020;594(3):466–476.

46.

, et al. MCM3AP-AS1/miR-876-5p/WNT5A axis regulates the proliferation of prostate cancer cells. Cancer Cell In. 2020;20:307.

47.

Liu

MCM3AP-AS1 KD Inhibits proliferation, invasion, and migration of PCa Cells via DNMT1/DNMT3 (A/B) methylation-mediated upregulation of NPY1 R. Mol Ther Nucleic Acids. 2020;20:265–278.

48.

Zhang

Luo

Zhang

LncRNA MCM3AP-AS1 regulates epidermal Growth factor receptor and autophagy to promote hepatocellular carcinoma metastasis by Interacting with miR-455. DNA. Cell Biol. 2019;38(8):857–864.

49.

Yang

. YY1-mediated overexpression of long noncoding RNA MCM3AP-AS1 accelerates angiogenesis and progression in lung cancer by targeting miR-340-5p/KPNA4 axis. J Cell Biochem. 2020;121(3):2258–2267.

50.

Liang

Jia

Chen

, et al. LncRNA MCM3AP-AS1 promotes proliferation and invasion through regulating miR-211-5p/SPARC axis in papillary thyroid cancer. Endocrine. 2019;65(2):318–326.

51.

Yang

Sun

Guo

Qin

Liu

. Long non-coding RNA MCM3AP-AS1 promotes growth and migration through modulating FOXK1 by sponging miR-138-5p in pancreatic cancer. Mol Med. 2019;25(1):55.

52.

Chen

Zhu

Feng

. LncRNA MCM3AP-AS1 promotes breast cancer progression via modulating miR-28-5p/CENPF axis. Biomed Pharmacother. 2020;128:110289.

53.

Kobayashi

Komatsu

Ichikawa

, et al. Overexpression of denticleless E3 ubiquitin protein ligase homolog (DTL) is related to poor outcome in gastric carcinoma. Oncotarget. 2015;6(34):36615–36624.

54.

Cui

Wang

Lei

, et al. DTL promotes cancer progression by PDCD4 ubiquitin-dependent degradation. J Exp Clin Cancer Res. 2019;38(1):350.

55.

Guo

Zeng

Huang

Wang

Xie

. MiR-508-3p inhibits cell invasion and epithelial-mesenchymal transition by targeting ZEB1 in triple-negative breast cancer. Eur Rev Med Pharmacol Sci. 2018;22(19):6379–6385.

56.

Zhai

Zhou

Zhao

, et al. Identification of miR-508-3p and miR-509-3p that are associated with cell invasion and migration and involved in the apoptosis of renal cell carcinoma. Biochem Biophys Res Commun. 2012;419(4):621–626.

57.

Wang

Zhao

Yang

, et al. An immune-related lncRNA signature for patients with anaplastic gliomas. J Neurooncol. 2018;136(2):263–271.

58.

Erlandsson

Carlsson

Lundholm

, et al. M2 macrophages and regulatory T cells in lethal prostate cancer. Prostate. 2019;79(4):363–369.

59.

Bianchi-Frias

Damodarasamy

Hernandez

Gil da Costa

. The aged microenvironment influences the tumorigenic potential of malignant prostate epithelial cells. Mol Cancer Res. 2019;17(1):321–331.

60.

Shen

Chen

, et al. In vitro-induced M2 type macrophages induces the resistance of prostate cancer cells to cytotoxic action of NK cells. Exp Cell Res. 2018;364(1):113–123.

61.

Jones

Sinder

Paige

, et al. Trabectedin reduces skeletal prostate cancer tumor size in association with effects on M2 macrophages and efferocytosis. Neoplasia. 2019;21(2):172–184.

62.

Zhang

Zhao

Liu

Jin

. Trichostatin a inhibits proliferation of PC3 prostate cancer cells by disrupting the EGFR pathway. Oncol Lett. 2019;18(1):687–693.

63.

Fortson

Kayarthodi

Fujimura

, et al. Histone deacetylase inhibitors, valproic acid and trichostatin-a induce apoptosis and affect acetylation status of p53 in ERG-positive prostate cancer cells. Int J Oncol. 2011;39(1):111–119.

64.

Makarević

Rutz

Juengel

, et al. Influence of the HDAC inhibitor valproic acid on the growth and proliferation of temsirolimus-resistant prostate cancer cells in vitro. Cancers (Basel). 2019;11(4):566.

65.

Lan

Yuan

, et al. Valproic acid (VPA) inhibits the epithelial-mesenchymal transition in prostate carcinoma via the dual suppression of SMAD4. J Cancer Res Clin Oncol. 2016;142(1):177–185.

66.

Jiang

Zheng

Huang

, et al. Role of SMAD4 in the mechanism of valproic acid’s inhibitory effect on prostate cancer cell invasiveness. Int Urol Nephrol 2014;46(5):941–946.

67.

Wang

, et al. Trichostatin A, a histone deacetylase inhibitor, reverses epithelial-mesenchymal transition in colorectal cancer SW480 and prostate cancer PC3 cells. Biochem Biophys Res Commun. 2015;456(1):320–326.

68.

Kim

Park

, et al. HDAC inhibitor, CG-745, enhances the anti-cancer effect of anti-PD-1 immune checkpoint inhibitor by modulation of the immune microenvironment. J Cancer. 2020;11(14):4059–4072.

69.

Shin

Kim

Lee

Park

Lee

. Trichostatin A sensitizes hepatocellular carcinoma cells to enhanced nk cell-mediated killing by regulating immune-related genes. Cancer Genomics Proteomics. 2017;14(5):349–362.

70.

Lai

, et al. Long non-coding RNA cox-2 prevents immune evasion and metastasis of hepatocellular carcinoma by altering M1/M2 macrophage polarization. J Cell Biochem. 2018;119(3):2951–2963.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.16 MB

0.05 MB

0.04 MB

0.02 MB

0.01 MB