Identification of topoisomerase 2A as a novel bone metastasis-related gene in liver hepatocellular carcinoma

TOP2A is associated with bone metastasis of LIHC

The mRNA expression data and clinical data of 48 matched tumoral and peritumoral tissues from LIHC patients, including 24 patients with bone metastasis (BM), were retrieved from the GSE27635 database. Firstly, volcano plot showed the DEGs between 48 tumoral and peritumoral tissues (Figure 1A). Meanwhile, the DEGs between LIHC patients with or without BM and their counterpart tissues were displayed, respectively (Supplementary Figure 1A, 1B). Interestingly, TOP2A, BIRC5, and SMARCA were all significantly up-regulated, while NAT2, PLG, and ESR1 were remarkably down-regulated. Among them, the upregulation of TOP2A was most significant in the 3 groups. Moreover, Figure 1B, 1C revealed that a total of 46 genes had differential expressions between LIHC patients with and without BM. Next, we used Venn diagrams (Figure 1D) to identify overlapping BM-related genes with the above DEGs from LIHC and paracancerous tissues, and DEGs from BM and non-bone metastasis (NBM) LIHC. There were 25 overlapping BM-related genes, including TOP2A. To further explore the function of the 25 overlapping BM-related genes, GO enrichment analyses and KEGG pathway analyses were conducted. As shown in Figure 1E, the overlapping BM-related genes were enriched in prostate bud formation, morphogenesis of embryonic epithelium, and homeostasis of number of cells in terms of biological process. In the category of cellular component, the above genes were associated with RNA polymerase II transcription factor complex, DNA repair complex, and transcription factor complex. Additionally, in the analysis of the molecular function, protein C-terminus binding, cytokine receptor binding, and receptor ligand activity were the most enriched terms mediated by the overlapping genes. KEGG pathways analysis revealed that the above BM-related genes were mainly related to prostate cancer, chronic myeloid leukemia, and cell cycle. Pathway analysis further illustrated the relationship between 25 overlapping BM-related genes and cancer-related pathways, showing that TOP2A was closely associated with the activity of EMT and cell cycle pathway (Supplementary Figure 1C). As TOP2A displayed great significance in expression and correlation with LIHC bone metastasis, we selected it for the next research.

Figure 1. TOP2A is associated with bone metastasis of liver hepatocellular carcinoma. (A) Volcano plot of differentially expressed genes (DEGs) between LIHC and paracancerous tissues from GSE27635 database, red indicated up-regulated genes, and blue indicated down-regulated genes. (B) Volcano plot of DEGs between LIHC with BM and NBM. (C) Heatmap of DEGs between LIHC with and without BM. (D) Venn diagram of DEGs from A and B. (E) GO and KEGG enrichment analysis of the above overlapping genes. P, peritumor; T, tumor; BM, bone metastasis; NBM, non-bone metastasis.

Aberrant expression of TOP2A in LIHC patients

In order to further explore the role of TOP2A in LIHC, we first identified its mRNA and protein levels in LIHC. The results showed that the transcriptional level of TOP2A was significantly elevated in LIHC tissues compared with normal liver tissues or their matched adjacent normal tissues (Figure 2A–2C). Moreover, consistent results appeared in 10 GEO datasets, confirming the high expression of TOP2A in LIHC (Table 1). Then, we further analyzed the protein expression level of TOP2A through the HPA database. Supplementary Figure 2A revealed that TOP2A protein expression level was upregulated in LIHC with high or medium staining, but not detected in normal liver tissues. Moreover, we purchased the LIHC tissue microarray, consisting of 16 normal liver samples and 79 cases of LIHC that included 42 cases in TNM stage II, 30 cases in stage III, and 2 cases in stage IV, while others were unknown. Immunohistochemistry (IHC) analysis showed that the protein expression level of TOP2A increased with the progression of clinical staging, with the highest expression in stage IV, higher expression in stage III, and lower expression in stage II, all of which were higher than those in the normal liver tissues (Figure 2D and Supplementary Figure 2B). Besides, correlation analysis between TOP2A expression and clinical characteristics of LIHC patients from TCGA revealed a remarkable association between high expression of TOP2A and clinical staging, histologic grade, and T stage in LIHC (Supplementary Figure 2C–2E). Taken together, these findings suggest that TOP2A is overexpressed at both mRNA and protein levels in LIHC patients.

Figure 2. The expression of TOP2A in LIHC. (A) The mRNA expression level of TOP2A of LIHC and normal tissues in TCGA database. (B) TOP2A mRNA expression of matched LIHC and normal tissues in TCGA database. (C) The mRNA expression level of TOP2A of LIHC and normal tissues in TCGA and GTEx database. (D) IHC of TOP2A protein expression in LIHC with different stage and normal liver tissues from tissue microarrays. (E) The mRNA expression level of TOP2A of LIHC, LIHC-BM and normal liver tissues in Tianjin Medical University Cancer Institute and Hospital. (F) IHC of TOP2A protein expression in three paired LIHC-BM lesions and counterpart primary LIHC tissues from Tianjin Medical University Cancer Institute and Hospital. Magnification: 200 ×, scale bars: 100 μm. ***P < 0.001.

Next, to further explore the relationship between TOP2A expression and bone metastasis in LIHC, we collected clinical samples from our center, including 34 normal liver tissues, 36 LIHC tissues, and 13 bone metastasis tissues. qRT-PCR assay illustrated that TOP2A displayed the highest expression in LIHC bone metastasis lesions, followed by LIHC tumor tissues, while the expression level was the lowest in normal liver specimens (Figure 2E), indicating that TOP2A was involved in the LIHC bone metastasis. Additionally, we collected three paired LIHC-BM lesions and counterpart primary LIHC tumors from our center for IHC analysis (Figure 2F). The results showed that compared with the corresponding primary tumors, the protein expression level of TOP2A was elevated in LIHC bone metastatic lesions. The above findings indicate that TOP2A expression is closely related to the carcinogenesis and progression of LIHC, including bone metastasis.

Prognostic value of TOP2A in LIHC patients

After exploring the characteristics of TOP2A expression levels, we evaluated the prognostic value of TOP2A mRNA expression in LIHC from the TCGA database using Kaplan–Meier plotter. Our results suggested that higher mRNA levels of TOP2A were associated with poorer OS, DSS, and PFI in LIHC patients (Figure 3A and Supplementary Figure 3A, 3B). In addition, the time-dependent ROC curves of TOP2A were applied to predict OS, DSS, and PFI. The AUC for 1-year survival was the largest, while it decreased with the years, for 3-year survival being the smallest, but all values exceeded 0.55, which identified a better survival prediction ability of TOP2A (Figure 3B and Supplementary Figure 3C, 3D). To further explore whether TOP2A was an independent prognostic factor in LIHC patients, a multivariable Cox proportional hazard model was conducted using the TIMER database. The results confirmed that TOP2A was an independent risk factor for clinical outcomes in LIHC patients (HR 1.233, 95%CI 1.044-1.458, P =0.014) (Table 2). We further validated the expression characteristics and prognostic value of TOP2A in LIHC from the ICGC database. Consistently, TOP2A was overexpressed in LIHC (Figure 3C). LIHC patients in the ICGC cohort were divided into high-risk (n = 120) and low-risk groups (n = 120) based on the medium cutoff point. The distribution of risk score, survival status, and expression of TOP2A of each LIHC patient was shown in Figure 3D. Kaplan-Meier curve showed that patients in the high-risk group had worse OS (Figure 3E, HR 7.171, 95% CI 3.021-17.02, P < 0.0001). Time-dependent ROC curves revealed a better predictive ability of TOP2A in LIHC (Figure 3F). Taken together, TOP2A is highly expressed in LIHC and indicates poor prognosis, which may serve as a potential prognostic biomarker.

Figure 3. The prognostic value of TOP2A in LIHC. (A) The overall survival curve of TOP2A in LIHC. (B) ROC curves of TOP2A to predict the sensitivity and specificity of 1-, 2-, and 3-year overall survival. (C) The mRNA expression level of TOP2A of LIHC and normal tissues in ICGC database. (D) Distribution of risk score, survival status, and TOP2A expression profile of LIHC in ICGC database. (E) Kaplan-Meier curve of overall survival in high-and low-risk groups. (F) ROC curve of the TOP2A risk score model. HR, hazard ratio. ***P < 0.001.

TOP2A expression levels and prognostic analysis in pan-cancer

To explore the role of TOP2A in cancers, we further evaluated TOP2A expression in various tumors and adjacent normal tissues. TCGA and GTEx data illustrated that TOP2A expression was significantly upregulated in most cancer types, such as BLCA, CESC, BRCA, CHOL, COAD, ESCA, GBM, HNSC, KICH, KIRC, KIRP, LAML, LIHC, LGG, LUAD, LUSC, OV, PCCG, PRAD, PEAD, SKCM, STAD, TGCT, THCA, UCEC, and UCS, suggesting that TOP2A might function as an oncogene in diverse tumors (Supplementary Figure 4A, 4B). Additionally, the differential expression levels of TOP2A were observed between various tumors and metastatic tumors compared with their matched normal tissues. As the tumor progressed, the expression level of TOP2A was also increased in LUAD, LIHC, SKCM, PRAD, BRCA, and KIRC tumors (Supplementary Figure 4C–4H). Next, we further investigated the prognostic value of TOP2A in pan-cancer using the TCGA database. The results showed that elevated expression of TOP2A was significantly associated with poorer OS in KIRP, LGG, KIRC, MESO, KICH, ACC, LIHC, PAAD, LUAD, and UVM (HR > 1, P-value <0.05). While TOP2A had no prognostic values in other types of cancers (HR ≤ 1, P-value > 0.05) (Supplementary Figure 5). Collectively, these findings indicate that TOP2A involved in tumor carcinogenesis and progression might be a prognosis biomarker in pan-cancer, especially in LIHC.

DNA methylation alterations of TOP2A and gene mutation associated with TOP2A in LIHC

It is well known that genetic and epigenetic changes play a key role in the occurrence and development of cancer. First of all, DNA methylation analysis of TOP2A in LIHC was explored using the TCGA database. As shown in Figure 4A, the DNA methylation level of the TOP2A gene promoter in tumor tissues was significantly lower than that in their adjacent normal tissues. Spearman correlation analysis further revealed a negative correlation between DNA methylation and TOP2A mRNA expression levels in LIHC (Figure 4B), implying that DNA hypomethylation might be an important reason for the high expression of TOP2A in LIHC. Moreover, survival analysis showed that lower DNA methylation of TOP2A was associated with poorer survival (Figure 4C). These results suggest that TOP2A is regulated by epigenetic DNA methylation, which is associated with tumor development in LIHC.

Figure 4. DNA methylation alterations of TOP2A and gene mutation analysis associated TOP2A in LIHC from TCGA database. (A) The methylation level of TOP2A in LIHC. (B) Spearman correlation between DNA methylation and mRNA expression of TOP2A in LIHC. (C) Progress-free survival (PFS) curve of TOP2A methylation in LIHC. (D) Gene mutations associated with TOP2A in LIHC. (E) TOP2A expression level in different TP53 mutant status in LIHC. (F) Overall survival (OS) curve of TP53 mutation in LIHC. ***P < 0.001.

Next, we found that LIHC patients with high TOP2A expression were mainly accompanied by 6 genes mutation, including TP53, CTNNB1, TTN, MUC16, ALB, and LRP1B. Among them, the mutation frequency of TP53 was the highest in LIHC with higher expression of TOP2A, in which missense mutation accounted for the main mutation type (Figure 4D). Therefore, we further explored the relationship between p53 mutation and TOP2A expression. The results showed that TOP2A mRNA expression levels were higher in TP53-mutated LIHC than those in TP53-wild LIHC (Figure 4E), implying that TOP2A as the key regulator in DNA replication and repair might promote tumor aggression by affecting the genetic integrity of TP53. Further survival analysis suggested that TP53 mutation was associated with a poor prognosis of LIHC (Figure 4F). In summary, TOP2A is up-regulated by gene promoter hypomethylation and associated with TP53 gene mutation, which leads to poorer prognosis in LIHC.

Potential function and mechanism analysis of TOP2A in LIHC

Since TOP2A played an important role in LIHC, we further explored its underlying biological function and molecular mechanism. First, co-expressed genes with TOP2A were identified in LIHC using the TCGA database, and the heatmap showed the top 53 genes positively or negatively related to TOP2A in LIHC (Figure 5A). Then the up-regulated or down-regulated genes were separately performed GO and KEGG analyses, to figure out the molecular mechanisms by which TOP2A regulated oncogenesis. KEGG pathway analysis revealed that the up-regulated genes were significantly associated with herpes simplex virus 1 infection, cell cycle, viral carcinogenesis, spliceosome, human T−cell leukemia virus 1 infection, Hippo-YAP signaling pathway, and adherens junction, etc. (Figure 5B, Left). Intriguingly, the Hippo-YAP pathway has been reported to regulate the balance between osteoblastic bone formation and osteoclast bone resorption, and promote bone metastasis [16, 17]. Additionally, the adherens junction was linked to EMT and metastasis [18]. The above results suggest that TOP2A may be involved in bone metastasis of liver cancer. Furthermore, GO analysis showed that these upregulated genes were primarily associated with signal transduction by p53 class mediator, regulation of DNA metabolic process, nuclear division, histone modification, covalent chromatin modification, cell cycle checkpoint, DNA replication, and so on, which all were related to tumor proliferation and epigenetic regulation (Figure 5B, Right). Meanwhile, KEGG analysis revealed that down-regulated genes were mainly enriched in complement and coagulation cascades, retinol metabolism, metabolism of xenobiotics by cytochrome P450, chemical carcinogenesis, and other metabolic-related pathways (Supplementary Figure 6A). GO enrichment analysis displayed that down-regulated genes were primarily connected with the small molecule catabolic process, carboxylic acid biosynthetic process, fatty acid metabolic process, organic acid biosynthetic process, carboxylic acid biosynthetic process, and alcohol metabolic process (Supplementary Figure 6B). GSEA further identified that gene sets related to metastasis, cell cycle, and DNA repair were enriched in LIHC with high expression of TOP2A, suggesting that TOP2A was positively correlated with metastasis and cell cycle in LIHC (Figure 5C and Supplementary Figure 6C).

Figure 5. Functional enrichment analysis of TOP2A in LIHC. (A) Heatmap of TOP2A co-expressed genes in LIHC using TCGA database. (B) KEGG and GO enrichment analysis for the up-regulated genes in LIHC. (C) Gene Set Enrichment Analysis (GSEA) for TOP2A in LIHC.

More importantly, correlation analysis showed a positive relationship between TOP2A expression and PI3K-AKT-mTOR pathway, EMT markers, and tumor proliferation signature (Figure 6A–6C). Along with that, TOP2A was positively correlated with metastasis-related genes like VEGFA, MMP2, and MMP9 (Figure 6D). Taken together, these results suggest that TOP2A may be involved in the occurrence and progression of LIHC through multiple pathways.

Figure 6. Correlation analysis between TOP2A expression and pathways and markers related to metastasis or proliferation using TCGA database. (A–C) Spearman correlation analysis between TOP2A and PI3K-AKT-mTOR pathway, EMT markers, and tumor proliferation signature. (D) Spearman correlation analysis between the expression level of TOP2A and VEGFA, MMP2, and MMP9.

Immune infiltration and immunotherapy related to TOP2A expression

Recent studies have shown that the immune microenvironment exhibits an important influence on tumor progression and prognosis [19]. Therefore, we performed a comprehensive analysis of the correlation between TOP2A expression and immune cell infiltration and response to immunotherapy in LIHC. Firstly, we studied the mRNA expression of TOP2A in different immune subtypes according to molecular typing of immune subtypes using the TISIDB database [20]. The result showed a significant difference in TOP2A expression across different immune subtypes (C1: wound healing, C2: IFN-γ dominant, C3: inflammatory, C4: lymphocyte deplete, and C6: TGF-β dominant) in LIHC (Figure 7A). The highest proportion of patients was observed in subtype C4 and the highest level of TOP2A expression was in subtype C1 in LIHC, which may suggest that TOP2A played an important role in the immune microenvironment and progress of LIHC. Correlation analysis between TOP2A expression and immune cells showed a positive relationship between TOP2A expression and the infiltration of activated CD4 T cell, central memory CD4 T cell, effector memory CD4 T cell, memory B cell, natural killer T cell, and type 2 T helper cell, while TOP2A was negatively correlated with activated CD8 T cell, CD56 bright natural killer cell, CD56 dim natural killer cell, effector memory CD8 T cell, eosinophil, mast cell, monocyte, neutrophil, and type 1 T helper cell (Figure 7B). The tumor microenvironment analysis revealed that TOP2A expression was negatively correlated with the Stromal score and ESTIMATE score, while there was no statistical significance between TOP2A expression and the Immune score (Figure 7C). Immune checkpoint inhibitor is a new tumor treatment strategy, which can effectively improve the prognosis of several malignant tumors [21]. Next, we further investigated the relationship between TOP2A expression and immune checkpoint molecules. Interestingly, Figure 8A exhibited higher expression levels of 8 immune checkpoint molecules (CD274 (PD-L1), CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT, and SIGLEC15) in LIHC with higher TOP2A expression, compared with lower TOP2A expression group and normal tissue group. Moreover, the radar map demonstrated that TOP2A expression was positively correlated with most immune markers and immune checkpoint molecules, including CD274 and CTLA4 (Figure 8B). Higher TIDE scores in LIHC patients with higher TOP2A expression suggest poorer efficacy of immune checkpoint inhibition therapy (Figure 8C). Intriguingly, the expression of TOP2A was higher in the non-response ICB immunotherapy group (Figure 8D). Collectively, these results indicate that TOP2A expression may be related to tumor immune tolerance in LIHC and plays a crucial role in tumor immunotherapy and prognosis.

Figure 7. Correlation analysis of TOP2A expression with immune infiltration and tumor microenvironment in LIHC. (A) TOP2A mRNA expression in different immune subtypes in LIHC via TISIDB. (B) Radar map of correlation between TOP2A expression and the abundance of immune cells. (C) Correlation between TOP2A expression and ESTIMATE score, Immune score, and Stromal score.

Figure 8. Association analysis of TOP2A expression with immunotherapy in LIHC. (A) 8 immune checkpoint genes associated with TOP2A expression. (B) Radar map of correlation between TOP2A expression and immune marker sets in LIHC. (C) Comparison of TIDE score between the higher (G1) and lower TOP2A expression (G2) groups of LIHC. (D) TOP2A expression in groups with a different ICB immunotherapy response status. **P < 0.01, ***P < 0.001.

Hippo-YAP signal pathway mediates TOP2A-driven LIHC cell growth, migration, and invasion

To evaluate the potential role of TOP2A in LIHC, preliminary in vitro experiments were performed. Firstly, we analyzed gene expression of TOP2A in a normal liver cell line (LO-2) and four liver cancer cell lines (HepG2, Hep3B, HCCLM3, and MHCC97H). The results of qRT-PCR and western blot assays indicated that TOP2A mRNA and protein levels were elevated in liver cancer cell lines compared to that in a normal liver cell line, respectively (Figure 9A and Supplementary Figure 7A), among which HCCLM3 cells with high metastatic characteristics showed the highest expression of TOP2A and Hep3B cells also exhibited higher TOP2A expression level. These findings suggested that the overexpression of TOP2A might play a critical role in tumorigenesis development. Therefore, we silenced endogenous TOP2A expression by using small interfering RNA targeting TOP2A (si-TOP2A) and conducted in vitro functional assays in HCCLM3 and Hep3B. The efficacy of TOP2A knockdown in HCCLM3 cells or Hep3B was tested at mRNA and protein levels (Figure 9B and Supplementary Figure 7B). As we observed that TOP2A may be associated with the Hippo-YAP signaling pathway, which was reported to be related to bone metastasis, we adopted western blot assay and verified that p-YAP protein was significantly increased when inhibiting the expression of TOP2A, and MST1/2 inhibitor XMU-MP-1 decreased the level of p-YAP in the TOP2A silencing group (Figure 9C and Supplementary Figure 7C), suggesting that TOP2A could stimulate the Hippo-YAP signaling. CCK-8 and colony formation assays were utilized to validate the effect of TOP2A knockdown on LIHC cell proliferation, revealing that the cell proliferation ability was significantly suppressed after TOP2A knockdown and could be partially rescued by XMU-MP-1 (Figure 9D, 9E). Moreover, transwell assays indicated that TOP2A knockdown could weaken the ability to migrate and invade in HCCLM3 and Hep3B cells, which could be restored by XMU-MP-1 (Figure 9F). Collectively, we conclude that TOP2A, mediating the Hippo-YAP signal pathway, could promote cancer cell growth and progression in LIHC.

Figure 9. Hippo-YAP signal pathway mediates TOP2A-driven LIHC cell growth, migration, and invasion. (A) The expression of TOP2A was detected by qRT-PCR and Western blot in normal liver cell line and four liver cancer cell lines. (B) The knockdown efficiency of TOP2A siRNA was determined by qRT-PCR and Western blot in HCCLM3 and Hep3B cells. (C) The expression level of p-YAP and YAP was detected by Western blot in the indicated groups of HCCLM3 and Hep3B cells. β-actin was used as a control. (D, E) the effect of TOP2A on cell proliferation was assessed by colony formation and CCK8 assays in HCCLM3 and Hep3B cells. (F) The effect of the TOP2A on cell migration and invasion was tested by transwell assays in HCCLM3 and Hep3B cells. *P < 0.05, ***P < 0.001, and ****P < 0.0001. Data presented as mean ± standard deviation (SD). Experiments were repeated at least three times.

TOP2A enhances bone-specific metastatic potential and tumor-induced osteolysis in LIHC

To further explore the role of TOP2A in LIHC bone metastasis, we used MG63 cells to imitate the bone microenvironment in vitro and detected the process of bone-specific metastasis. Through chemotactic migration and heterogeneous cell-cell adhesion assays, we revealed that knockdown of TOP2A could inhibit the chemotaxis migration ability of HCCLM3 and Hep3B cells toward MG63 cells (Figure 10A), and also heterogeneity adhesion to MG63 cells (Figure 10B), which could be partially rescued by XMU-MP-1. Moreover, since LIHC bone metastasis is typically osteolytic, involving the activation of osteoclasts, we next employed the in vitro osteoclastogenesis assay to confirm the functional role of TOP2A in LIHC bone metastasis. Compared to the control groups, fewer TRAP+ multi-nucleic cells were observed among mouse primary bone marrow monocyte (BMM) cells cultured with conditional medium (CM) from HCCLM3 and Hep3B cells with TOP2A silencing (Figure 10C). Simultaneously, XMU-MP-1 can restore the number of TRAP+ osteoclasts induced by the CM from TOP2A-silencing LIHC cells. qRT-PCR further showed that some osteoclast markers including TRAP, CTSK, and NFACT1 would be decreased in the BMM cells treated with CM from LIHC cells with TOP2A knockdown, while which can be rescued by XMU-MP-1(Figure 10D). These results suggested that TOP2A was involved in promoting bone-specific metastatic potential and tumor-induced osteolysis of LIHC.

Figure 10. TOP2A enhances bone-specific metastatic potential and tumor-induced osteolysis in the bone microenvironment. (A) Chemotactic migration of HCCLM3 and Hep3B cells toward MG63 cells was evaluated by transwell assays. (B) The adhesion ability of HCCLM3 and Hep3B cells to MG63 cells was assessed by heterogeneous cell-cell adhesion assay. (C) TRAP staining of bone marrow monocyte (BMM) cells treated with conditioned medium (CM) from HCCLM3 and Hep3B cells. Scale bar, 200 μM. Quantification of TRAP+ osteoclast number. (D) qRT-PCR verified the expression of TRAP, CTSK, and NFACT1 in the indicated groups of BMM cells. GAPDH as the loading control. (E) 2.5×10⁵ control or TOP2A knockdown HCCLM3 and Hep3B cells were injected directly into the tibias of the mice. Representative bioluminescent images taken on day 28 after inoculation are shown. (F) Representative micro-CT images of tibia osteolytic lesions. (G) H&E and TRAP staining of tibia lesions of mice. Quantification of TRAP+ osteoclasts. Scale bar, 100 μm. *P < 0.05, **P < 0.01. B, bone; M, bone marrow; T, tumor. Data presented as mean ± SD. Experiments were repeated at least three times.

Then, we further explored the effect of TOP2A in vivo. Four weeks after injecting control or TOP2A-knockdown HCCLM3 cells directly into the bone marrow cavity of mice tibia, a small animal bioluminescent imaging device detected that TOP2A depletion reduced tumor expansion in the tibia (Figure 10E). Besides, micro-CT of the tibia showed that mice injected with control cells had more severe osteolysis compared to mice injected with TOP2A-knockdown HCCLM3 cells (Figure 10F). Histological H&E staining confirmed a decrease in tibia tumor burden in mice treated with TOP2A-knockdown HCCLM3 cells. TRAP staining found that TOP2A knockdown decreased the number of TRAP+ osteoclasts at the bone-tumor interface (Figure 10G). Therefore, these findings demonstrate that TOP2A is responsible for osteolytic metastasis and tumor-induced osteolysis in LIHC.

Data source and acquisition

The mRNA expression data and clinical information of 48 LIHC patients with and without BM, were obtained from Gene Expression Omnibus (GEO) database (GSE27635) [40]. The other ten public LIHC microarray datasets were also extracted from the GEO database (http://www.ncbi.nlm.nih.gov/geo), including GSE22058, GSE25097, GSE36376, GSE14520, GSE10143, GSE46444, GSE54236, GSE63898, GSE64041, GSE76427. The mRNA data and corresponding clinical information of pan-cancer were derived from the UCSC Xena database (http://xena.ucsc.edu/) [41]. RNA sequencing data of LIHC were acquired from the TCGA database (https://portal.gdc.cancer.gov/), which included gene expression, genomic mutation, DNA methylation, and clinical data [42]. The RNA expression information of normal liver tissue came from Genotype-Tissue Expression (GTEx) database (https://commonfund.nih.gov/GTEx/) [43]. Transcriptome RNA-seq data and clinical follow-up data of LIHC were retrieved from the ICGC database (https://dcc.icgc.org/) [44]. The cell lines mRNA expression of LIHC was derived from the CCLE dataset (https://portals.broadinstitute.org/ccle/about) [45].

Screening of differentially expressed genes (DEGs) associated with bone metastasis of LIHC

After preprocessing RNA-seq data of LIHC from the GSE27635 dataset, the R package “Limma” was used to screen DEGs between LIHC with BM and without BM, or LIHC and paracancerous normal tissue. |log2 fold change (FC)| >1.5 and P-Value <0.05 were regarded as the thresholds to identify DEGs. The results were presented with volcano plot and heatmap. The intersection of DEGs was obtained by using the “VennDiagram” package. The overlapping BM-related genes were further subjected to functional enrichment analyses.

Analysis of TOP2A expression in pan-cancer and LIHC

The RNA expression levels of TOP2A in pan-cancer or LIHC and corresponding normal tissues were obtained from the TCGA/GTEx/GEO/ICGC database. The R v4.0.3 software and package ggplot2 (v3.3.3) were used for the above analyses and visualization. Correlations between TOP2A expression and clinical characteristics of LIHC were investigated using the TCGA database. The protein expression levels of TOP2A in LIHC and normal liver tissues were evaluated using the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/), which is an open-access resource to search for specific genes/proteins for human proteins [46]. The images were available from v23.proteinatlas.org via https://www.proteinatlas.org/ENSG00000131747-TOP2A/tissue/liver#img and https://www.proteinatlas.org/ENSG00000131747-TOP2A/pathology/liver+cancer#img.

Prognosis analysis

The survival data of pan-cancer were acquired from the TCGA database. The “survival” and “survminer” R packages were used to construct the Kaplan–Meier curves with log-rank test, and the forest maps with P-value, HR, and 95% CI were constructed by the R “forestplot” package. HR > 1 and log-rank P < 0.05 were considered to be a risk factor for prognosis and statistically significant. The multivariable Cox proportional hazard model was conducted using the TIMER database (https://cistrome.shinyapps.io/timer/) to further explore whether TOP2A was an independent factor in LIHC patients. The receiver operating characteristic (ROC) curves for 1-, 2-, and 3-year were generated with R package “timeROC” to evaluate the survival predictive power, and areas under the curve (AUC) were calculated.

Analysis of DNA methylation and genomic alterations

For DNA methylation analysis, the DNA methylation data were downloaded from the TCGA database. The R package “methylKit” was used to assess DNA methylation levels in LIHC and normal liver tissues. Correlation between DNA methylation level and TOP2A RNA expression, and survival of LIHC were further explored with corresponding R package. The somatic mutation data of LIHC from TCGA were visualized using the R package “maftools”. Mutant types included missense mutation, in frame del, nonsense mutation, splice site, in frame ins, frame shift del, frame shift ins, and multi hit. The R package was applied to further explore the relationship between TP53 mutation and TOP2A expression, and survival of LIHC.

Functional enrichment analysis

According to the median value of TOP2A expression, they were divided into TOP2A high and low expression groups in LIHC. |log2 fold change (FC)| >2 and P-Value <0.05 were regarded as the thresholds to identify DEGs. The heatmap constructed by the “pheatmap” R package showed the top 50 co-expressed genes with TOP2A in LIHC. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed with the R package “Cluster Profiler” to analyze the potential function and molecular mechanism of TOP2A in LIHC. The R package “ggplot2” was used to perform gene set enrichment analysis (GSEA) analysis. The spearman correlations between TOP2A expression of LIHC and pathways or biomarkers related to metastasis or proliferation were analyzed with R software GSVA package or TIMER web server (https://cistrome.shinyapps.io/timer/) [47].

Analysis of immune infiltration and immunotherapy

TOP2A gene expression across diverse immune subtypes in LIHC was analyzed through TISIDB (http://cis.hku.hk/TISIDB/), which is a web portal for tumor and immune system interaction [48]. The immune subtypes included C1 (wound healing), C2 (IFN-γ dominant), C3 (inflammatory), C4 (lymphocyte deplete), and C6 (TGF-β dominant) subtypes. The “fmsb” R package was utilized to analyze the relationship between the expression of TOP2A and the abundance of tumor immune cell infiltration. ESTIMATE provides researchers with scores for tumor purity, the level of stromal cells, and the infiltration level of immune cells in tumor tissues based on expression data (https://bioinformatics.mdanderson.org/estimate/index.html), which was used to evaluate the relationships between the expression of TOP2A and ESTIMATE score, immune score, and stromal score. The heatmap was applied to analyze the expression distribution of immune checkpoints gene in tumor tissues and normal tissues using the “ggplot2” R package and “pheatmap” R package. The relationship between TOP2A expression and immune markers set was assessed with CIBERSORT algorithm and visualized with “fmsb” R package. Potential immune checkpoint blockade (ICB) response was predicted with Tumor Immune Dysfunction and Exclusion (TIDE) (http://cis.hku.hk/TISIDB/index.php) and ImmuCellAI (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) algorithm.

Patient samples

The Fresh tissue samples of 34 normal liver, 36 LIHC, and 13 liver cancer bone metastasis lesions were obtained at Tianjin Medical University Cancer Institute and Hospital. These tissues were used to detect the RNA expression level of TOP2A. The tissue microarray was purchased from Xi’an Biotech Co., Ltd (cat.no. D950601, Xi’an, China), consisting of 79 LIHC tissues and 16 liver tissues. Among them, 79 cases of LIHC included 42 cases in TNM stage II, 30 cases in stage III, and 2 cases in stage IV, others were unknown. 3 paired primary LIHC and counterpart LIHC bone metastasis sections were collected from Tianjin Medical University Cancer Institute and Hospital. The tissue microarray and sections were used to detect protein expression levels of TOP2A.

Immunohistochemistry (IHC)

The expression of TOP2A protein in paraffin-embedded tissue from patients was evaluated using IHC. In short, 4-μm tissue slices were sequentially deparaffinized by xylene and rehydrated in gradient ethanol, followed by antigen repair with citrate antigen repair buffer (PH6.0). The sections were then exposed to 3% hydrogen peroxide solution to block endogenous peroxidase activity. After blocking with serum for 30 minutes at room temperature, the slices were incubated with primary antibodies against TOP2A (dilution 1:2000; cat no. ab52934, Abcam, USA) overnight at 4° C. Next day, the sections were incubated with secondary antibodies at room temperature for 1 h, and then stained with 3,3’ -diaminobenzidine (ZSGB bio, Beijing, China). The IHC score of protein expression was calculated by multiplying the intensity score and the percentage score. The staining intensity score was designated as follows: 0 (no staining), 1 (light yellow), 2 (brown) and 3 (brownish red). The percentage score was graded as follows: 0 (no positive), 1 (positive<10%), 2 (10%≤positive<50%), 3 (50%≤positive<80%) and 4 (positive≥80%).

Cell culture, transfection, and treatment

The normal human liver cell line LO2 and four liver cancer cell lines (HepG2, Hep3B, MHCC97H, and HCCLM3) were obtained from Professor Zhang Xiaodong. HEK-293T cell lines were purchased from American Type Culture Collection (ATCC, Manassas, VA, USA). Four liver cancer cell lines and HEK-293T cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Cellmax, Beijing, China) containing 10% fetal bovine serum (FBS; Cellmax, Beijing, China) at 37° C in an environment with 5% CO2. LO2 was maintained in 1640 medium containing 10% FBS. Small interfering RNA (siRNA) and small hairpin RNA (shRNA) targeting TOP2A and corresponding negative control siRNA and shRNA were purchased from Shenggong (Shanghai, China). The sequence of TOP2A siRNA was as follows: CCCAACTTTGATGTGCGTGAA. We used a transfection reagent (Zeta Life, Menlo Park，CA, USA) to transfect siRNA according to the manufacturer’s instructions. After transfection 48 h, cells were harvested for RNA and protein extraction and processed for functional assays. For in vivo knockdown of TOP2A in xenograft experiments, the shRNA sequence was cloned into pLKO.1 lentiviral vector and used to construct a stable knockdown TOP2A HCCLM3 cell line. For rescue experiments, TOP2A silencing cells, including HCCLM3 and Hep3B, were treated with the MST1/2 inhibitor XMU-MP-1 (3 μM, MedChemExpress, Shanghai, China) for 24 h.

qRT-PCR

Total RNA was isolated using TRIzol reagent (Yeasen, Shanghai, China), and cDNA was synthesized using the cDNA Synthesis Kit (TransGen, China). Relative quantification of gene expression level was conducted using FastStart Universal SYBR Green Master Mix (Yeasen, China) with GAPDH mRNA as an endogenous control. The expression values of target genes were calculated using the 2 ^−ΔΔCT method. The sequences of gene-specific primers were as follows: GAPDH (human) forward primer, 5’-CATGTACGTTGCTATCCAGGC-3’; GAPDH (human) reverse primer, 5’-CTCCTTAATGTCACGCACGAT-3’. TOP2A (human) forward primer, 5’- ACCATTGCAGCCTGTAAATGA-3’; TOP2A (human) reverse primer, 5’- GGGCGGAGCAAAATATGTTCC-3’. TRAP (mouse) forward primer, 5’- GCAACATCCCCTGGTATGTG-3’; TRAP (mouse) reverse primer, 5’- GCAAACGGTAGTAAGGGCTG-3’. CTSK (mouse) forward primer, 5’- ATGTGAACCATGCAGTGTTGGTGG-3’; CTSK (mouse) reverse primer, 5’- ATGCCGCAGGCGTTGTTCTTATTC-3’. NFATC1 (mouse) forward primer, 5’- CAGTGTGACCGAAGATACCTGG-3’; NFATC1 (mouse) reverse primer, 5’- TCGAGACTTGATAGGGACCCC-3’. GAPDH (mouse) forward primer, 5’- TGTTTCCTCGTCCCGTAG-3’; GAPDH (mouse) reverse primer, 5’- CAATCTCCACTTTGCCACT -3’.

Western blot

Total protein was isolated from cells using RIPA reagent (Solarbio, Beijing, China) containing 10% PMSF inhibitor. Protein samples were separated using 10% SDS-PAGE for 2 h and transferred into PVDF membranes at 250 mA for 1.5 h. After blocking with 10% non-fat milk for 2 h, primary antibodies against TOP2A (dilution 1:1000; cat no.ab52934, Abcam, USA), YAP (dilution 1:1000; cat no.ab52771, Abcam, USA), p-YAP (dilution 1:1000; cat no. ab76252, Abcam, USA) and β-actin (dilution 1:12000; cat no. AF7018, Affinity, China) were incubated at 4° C overnight. The membranes were then incubated with secondary antibodies after washing three times with TBST. Finally, the ECL reagent was used to visualize the bands, and β-actin was used as a loading control to calculate the relative expression of target proteins.

Cell count kit-8 (CCK-8) assay

A total of 3×10³ HCCLM3 or Hep3B cells were plated into the 96-well plate. After 24 h, 48 h, and 72 h, 10 ul of CCK-8 reagent (Yeasen, Shanghai, China) was added to each well. After incubation for 1 h at 37° C, the absorbance of each well at 450 nm was detected.

Transwell assay

The cell migration and invasion ability were evaluated using the transwell assay. 5 x 10^5 HCCLM3 or 1 x 10^5 Hep3B cells suspended in 200 μL serum-free DMEM were added to each transwell cell culture insert (LABSELECT, China; pore size: 8 μm), whereas 600 μL DMEM containing 20% FBS was placed in the lower chamber in 24-well dishes. After incubation for 24 h at 37° C, migrated or invasive cells through the membrane were fixed with methanol for 10 min and stained with 0.1% crystal violet for 30 min. After staining, cells were calculated and imaged under the microscope.

Colony formation assays

Transfected cells were seeded in 6 well plates at a density of 500 cells/well and then cultured for 14 days. The formed colonies were fixed with methanol and then stained with 0.1% crystal violet.

Chemotactic migration assay

In vitro Chemotactic migration assay of HCCLM3 or Hep3B cells toward a mimic bone microenvironment was performed using transwell cell culture inserts. Firstly, the osteogenic MG63 cells were pre-seeded into the bottom chambers of the 24-well plate. Then, 5 x 10^5 HCCLM3 or 1 x 10^5 Hep3B cells suspended in 200 μL serum-free DMEM were added to the upper chambers when MG63 cells reached about 80% confluent. After 24 h, the migrating cells were fixed with methanol, stained with 0.1% crystal violet, and counted under a microscope.

Heterogeneous cell-cell adhesion assay

1 x 10^5 HCCLM3 or Hep3B cells labeled with Green fluorescent protein (GFP) were inoculated into 24-well plate pre-inoculated with MG63 at a confluence of nearly 100%. After 30 minutes, non-adherent cancer cells were washed away by PBS and counted under a microscope.

Osteoclastogenesis assay in vitro

Primary BMM cells were isolated from femur and tibia of 6-week-old C57BL/6 mice and cultured in basal culture medium (α-MEM supplemented with 10% FBS and 1% penicillin/streptomycin) overnight. The next day, non-adherent cells were plated in 24-well plates at 1×105 per well in a 4:1 mixture of basal culture medium/filtered conditional medium (collected from Hep3B and HCCLM3 cells) supplemented with recombinant murine M-CSF (40 ng/mL, Amizona Scientific, LLC, Birmingham, AL, USA). The day 3, recombinant murine sRANKL (50 ng/mL, Amizona Scientific, LLC, Birmingham, AL, USA) was added and replaced medium every 3 days. The tartrate-resistant acid phosphatase (TRAP) staining was conducted on day 7 according to the manufacturer’s instructions (Amizona Scientific, LLC, Birmingham, AL, USA). TRAP+ multinucleated osteoclasts (multinucleated cells ≥3 nuclei) were quantified as mature osteoclasts and counted under a microscope for analysis. qRT-PCR examined the expression of TRAP, CTSK, and NFACT1 (osteoclast markers) in the above BMM.

Animal experiments

Female BALB/c-nude mice (4-6 weeks old) were purchased from Beijing SPF Biotechnology Co., Ltd (China) and maintained under a specific-pathogen free (SPF) facility. All animal procedures were approved by the Animal Care Committee of the Tianjin Medical University Cancer Institute and Hospital and handled in accordance with the Guide for the Care and Use of Laboratory Animals. For intra-tibia injections, mice were anesthetized using isoflurane inhalation. HCCLM3-luc-sh-NC cells and HCCLM3-luc-sh-TOP2A cells (2.5 × 10⁵ in 25 ul PBS) were injected into the tibia bone marrow cavity of 5 mice in each group respectively. Bioluminescence imaging acquired with the IVIS system (Caliper Life Sciences, Hopkinton, MA, USA) was used to detect cancer cell inoculation and progression. The osteolytic lesions caused by the tumor on the tibia were visualized using micro-CT imaging (Inviscan IRIS, Strasbourg, France). After 4 weeks, the hindlimb bones were excised, fixed in 10% neutral-buffered formalin, decalcified with 10% EDTA for 3 weeks, embedded in paraffin for hematoxylin and subjected to eosin (H&E) and TRAP staining. Mature osteoclast number was identified by multinucleated TRAP+ cells along the tumor–bone interface on TRAP-stained sections and quantified as osteoclast number per millimeter interface.

Statistical analysis

All bioinformatics statistical analyses and R packages were performed using R software version 4.0.3. The experimental data were presented as Mean ± SD. Student’s t-test and one-way ANOVA were used to evaluate the difference between groups. Each assay should be conducted independently at least three times. All statistical analyses of experiments were performed with GraphPad Prism 5.0. Unless otherwise specified, a P-value <0.05 was considered statistically significant.