Unveiling Cellular Traits of Osteosarcoma Cancer Stem Cells via Bulk and Single-Cell RNA-Sequencing Analysis

Zhao, Zhen; Niu, Shun; Yu, Zhixiang; Zhang, Hongtao; Liang, Lizhuo; Xu, Kui; Dong, Chuan; Long, Hua

doi:https://doi.org/10.1155/2023/2923728

Stem Cells International

On this page

Abstract Introduction Materials and Methods Results Discussion Conclusion Data Availability Additional Points Conflicts of Interest Authors’ Contributions Acknowledgments Supplementary Materials References Copyright Related Articles

Special Issue

Stem Cell Therapy and its Application in Precision Medicine

View this Special Issue

Research Article | Open Access

Volume 2023 | Article ID 2923728 | https://doi.org/10.1155/2023/2923728

Unveiling Cellular Traits of Osteosarcoma Cancer Stem Cells via Bulk and Single-Cell RNA-Sequencing Analysis

Zhen Zhao,¹Shun Niu,¹Zhixiang Yu,^2,3Hongtao Zhang,¹Lizhuo Liang,¹Kui Xu,¹Chuan Dong,¹and Hua Long¹

Academic Editor: Katia Mareschi

Received22 Jul 2022

Revised15 Aug 2023

Accepted07 Oct 2023

Published22 Nov 2023

Abstract

Osteosarcoma (OS) is a prototypical sarcoma, predominantly affecting adolescents. The hallmark of cancer stem cell (CSC) behavior pervades OS, invariably signifying an unfavorable prognosis. However, the intricacies underlying OS metastasis remain incompletely comprehended. As the frontiers of the scientific investigation push forward, encompassing both bulk sequencing and single-cell sequencing (scRNA), the domain of bioinformatics finds multifaceted utility. In this study, we combined scRNA and bulk sequencing with the clinical metadata to excavate the latent molecular substrates governing metastatic propensities within OS. Our scRNA analysis indicated that cell-stemness-related pathways might play vital roles in OS metastasis. Subsequently, an autonomous reservoir of bulk sequencing data set was subjected to weighted gene co-expression network analysis (WGCNA) to identify 10 gene clusters. After analyzing the clinical data, we were able to identify two hub genes, HDAC2 and HSPA4, which are strongly linked to the cancer cell stemness. Moreover, we performed transwell assays to validate the regulation of metastatic behaviors by miR-1-3p, HDAC2, and HSPA4 in OS cells. Significantly, our study was fortified by immunohistochemistry (IHC) analyses performed on tumor tissues acquired from OS patients, thereby accentuating the clinical import of our experimental endeavor. Notably, we unveiled the suppressive influence of miR-1-3p on both HDAC2 and HSPA4, with miRTarBase substantiating lower expression in tumor tissue relative to the normative cohort. Insights gleaned from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) program (https://ocg.cancer.gov/programs/target) further augmented the clinical significance of the identified hub genes. In conclusion, our findings highlight the significant role of cell junctions in governing OS stemness and metastasis, underscored by the integration of single-cell sequencing and bulk-sequencing analyses. Moreover, our foundational experiments identified three hub molecules closely associated with the metastatic behaviors. Our findings provide novel insights into the clinical treatment and fundamental understanding of OS.

1. Introduction

Osteosarcoma (OS) represents a prevalent form of malignant bone-derived tumor, often impacting adolescents [1]. Analogous to the other sarcomatous cells, OS cells comprise a cohort of malignant mesenchymal cells. As postulated by Raimondi et al. [2], the bone microenvironment plays a pivotal role in fostering OS progression. Concurrently, OS orchestrates the fabrication of an osteoid matrix and fibrillary stroma, thereby promoting metastasis and invasion. Predominantly, OS manifests at the extremities of long bones, encompassing the femur, tibia, fibula, and humerus. Song et al. [3] have observed a steady annual increase of 0.3% in OS incidence over the past decade. Despite therapeutic interventions such as perioperative neoadjuvant chemotherapy and limb salvage, the 5-years survival rate for individuals with non-metastatic OS stands at 60%–70% [4]. In stark contrast, the 5-years survival rate for those affected by metastatic OS plummets to less than 30% [5], signifying considerably shorter survival durations in the later cohort.

Contemporary investigations underscore metastasis as a salient hallmark of osteosarcoma, integrally linked with cancer stem cells (CSCs). CSCs, a distinct cell type endowed with the capacity to differentiate into diverse lineages, bear pivotal significance in OS regeneration. Extensive evidence points toward the pivotal role of tumor cell stemness in the tumor initiation, progression, and chemoresistance. In recent years, CSCs have assumed the mantle of intense research scrutiny, marked by traits akin to stem cells, including self-renewal, inexhaustible proliferation, and potent differentiation prowess. These attributes enable osteosarcoma stem cells to evade the toxic effects of chemotherapy through self-renewal, unlimited proliferation, and other characteristics during chemotherapy, causing multidrug resistance to chemotherapy, leading to recurrence and metastasis. However, the intricate mechanisms underpinning the sustenance of self-renewal and stemness within tumor stem cells remain enigmatic.

In the realm of OS research, the advent of single-cell RNA sequencing (scRNA) in recent years has provided fresh insights. Zhou et al. [6] unveiled an immunosuppressive microenvironment within advanced osteosarcoma, while Liu et al. [7] illuminated the intricate tapestry of the tumor microenvironment. Currently, harnessing the full potential of microarray technology and bioinformatics analysis offers the promise to discover new clues pertaining to hub genes. However, conventional bulk-sequencing analysis suffers from glaring deficiencies—namely, the dearth of clinical contextualization and the oversight of cellular heterogeneity. These limitations inherently curtail the clinical applicability of bulk sequencing. In this study, we applied a methodology of weighted gene co-expression network analysis (WGCNA) combined with scRNA. WGCNA is an information reduction method and unsupervised classification strategy [8] while scRNA-seq is a powerful method to comprehensively characterize the cellular differentiation within OS tissues [9]. The combination of scRNA and bulk sequencing analyses in delineating gene modules and their in-depth correlation with phenotypic traits has gained widespread traction, offering to uncover potential biomarkers for early cancer diagnosis, thereby fostering a heightened comprehension of hub genes in tumor progression and paving the way for personalized therapeutic targets.

In this study, we procured a multitude of datasets pertaining to OS patients from the Gene Expression Omnibus (GEO) database. These datasets included scRNA data, mRNA expression matrices, and comprehensive clinical profiles. Subsequently, we conducted a meticulous gene co-expression network analysis to identify potential genes linked to OS metastasis. To assess the interconnectivity within this repertoire of promising genes, we harnessed the power of Weighted Gene Co-expression Network Analysis (WGCNA). Based on the clinical information and microarray data from GSE33382, GSE14827, and GSE28423 datasets, we selected the candidate genes (YWHAG, SACM1L, HDAC2, and HSPA4) involved in OS metastasis, and identified two hub genes—HDAC2 and HSPA4. Morisaki et al. [10] reported that HSPA4 was a potential cancer stem cells (CSCs) marker of gastric cancer. Jamaladdin et al. [11] and Wu et al. [12] found that HDAC2 is essential for cellular proliferation and the self-renewal of CSCs by maintaining the expression of the key pluripotency transcription factors. Then we explored the miRNA associated with the key genes we screened. Further, we applied RT-qPCR, western blot, transwell and IHC assays to verify the roles of metastasis in the OS progression of these genes.

2. Materials and Methods

2.1. Data Collection

In this endeavor, we obtained OS patient data from the Gene Expression Omnibus (GEO) database and Therapeutically Applicable Research to Generate Effective Treatments (TARGET) program (https://ocg.cancer.gov/programs/target). This comprehensive data amalgamation included scRNA data, mRNA expression matrices, and clinical sample information. Our data acquisition process rigorously adhered to the ethical and procedural tenets set forth by each dataset.

2.2. Single-Cell Sequencing Data Analysis

Downstream data analysis was performed on publicly available scRNA data using the “Seurat” package (version 4.0.5). First, we filtered out cells having expression data for less than 250 genes. After normalizing gene expression profiles, 2,000 highly variable genes (HVGs) were selected. For identifying marker genes in each sample, we used the “FindAllMarkers” function with adjusted and |log2 fold change| > 0.25 as cutoffs [13].

2.3. Enrichment Analysis

To better understand the potential mechanism of the metastasis in OS, gene set enrichment was performed with the help of KOBAS [14]. The cutoff for the pathway selection was .

2.4. Construction of the Gene Co-Expression Network

We performed gene co-expression network analysis to identify the potential genes related to OS metastasis. Genes and sample profiles were filtered to ensure appropriate data quality, and outlier samples were identified and excluded. WGCNA was performed using the “WGCNA” package in R software (3.6.3) to evaluate the connectivity within these potential genes.

The selection of an optimal soft-thresholding power β for ideal correlation coefficients ensured the construction of a coherent and well-scaled co-expression network. Ultimately, the co-expression patterns were distilled into distinct clusters through cluster analysis. The genes sharing a similar expression pattern were assigned to the same module. Each module contained more than 30 genes. We set 0.25 as a height cutoff to merge several modules for further gene module analysis. We developed more details of the WGCNA in the supplementary materials.

2.5. Identification of Clinically Significant Modules

Gene significance (GS) measures the degree of correlation between a gene and clinical sample information. A higher GS value indicates significant biological relevance, whereas a GS value of zero denotes no specific functional association. Modules with P_GS < 0.05 were considered significant.

2.6. Identification of Hub Genes

Utilizing clinical insights and microarray data extracted from GSE33382, GSE14827, and GSE28423 datasets, we meticulously selected the candidate genes involved in OS metastasis. In accordance with the principles of WGCNA module–trait relationships, they exhibited (1) high within-module connectivity (ModuleMembership > 0.9) and (2) a robust correlation with specific clinical traits (GeneSignificance > 0.3). Bolstering the credibility of our findings, an independent dataset (GSE14827) was subjected to a differential gene expression (DEGs) analysis. Subsequently, the intersection between the genes from the WGCNA module within GSE33382 and DEGs from GSE14827 emerged as the cohort of hub genes. Further, we tried to find some shared attributes in these hub genes at miRTarBase [15]. To explore the hub genes’ clinical value, we performed a series of verification steps within the TARGET dataset using GraphPad Prism 8.

2.7. Stemness Characters Analysis

RNA-sequencing expression profiles and corresponding clinical information for sarcoma were downloaded from the TCGA dataset (https://portal.gdc.com). Use the OCLR algorithm to calculate mRNAsi which constructed by Malta et al. [16]. Based on the mRNA expression signature, the gene expression profile contains 11,774 genes. We used the same Spearman correlation (RNA expression data). The minimum value was subtracted, and the result was divided by the maximum maps the dryness index to the range [17]. For further research, the LinkedOmics (http://www.linkedomics.org/admin.php) was employed to analyze the KEGG pathway of hub genes in sarcoma [18].

2.8. Patients and Cell Line

The study design was approved by the institutional review board and all methods were performed following the relevant guidelines and regulations (ethical approval reference number No.TDLL-202211-01). Written informed consent was obtained before the clinical sample collection. To investigate gene functions of hub genes, we collected metastatic and preinvasive OS samples from 11 patients. A pathologist reviewed the histology of all tumor samples. At the time of sample collection, none of the patients had received therapeutic medications or previous surgical interventions.

The MG63 OS cells, acquired from ATCC, were cultured in high-glucose Dulbecco’s modified Eagle’s Medium (DMEM; Gibco, Carlsbad, CA, USA) with 10% fetal bovine serum (FBS; Gibco). Cells were cultured at 37°C with 5% CO₂.

2.9. Cell Transfection Assay

HDAC2 shRNA, HSPA4 shRNA, miR-1-3p mimics, and their corresponding negative control (NC) were designed and synthesized by GeneChem Co. (Shanghai, China). We followed the manufacturer’s protocol using Lipofectamine 3000 reagent (Invitrogen, CA, USA) for cell transfections.

2.10. RT-qPCR

Total RNA was extracted from transfected cells using an extraction kit (YEASEN Biotech Co. Ltd., Shanghai, China). For detecting miRNA-1-3p expression, 10 ng RNA extract was reverse transcribed with MicroRNA Reverse Transcription Kit (Sangon Biotech Shanghai, China). We then performed RTqPCR using the ABI-Prism 7500 real-time PCR system (Applied Bioscience, USA) using an RT–qPCR kit (Sangon Biotech Shanghai, China). U6 expression was used for normalizing miRNA expression. For detecting HSPA4 and HDAC2, 10 µg of total RNA was used to synthesize first-strand cDNA using a First-Strand cDNA Synthesis Kit (YEASEN Biotech Co. Ltd., Shanghai, China). RT–qPCR was performed using ABI-Prism 7500 real-time PCR system (Applied Bioscience, USA). We used the SYBR green method to perform the qPCR assay (YEASEN Biotech Co. Ltd., Shanghai, China) and used GAPDH as housekeeping control. Expression levels of genes were calculated using the 2^−ΔΔt method.

2.11. Cell Counting Kit-8 (CCK-8) Assay

To evaluate cell viability, a CCK-8 assay was performed. In brief, 0.1 mL culture medium was added and ∼8,000 cells were seeded in each well of a 96-well plate. The cells were then incubated for 2 hr at 37°C with 5% CO₂. A Microplate Reader (OLYMPUS, Japan) was used to measure the absorbances at 450 nm.

2.12. Transwell Migration and Invasion Assay

Transwell migration and invasion assay were applied to evaluate the behaviors of migration and invasion. In the migration assay, transfected cells were resuspended in 200 µL serum-free medium with 2 × 10⁴ cells in the upper chamber and 600 µL medium with 10% serum in the lower chamber. After 24 hr of incubation, cells traversing the membrane to access greater nutrition in the lower chamber were delicately removed. The cells that successfully negotiated the membrane were subjected to staining with crystal violet and subsequently quantified under a microscope.

To evaluate the invasion ability of OS cells, the transwell chambers were coated with Matrigel solution (500 ng/µL; BD, Franklin Lakes, NJ, USA), and the process described above was followed.

2.13. Immunohistochemistry

We performed immunohistochemical staining on 4-µm thick tissue microarray (TMA) slides. After deparaffinization, we applied 1-mM EDTA, pH 9.0 in a steam cooker for antigen retrieval. Subsequently, we applied the HDAC2 and HSPA4 primary antibodies at a concentration of 1 : 100 diluents for 1 hr. Then, the slides were incubated with HRP-conjugated secondary antibody for 15 min. We took diaminobenzidine before mounting, and the slides were counterstained with hematoxylin.

2.14. Statistical Analysis

All experiments were performed in a minimum of triplicate. GraphPad Prism 8 software was used to perform the ROC analysis, KM plot, and Wilcoxon test. Value was considered to indicate a statistically significant difference.

The analysis flow of our research is shown in Figure S1 (Supplementary 1).

3. Results

3.1. Single-Cell Sequencing Analysis of Primary and Metastatic OS Reveals Distinct Metastasis-Associated Changes

Leveraging a primary OS sample (BC16) and a corresponding metastatic OS sample (BC10) from GSE152048, we isolated a total of 2,7691 cells from the two specimens. This comprised 1,7481 cells from BC10 and 1,0210 cells from BC16, post meticulous cell quality control measures. We applied UMAP and TSNE to visualize the cells in cell projection clustering (Figure 1).

(a)

(b)

To elucidate the characteristic changes of metastatic OS, we used the “FingMarkers” function to identify the differentially expressed genes (DEGs) between primary and metastatic OS. There were a total of 1,761 DEGs. KOBAS (http://kobas.cbi.pku.edu.cn/genelist/) [14] analysis of these DEGs indicated significant enrichment of focal adhesion, gap junction, and adherent junction pathways. These findings indicate that the regulation of cell adhesions is a crucial feature of OS metastasis (Table 1).

3.2. Construction of Weight Co-Expression Network

We listed all gene expression variances between 34 metastatic and 19 nonmetastatic OS samples from GSE33382 (Table 2). The Baseline characters of GSE33382 is shown in Supplementary 2. Then we picked the top 25% different genes by analysis on variance (6,244 genes) as DEGs for the next WGCNA analysis. First, the top 25% of DEGs were used to construct the network (Figure 2(a)). We then clustered the samples using the average linkage and Pearson’s correlation methods. We picked the largest sample cluster for further analysis, which included 53 OS samples. We defined the soft-thresholding parameter β = 3 to build the free-scale system (Figure 2(b)), which yielded a high-quality network with an R² value of 0.96 (Figure 2(c)). Average linkage hierarchical clustering identified 10 modules, which were marked in different colors for facilitating visualization. The genes that could not be clustered were assigned to the gray module and were not analyzed further (Figure 3(a)).

(a)

(b)

(c)

Figure 2

(a) Clustering was based on the top 25% variance of the genes in GSE33382. The red color represents the patients’ gender, OS histological subtype, tumor, grade, and metastasis location. Color is related to the different classifications of each clinical phenotype. (b) Determination of the soft-thresholding power (β) with the cutoff at 0.9 in WGCNA and mean connectivity for various soft-thresholding powers. These two figures show that β = 3 may output a high mean number of connections and a scale-free network. (c) Analysis of the scale-free topology model fitting index, R² = 0.96 means the scale-free model is ideal.

(a)

(b)

(c)

(d)

3.3. Identification of Key Modules and Gene Enrichment within the Turquoise Module

Analysis of the heatmap unveiled the relative independence of the module under scrutiny from the other modules (Figure 3(b)). The turquoise module exhibited a significant negative correlation with metastasis (Figure 3(c)), thereby warranting its classification as a key module with the clinical relevance (Figure 3(d)). To further analyze the genes in the modules, we performed enrichment analysis using Metascape (http://metascape.org/). We found that HDAC2 and HSPA4 were enriched in the host cell factor (HCF-1) complex with SAP30 (Table S1 (Supplementary 1)). This indicated that the HCF-1 complex might play a vital role in the metastatic behaviors of OS.

3.4. Identification of Hub Genes and miRNAs

Relying on the aforementioned criteria (|MM| > 0.9, |GS| > 0.3), 157 genes were selected. “However, not all of them exhibited significant differential expression in another independent dataset (Figure 4(a)). The final cohort of hub genes was derived as the intersection between genes from WGCNA modules and DEGs (Figure 4(b)). The rigorous approach yielded four hub genes, namely YWHAG, SACM1L, HDAC2, and HSPA4. To investigate the potential interactions between hub genes and miRNAs, we used miRTarBase and identified HDAC2 and HSPA4 as the potential targets of hsa-miR-1-3p. This indicates a possible inhibitory regulation by hsa-miR-1 in OS metastasis. To validate this, we further selected a miRNA array-sequencing dataset of OS samples (GSE28423) and screened the differentially expressed miRNAs between nonmetastatic and metastatic groups. Interestingly, we found higher expression of hsa-miR-1 in the nonmetastasis group (Figure 4(c)), suggesting its suppressive effect on OS’s metastatic behaviors.

(a)

(b)

(c)

3.5. Hub Genes and miRNA Validation

Within the TARGET dataset, we observed a significant upregulation of HDAC2 and HSPA4 in metastatic groups compared to nonmetastasis groups (Figures 5(a) and 5(b)). The ROC analysis underscored their potential as predictive factors for OS metastasis () (Figures 5(e) and 5(f)). The prognostic analysis through the KaplanMeier Plot highlighted HDAC2’s robust performance (Figure 5(c)), whereas HSPA4’s association with prognosis was nonsignificant (Figure 5(d)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5

(a) HDAC2 expression between metastasis and non-metastasis groups in TARGET database; (b) HSPA4 expression between metastasis and non-metastasis groups in TARGET database; (c) KM-plot for the low HDAC2 expressed and high HDAC2 expressed patients survival with the clinical information from TARGET; (d) KM-plot for the low HSPA2 expressed and high HSPA4 expressed patients survival with the clinical data from TARGET; (e) The ROC for the HDAC2 to be a diagnostic factor for metastatic behaviors of OS; (f) The ROC for the HDAC2 to be a diagnostic factor for metastatic behaviors of OS; means .

Further validation ensued to corroborate miR-1-3p’s downregulatory impact on HDAC2 and HSPA4 expression through in vitro RT–qPCR assays. Initial transfection assays with miR-1-3p mimics in MG63 cells established successful transfection, verified by RT–qPCR (Figure S2 (Supplementary 1)). Subsequent RT–qPCR assays revealed that the upregulation of miR-1-3p significantly downregulated HDAC2 and HSPA4 in MG63 OS cells (, Figures 6(a) and 6(b)).

(a)

(b)

3.6. Stemness Featured of HDAC2 and HSPA4

To depict stemness featured of HDAC2 and HSPA4 in sarcoma, mRNA expression-base stemness index (mRNAsi) and were evaluated, and significant differences were shown for HSPA4 (Supplementary 3). We wished to explore if “HDAC2 and HSPA4 axis” play an important role in stemness. We used the LinkedOmics database and GSEA to analyze which signaling pathways were enriched in sarcoma. High expression of HDAC2 was mainly involved in “cell cycle G2/M phase transition” and “mitotic cell cycle phase transition”. High expression of HSPA4 was mainly involved in the “translational elongation” and “RNA 3’-end procession” (Supplementary 3). “HDAC2 and HSPA4 axis” was significantly enriched in promoting cell cycle which indicated that HDAC2 and HSPA4 play key roles in stemness.

3.7. Inhibition of HDAC2\HSPA4 and Promotion of miR-1-3p Suppresses Metastasis in MG63 Cells

To further investigate the potential regulatory roles of the hub genes HDAC2 and HSPA4 on OS cell metastasis, we conducted transwell assays on MG63 cells transfected with HSPA shRNA and HDAC2 shRNA, respectively (verified via PCR result in Figure S2 (Supplementary 1). As shown in Figures 7 and 8, MG63 cells in the HDAC2-KD, HSPA4-KD, and miR-1-3p-UP groups displayed lower invasion and migration abilities than their counterparts.

These results indicated that HDAC2 shRNA or HSPA4 shRNA, as well as miR-1-3p mimic, showed a decrease in migration and invasion than the negative control groups.

3.8. Inhibition of HDAC2 and Promotion of miR-1-3p Would Accelerate MG63 Cell Death

We performed a CCK-8 assay to evaluate the effect of HDAC2 or HSPA4 and miR-1-3p on cell growth (Figure 9). We found that MG63 transfected with HDAC2 shRNA or miR-1-3p mimics significantly suppressed MG63 cell proliferation (Figures 9(a) and 9(b)). However, low expression of HSPA4 did not inhibit MG63 survival (Figure 9(c)).

(a)

(b)

(c)

3.9. Metastatic OS Cells Express Higher Levels of HDAC2 and HSPA4

Metastatic and nonmetastatic OS tissue samples were obtained from Tangdu hospital. We observed higher expression of HSPA4 and HDAC2 in patients with metastatic OS (Figure 10(a)–10(d)). IHC analysis confirmed the elevated expression. Besides, the clinical scRNA samples showed differential expression of HSPA4 and HDAC2 between primary and metastatic OS (Figure 11).

(a)

(b)

4. Discussion

Metastasis and cell stemness represent pivotal hallmarks across various malignant tumors, including OS. A marked prognostic distinction between metastatic and nonmetastatic OS has been well-established [19]. Necessitating further exploration of the potential hub genes underlying metastatic OS. However, a comprehensive understanding of key molecules involved in the OS remains elusive [20].

To discover novel biomarkers or treatment targets for metastatic OS, we used the mRNA expression matrices and clinical sample information from the GSE33824 dataset. We performed WGCNA and, in conjunction with the clinical data, identified four hub genes—YWHAG, SACM1L, HDAC2, and HSPA4. HDAC2 and HSPA4 have demonstrated significant association with the cancer progression. Subsequently, we subjected these hub genes to experimental validation within the TARGET dataset, thereby affirming HDAC2 and HSPA4 as pivotal actors in the OS metastatic pathways. Leveraging miRTarBase [15] (http://mirtarbase.mbc.nctu.edu.tw/index.html), we identified HDAC2 and HSPA4 as potential targets of miR-1-3p. Our validation efforts encompassed transwell assays, IHC, and RT–qPCR, collectively showcasing miR-1-3p’s role as an OS metastasis suppressor.

Previous studies by Novello et al. [21] underscored miR-1’s ability to suppress OS cell proliferation and regulate the cell cycle. While a significantly lower expression of miR-1 was observed in OS tissues compared to healthy tissues, its expression did not differ between high- and low-grade OS. Moreover, Niu et al. [22] demonstrated that miR-1 inhibits OS cells’ growth, migration, and invasion by downregulating VEGFA. Moreover, using two OS cell lines, namely SAOS-2 and U2OS, they found that downregulation of VEGFA by miR-1 led to the inhibition of cell growth, migration, and invasion in U2OS cells. Moreover, Niu et al. [22] reported that miR-1 could activate p53 independently through PAX3. Together, these findings strongly suggest an inhibitory role of miR-1 in OS.

Consistent with this, our study shows that miR-1-3p could have a suppressive effect on OS cell viability, migration, and invasion via downregulation of HDAC2 and HSPA4, which could be potential therapeutic targets for OS intervention in the future. Several researchers [23–25] have reported that in Duchenne muscular dystrophy, miR-1 is involved in muscle homeostasis in the Dys-nNOS-HDAC2 pathway. Han et al. [26] discovered and explored the miR-1/HSPA4 axis in sepsis mice models. Our study indicated that miR-1-3p was involved in OS metastasis via the regulation of HDAC2 or HSPA4 and validated the regulatory interaction between miR-1-3p and HDAC2/HSPA4.

Histone deacetylases (HDACs) participate in histone acetyl group removal, inducing chromatin compaction [27]. Watanabe et al. [28] reported that HDACs could accelerate OS cell growth. It has been shown that HDAC inhibitors could prevent OS metastasis in vitro [29]. Of significance, HDAC2 stands as a pivotal HDAC family member. Li et al. [30] reported that HDAC2 triggered the migration of OS cells through the upregulation of IL-6. Our findings furnish strong evidence for HDAC2’s ability to facilitate the metastasis of OS cells. Moreover, the observed link between HDAC2 and miR-1-3p offers a potential avenue for OS treatment, unearthing a potential mechanism for miR-1-3p’s anti-metastatic function.

Banerjee et al. [31] highlighted the prominent role of heat shock 70-kDa protein 4 (HSPA4) in triple-negative breast cancer metastasis through network analysis. Additionally, Gu et al. [32] demonstrated that targeting HSPA4 with IgG could help tumor-educated B cells selectively promote breast cancer lymph node metastasis. However, only a few studies have explored the role of HSPA4 in OS progression. In our research, we proved that HSPA4 was associated with the metastasis feature of OS cells. Our transwell assays revealed that HSPA4 could facilitate MG63 cells’ migration and invasion. However, the inhibition of HSPA4 did not affect MG63 cell growth and OS patients’ prognosis, which suggest the presence of an alternative metastatic mechanism in an HSPA4-independent manner. Furthermore, both HSPA4 and HDAC2 emerged as potential targets of miR-1-3p. Subsequent investigations are essential to elucidate the intricate roles played by HSPA4 in OS metastasis.

5. Conclusion

In this study, first we initiated a comprehensive analysis by integrating scRNA and WGCNA. Our investigation led us to identify four hub genes—YWHAG, SACM1L, HDAC2, and HSPA4—and subsequently validate their significance through fundamental experimental studies of metastatic OS and cell stemness. Moreover, we explored HDAC2 and HPSA4 in an independent dataset (TARGET) and verified these two genes’ expression in our samples for the clinical value. The genes proved valuable in metastatic OS and cell stemness research and deserve more research in the future. Intriguingly, our enrichment analysis hinted at the possible implication of the HCF-1 complex in the metastatic behaviors of OS.

Data Availability

Public data can be got following the instructions mentioned in methods part.

Additional Points

Limitation. While our study successfully demonstrated the downregulation of HDAC2/HSPA4 by miR-1-3p, the underlying regulatory mechanisms of miR-1-3p remain unclarified.

Conflicts of Interest

The author(s) declare(s) that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

This research was a collaborative effort by Zhen Zhao and Chuan Dong. Zhen Zhao performed all WGCNA analyses, while both Zhen Zhao and Lizhuo Liang contributed to certain aspects of the fundamental experiment studies. Clinical information was graciously provided by Kui Xu and Hongtao Zhang and Zhixiang Yu. Professors Shun Niu and Hua Long provided guidance and collected clinical samples. Zhen Zhao, Shun Niu, and Zhixiang Yu contributed equally to this work.

Acknowledgments

We thank Bullet Edits Limited for the linguistic editing and proofreading of the manuscript. This work was supported by the Tang Du Hospital Discipline Innovation Development Program (2021LCYJ047), First Affiliated Hospital of AFMU Neck and Back Injury Special Project (FXRYJYT04), and Shaanxi Province Innovation Capability Support Project (2021TD-45).

Supplementary Materials

Supplementary 1. Table S1: part of the enrichment results on Metascape (http://metascape.org/). Figure S1: the flow of analysis of this research. Figure S2: (a) the RT–qPCR results for the HDAC2 expression level in MG63 with HDAC2 shRNA transfection (b) the RT–qPCR results for the HSPA4 expression level in MG63 with HSPA4 shRNA transfection.

Supplementary 2. Baseline characters of GSE3338.

Supplementary 3. The distribution of OCLR scores in HDAC2 and HSPA4 and analysis of the signaling pathway for HDAC2 and HSPA4 in human sarcoma.

References

C. Hiller, J. L. Wegler, and C. P. Forest, “Osteosarcoma: accurately diagnosing this bone-chilling disease,” JAAPA, vol. 29, no. 12, pp. 29–35, 2016.
View at: Publisher Site | Google Scholar
L. Raimondi, A. De Luca, A. Gallo et al., “Osteosarcoma cell-derived exosomes affect tumor microenvironment by specific packaging of microRNAs,” Carcinogenesis, vol. 41, no. 5, pp. 666–677, 2020.
View at: Publisher Site | Google Scholar
K. Song, J. Song, K. Lin et al., “Survival analysis of patients with metastatic osteosarcoma: a surveillance, epidemiology, and end results population-based study,” International Orthopaedics, vol. 43, no. 8, pp. 1983–1991, 2019.
View at: Publisher Site | Google Scholar
W. Yi, W. Yujiao, and H. Weimin, “Multiple recurrence of primary orbital synovial sarcoma: report of two cases and literature review,” Oncologie, vol. 24, no. 4, pp. 927–935, 2022.
View at: Publisher Site | Google Scholar
A. Misaghi, A. Goldin, M. Awad, and A. A. Kulidjian, “Osteosarcoma: a comprehensive review,” SICOT-J, vol. 4, p. 12, 2018.
View at: Publisher Site | Google Scholar
Y. Zhou, D. Yang, Q. Yang et al., “Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma,” Nature Communications, vol. 11, no. 1, p. 6322, 2020.
View at: Publisher Site | Google Scholar
Y. Liu, W. Feng, Y. Dai et al., “Single-cell transcriptomics reveals the complexity of the tumor microenvironment of treatment-naive osteosarcoma,” Frontiers in Oncology, vol. 11, p. 709210, 2021.
View at: Publisher Site | Google Scholar
M. Niemira, F. Collin, A. Szalkowska et al., “Molecular signature of subtypes of non-small-cell lung cancer by large-scale transcriptional profiling: identification of key modules and genes by weighted gene co-expression network analysis (WGCNA),” Cancers, vol. 12, no. 1, p. 37, 2020.
View at: Publisher Site | Google Scholar
X. Li, L. Liu, G. J. Goodall et al., “A novel single-cell based method for breast cancer prognosis,” PLOS Computational Biology, vol. 16, no. 8, Article ID e1008133, 2020.
View at: Publisher Site | Google Scholar
T. Morisaki, M. Yashiro, A. Kakehashi et al., “Comparative proteomics analysis of gastric cancer stem cells,” PLoS ONE, vol. 9, no. 11, Article ID e110736, 2014.
View at: Publisher Site | Google Scholar
S. Jamaladdin, R. D. Kelly, L. O’Regan et al., “Histone deacetylase (HDAC) 1 and 2 are essential for accurate cell division and the pluripotency of embryonic stem cells,” in Proceedings of the National Academy of Sciences of the United States of America, vol. 111 of 27, pp. 9840–9845, 2014.
View at: Publisher Site | Google Scholar
J. Wu, P. Zhu, T. Lu et al., “The long non-coding RNA LncHDAC2 drives the self-renewal of liver cancer stem cells via activation of Hedgehog signaling,” Journal of Hepatology, vol. 70, no. 5, pp. 918–929, 2019.
View at: Publisher Site | Google Scholar
T. Lina, H. Haiyan, Z. Yan et al., “Expression and clinical significance of ACTA2 in osteosarcoma tissue,” Oncologie, vol. 24, no. 4, pp. 913–925, 2022.
View at: Publisher Site | Google Scholar
D. Bu, H. Luo, P. Huo et al., “KOBAS-i: intelligent prioritization and exploratory visualization of biological functions for gene enrichment analysis,” Nucleic Acids Research, vol. 49, no. W1, pp. W317–W325, 2021.
View at: Publisher Site | Google Scholar
H.-Y. Huang, Y.-C. Lin, J. Li et al., “HD Huang. miRTarBase 2020: updates to the experimentally validated microRNA-target interaction database,” Nucleic Acids Research, vol. 48, no. D1, pp. D148–D154, 2020.
View at: Publisher Site | Google Scholar
T. M. Malta, A. Sokolov, A. J. Gentles et al., “Machine learning identifies stemness features associated with oncogenic dedifferentiation,” Cell, vol. 173, no. 2, pp. 338–354.e15, 2018.
View at: Publisher Site | Google Scholar
H. Lian, Y.-P. Han, Y.-C. Zhang et al., “Integrative analysis of gene expression and DNA methylation through one-class logistic regression machine learning identifies stemness features in medulloblastoma,” Molecular Oncology, vol. 13, no. 10, pp. 2227–2245, 2019.
View at: Publisher Site | Google Scholar
S. V. Vasaikar, P. Straub, J. Wang, and B. Zhang, “LinkedOmics: analyzing multi-omics data within and across 32 cancer types,” Nucleic Acids Research, vol. 46, no. D1, pp. D956–D963, 2018.
View at: Publisher Site | Google Scholar
C. Meazza and P. Scanagatta, “Metastatic osteosarcoma: a challenging multidisciplinary treatment,” Expert Review of Anticancer Therapy, vol. 16, no. 5, pp. 543–556, 2016.
View at: Publisher Site | Google Scholar
W. Dongqing, C. Heying, and H. Yi, “Polarized autologous macrophages (PAM) can be a tumor vaccine,” Oncologie, vol. 24, no. 3, pp. 441–449, 2022.
View at: Publisher Site | Google Scholar
C. Novello, L. Pazzaglia, C. Cingolani et al., “miRNA expression profile in human osteosarcoma: role of miR-1 and miR-133b in proliferation and cell cycle control,” International Journal of Oncology, vol. 42, no. 2, pp. 667–675, 2013.
View at: Publisher Site | Google Scholar
J. Niu, Y. Sun, Q. Guo, D. Niu, and B. Liu, “miR-1 inhibits cell growth, migration, and invasion by targeting VEGFA in osteosarcoma cells,” Disease Markers, vol. 2016, Article ID 7068986, 8 pages, 2016.
View at: Publisher Site | Google Scholar
V. Cazzella, J. Martone, C. Pinnarò et al., “Exon 45 skipping through U1-snRNA antisense molecules recovers the Dys-nNOS pathway and muscle differentiation in human DMD myoblasts,” Molecular Therapy, vol. 20, no. 11, pp. 2134–2142, 2012.
View at: Publisher Site | Google Scholar
D. Cacchiarelli, J. Martone, E. Girardi et al., “MicroRNAs involved in molecular circuitries relevant for the Duchenne muscular dystrophy pathogenesis are controlled by the Dystrophin/nNOS pathway,” Cell Metabolism, vol. 12, no. 4, pp. 341–351, 2010.
View at: Publisher Site | Google Scholar
D. H. Nguyen-Tran, N. C. Hait, H. Sperber et al., “Molecular mechanism of sphingosine-1-phosphate action in duchenne muscular dystrophy,” Disease Models & Mechanisms, vol. 7, no. 1, pp. 41–54, 2014.
View at: Publisher Site | Google Scholar
Y. Han, Y. Cai, X. Lai et al., “lncRNA RMRP prevents mitochondrial dysfunction and cardiomyocyte apoptosis via the miR-1-5p/hsp70 axis in LPS-induced sepsis mice,” Inflammation, vol. 43, no. 2, pp. 605–618, 2020.
View at: Publisher Site | Google Scholar
S. Ropero and M. Esteller, “The role of histone deacetylases (HDACs) in human cancer,” Molecular Oncology, vol. 1, no. 1, pp. 19–25, 2007.
View at: Publisher Site | Google Scholar
K. Watanabe, K. Okamoto, and S. Yonehara, “Sensitization of osteosarcoma cells to death receptor-mediated apoptosis by HDAC inhibitors through downregulation of cellular FLIP,” Cell Death & Differentiation, vol. 12, no. 1, pp. 10–18, 2005.
View at: Publisher Site | Google Scholar
X. Mu, D. Brynien, and K. R. Weiss, “The HDAC inhibitor vorinostat diminishes the in vitro metastatic behavior of osteosarcoma cells,” BioMed Research International, vol. 2015, no. 9, Article ID 290368, 6 pages, 2015.
View at: Publisher Site | Google Scholar
L.-B. Li, D. Y. M. Leung, R. J. Martin, and E. Goleva, “Inhibition of histone deacetylase 2 expression by elevated glucocorticoid receptor β in steroid-resistant asthma,” American Journal of Respiratory and Critical Care Medicine, vol. 182, no. 7, pp. 877–883, 2010.
View at: Publisher Site | Google Scholar
S. Banerjee, S. R. Kalyani Yabalooru, and D. Karunagaran, “Identification of mRNA and non-coding RNA hubs using network analysis in organ tropism regulated triple negative breast cancer metastasis,” Computers in Biology and Medicine, vol. 127, p. 104076, 2020.
View at: Publisher Site | Google Scholar
Y. Gu, Y. Liu, L. Fu et al., “Tumor-educated B cells selectively promote breast cancer lymph node metastasis by HSPA4-targeting IgG,” Nature Medicine, vol. 25, no. 2, pp. 312–322, 2019.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2023 Zhen Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

187

Downloads

196

Citations