The four OS tissue specimens along with their adjacent non-cancerous tissues used for RNA sequencing were taken from the primary tumor sites of patients diagnosed with OS by Qilu Hospital of Shandong University, and their clinical chart information was shown in Table 1. All participants provided written informed consent, and the study received ethical approval from the review board of the participating institution.

TRIzol (Accurate Biology, Hunan, China) was utilized for RNA isolation and purification of the samples. Subsequently, NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) was employed to perform quality control on the total RNA in terms of quantity and purity. The mRNA containing PolyA (polyadenylate) was specifically captured via Dynabeads Oligo (dT) beads (Thermo Fisher, CA, USA), which underwent two series of purification. To fragment the captured transcript under high-temperature conditions (94°C for 5–7 min), we used the NEBNext^® Magnesium RNA Fragmentation Module (NEB, cat. E6150, USA). Following this step, cDNA synthesis from the fragmented RNA was carried out using Invitrogen SuperScript™ II Reverse Transcriptase (CA, USA). The second strand cDNA synthesis involved DNA polymerase I along with RNase H and DNTP in a buffer solution; terminal repair and poly(A) were also performed. Magnetic beads were then utilized to screen and purify fragments based on their sizes. UDG enzyme digestion followed by PCR assay yielded double-stranded DNAs with a fragment size of 300bp ± 50bp. Finally, LC-Bio Technologies (Hangzhou) Co., Ltd conducted sequencing and data analysis while R package edgeR or DESeq2 analyzed significant differences among the samples.

DEGs were identified by comparing gene expression levels between OS tissues and adjacent non-cancerous tissues based on the criteria of |log fold-change (FC)| >1 .00 and p-value <0 .05. Venn diagrams were generated to determine DEGs that overlapped among the four pairs of specimens (http://bioinformatics.psb.ugent.be/webtools/Venn/). Additionally, heatmaps, histograms, and volcano plots were used to visualize gene expression levels in OS tissues versus adjacent non-cancerous tissues.

The functions of the DEGs were explored using Gene Ontology (GO), a bioinformatics resource that provides information on biological process (BP), cellular component (CC), and molecular function (MF). In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) was utilized for the identification of interactions and response networks involving DEGs.

The Search Tool for the Retrieval of Interacting Genes (STRING) database was utilized to analyze PPIs among the DEGs. STRING provides valuable information on protein interactions, including confidence scores, protein domains, and 3D structures. To identify key molecules involved in the tumorigenesis and development of OS, a PPI network was generated with STRING and was visualized with Cytoscape software, a tool specifically designed for PPI network visualization. We used the MCODE plugin to identify DEGs with high node scores and strong connectivity within the network. Additionally, we employed the Cytohubba plugin to identify hub genes using various topological analysis algorithms.

To elucidate the prognostic performance of the identified hub genes, we conducted survival analysis using the Kaplan-Meier (KM) method. OS samples from the GEO, TCGA, and EGA databases were classified into elevated or reduced expression cohorts according to the median expression of the hub genes. Subsequently, overall survival (OS) and recurrence-free survival (RFS) curves were plotted along with hazard ratio (HR) and 95% confidence interval (CI) calculations to assess associations between hub gene expression profiles and patient outcomes.

Total RNA was obtained from exponentially growing cells with TRIzol (Invitrogen, Carlsbad, CA, USA). The cDNA synthesis was performed with the Evo M-MLV RT Mix Kit with gDNA Clean using RT-qPCR Ver.2 from Accurate Biotechnology (Hunan) Co., Ltd (ChangSha, China). The primer sequences can be found in Table 2. Each sample underwent three replicate runs of RT-qPCR via SYBR Green Premix Pro Taq HS RT-qPCR kit [Accurate Biotechnology (Hunan) Co., Ltd (ChangSha, China)] on a LightCycler 96 instrument (Roche, Basel, Switzerland). GAPDH served as the endogenous control for normalization and the relative gene expression were computed utilizing the 2^−ΔΔCT formula.

Protein samples were prepared from collected cells and subjected to SDS gel electrophoresis before being transferred onto nitrocellulose membranes, which were then blocked at room temperature (RT) for 2 h in 5% skim milk, before overnight incubation with primary antibodies obtained from ImmunoWay Biotechnology Company (Plano, TX, USA). Subsequently, they underwent three washes with Tween 20 in Tris‐buffered saline (TBST) before a 1–2 h incubation at RT in secondary antibodies. After an additional three washes of the membranes, chemiluminescence reagents were used to visualize protein expression levels. The following antibodies were used:anti-β-actin (1:10000; Proteintech, 66009-1-Ig); anti-Vinculin (1:10000; Abways, CY5164); anti-SLC25A4 (1:1000; ABclonal, A15027); anti-CASQ1 (1:1000; ABclonal, A19640); anti-CASQ2 (1:2000; Proteintech, 18422-1-AP).

The cancer genomic atlas (TCGA) database, known for its comprehensive functionalities such as analyzing differential expression, survival rates, correlations, and identifying similar genes, was utilized to assess hub gene profiles in several tumor tissues and their matched healthy tissues based on patient data.

The HPA database offers valuable information regarding the distribution of proteins within different cancer types based on cellular and histopathological characteristics. This resource was employed for bioinformatics data mining purposes to present confocal images showcasing the cellular localization of hub genes.

To explore genetic alterations, mutations, and associated prognosis related to SLC25A4, we utilized the cBioPortal for Cancer Genomics (http://cbioportal.org). This web resource allows for multidimensional analysis and visualization of genomics data alongside clinical features from diverse cancer samples. For our study, we chose the whole genome pan-cancer analysis (ICGC/TCGA, Nature 2020). We obtained information about copy number variations (CNVs) in hub genes found in both cancerous tissues and cell lines.

Transcriptional data from GSE238110 was downloaded and the subjects were divided into elevated and reduced expression cohorts according to the median SLC25A4 contents, and then the inter-cohort differences were assessed by limma. Next, we employed GSEA to identify the DEG-associated signaling networks that are significant between the elevated- and reduced-SLC25A4 expression cohorts. (Subramanian et al., 2005).

The human OS cell lines MG-63, HOS, and Saos-2 were purchased from Genechem Co., Ltd (Shanghai, China), Procell Life Science & Technology Co., Ltd. (Wuhan, China), and BOSTER Biological Technology Co., Ltd (Wuhan, China), respectively. The human osteoblast cell line hFOB1.19 was obtained from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences (Shanghai, China). MEM (Procell, Wuhan, China) was used for HOS and MG-63 cell culture, and McCoy’s 5A medium (BOSTER, Wuhan, China) was used for the culture of Saos-2 cells. Moreover, DMEM/F12 medium (Gibco, MA, USA) was employed for the hFOB1.19 cell culture. All media contained 10% fetal bovine serum (FBS, Gibco) and were grown in an atmosphere containing 5% CO₂. OS cells were maintained at 37°C, while the hFOB1.19 cells were cultured at 34°C. Transient transfection of siRNA (RiboBio, Guangzhou, China), plasmids and the control was conducted with Lipofectamine 2000 (Invitrogen) following the provided directions. Both RT–qPCR and Western blotting were used to verify the transfection efficiency.

CCK-8 assays (Yeasen Biotechnology, Shanghai, China) were used to assess the proliferation of OS cells. HOS and MG-63 cells (1 × 10⁵ per well) were seeded in 96-well plates and incubated for 24, 48, 72, 96, or 120 h. CCK-8 solution (10 μL per well) was then added and incubated at 37°C for 1 h. Absorbances at 450 nm were read in an Epoch microplate spectrophotometer (SpectraMax i3x).

Cell proliferation rates were assessed by the Cell Light™ EdU Apollo567 In Vitro Kit (RiboBio). Briefly, cells cultured in 96-well plates were transfected for 48 h. Next, EdU (50 μM) was added to the cells and incubated (2 h, 37°C). The cells were then fixed in 4% formaldehyde for 0.5 h, followed by incubation with glycine (2 mg/mL) for 5 min and a 10-min cell permeabilization with 0.5% Triton X-100, before washing in PBS, incubation with the Apollo reaction cocktail for 30 min incubation, and washed twice with 0.5% Triton X-100. The cells were then stained with Hoechst 33342 for 30 min and evaluated under fluorescence microscopy (Olympus, Japan).

After transfection, cells were plated into 6-well plates and grown to 100% confluency. Wounds were operated with 200 μL tips in the middle of the wells, rinsed three times using PBS, and then the medium was changed with serum-free medium immediately. Image capture utilized an inverted phase-contrast microscope (OLYMPUS, Japan) at 0 h and 24 h, respectively.

Transwell chambers (8 µm pore size; Corning, NY, USA) were used to examine migration and invasion in OS cells. Cells (1 × 10) were placed in the top chamber which was precoated with 50 µL 1:8 diluted Matrigel (Corning) (to test invasion) or not (to test migration). Cells in the top chamber were in 200 µL of serum-free culture medium whereas cells in the bottom chamber were cultured in medium with 10% FBS. Following 37 °C incubation for 24 h, the cells in the lower chamber were fixed with 4% paraformaldehyde followed by staining with crystal violet for 1 h at RT. Cell numbers in five randomly selected fields were recorded. All the above experiments were repeated three times.

Apoptosis was assessed in treated OS cells after staining with Annexin V-FITC/PI. After cell precipitation and resuspension, the cells were exposed to 5 μL Annexin V-FITC and 10 μL PI Staining solution, and incubated in the dark at RT for 10–15 min. Then, samples were resuspended in 400 μL 1×Binding Buffer, placed on ice, and evaluated by flow cytometry (BD LSRFortessa) within 1 h.

All data were statistically analyzed using GraphPad Prism eight software. Analysis of the differences between the two groups was determined by Student’s t-test. A p-value of <0.05 was considered statistically significant.

In this study, we combined RNA sequencing, bioinformatics, and experimental results for multilevel verification of the identification of key genes in OS. The study design is shown in Figure 1. The RNA-sequencing results of OS and adjacent non-cancerous tissues were analyzed and DEGs were identified using the criteria of |Log2 FC| >1 and p-value <0.05 (Figure 2A). A total of 8353 DEGs were identified between OS and normal tissues, of which 5,056 were upregulated and 3,297 were downregulated (Supplementary Figure S1). There were 4,706 DEGs in OS 1, including 3,348 upregulated and 1358 downregulated genes. In all, 4,173 DEGs were identified in OS 2, of which 2,819 were upregulated and 1354 were downregulated. Additionally, 5,228 DEGs were identified in OS 3, including 3,642 upregulated and 1586 downregulated genes, while 5,662 DEGs were found in OS 4, with 4,207 upregulated and 1455 downregulated (Figure 2B). The top 100 DEGs were visualized using a heatmap, with red color indicating upregulation and blue downregulation (Figure 2C). The Venn diagrams showed that among the top 100 DEGs, there were 81 genes that were common to the four pairs of specimens (Figure 2D).

The biological roles of DEGs in OS tumorigenesis and metastasis were explored using KEGG and GO enrichment analyses. The BP, CC, and MF were analyzed using GO enrichment analyses (Figure 3A; Table 3). In BP analysis, DEGs showed predominant enrichment in muscle contraction, cell adhesion, and muscle filament sliding (Figure 3B). The CC analysis indicated that the DEGs played a role in Z disc, collagen-containing extracellular matrix, and sarcoplasmic reticulum (Figure 3C). Moreover, the MF analysis revealed that the DEGs showed major contribution in protein interaction, extracellular matrix structural component, and actin-binding (Figure 3D). KEGG network enrichment assessment revealed that most DEGs contributed to the calcium, cancer-related, Rap1 and PI3K-Akt networks (Figure 3E; Table 4). The above results indicated that these DEGs were potential regulators of OS initiation and progression.

The interactions among the DEGs were further explored via the STRING database. A PPI of three upregulated and 94 downregulated genes was constructed; this consisted of 97 nodes and 377 edges (Figure 4A). Various algorithms for topological analysis, including MCC, DMNC, MNC, degree, EPC, bottleNeck, EcCentricity, closeness, radiality, and betweenness, were used to predict and identify key nodes in the PPI network using Cytohubba. This led to the identification of 50 top genes (Figure 4B). The interactions between these genes were then evaluated, identifying 22 hub genes, including SLC25A4, ACTN2, TTN, MYL3, MYH7, MYOM2, MYOZ1, UNC45B, SYNPO2L, KLHL31, MYLPF, NRAP, CASQ1, PDLIM3, MYL6B, FLNC, MYH7B, MYO18B, SRL, TXLNB, MYOZ3, and HSPB7 (Figure 4C; Supplementary Figure S2). Hub genes were recognized via the MCODE plugin in Cytoscape, with the most important module containing 16 hub genes (Figure 4D). Among them, seven genes coincided with the hub genes in the Cytohubba plugin (Figure 4E).

K-M survival curves were used to confirm the associations between the hub genes and OS patient survival outcomes (Figure 5; Supplementary Figure S3). The findings showed that increased expression of CASQ1 (HR = 0.63 [0.42–0.94], p = 0.024), CASQ2 (HR = 0.45 [0.3–0.69], p = 0.00017), DES (HR = 0.53 [0.33–0.84], p = 0.0064), PDLIM3 (HR = 0.62 [0.4–0.95], p = 0.026), and SLC25A4 (HR = 0.56 [0.34–0.92], p = 0.02) had better overall survival, while those with higher expression of KLHL31 (HR = 1.68 [1.1–2.58], p = 0.016), MYBPC2 (HR = 1.7 [1.14–2.54], p = 0.008), SYNPO2L (HR = 1.9 [1.12–3.21], p = 0.015), TNNT1 (HR = 1.71 [1.11–2.64], p = 0.015), and MYL3 (HR = 1.5 [1.01–2.23], p = 0.044) had worse overall survival (Figure 5). The abbreviations, full names, and functions of these 10 genes are provided in Supplementary Table S1. In conclusion, these results suggested that the expression of these genes can be indicative of prognosis in OS patients and might act as significant prognostic biomarkers for OS patients.

Table 5 showed the transcript profiles of the ten genes obtained by RNA sequencing. Differences among hub gene transcript expressions in OS were further explored by RT-qPCR. The results showed that DES, CASQ1, CASQ2, MYBPC2, SYNPO2L, TNNT1, SLC25A4, and MYL3 had reduced expression in OS tissues, which consistent with our RNA sequencing results. Whereas PDLIM3 and KLHL31 were significantly upregulated at the expression levels in OS tissues (p < 0.05; Figure 6A). Moreover, we also validated the ten gene expression profiles in OS cell lines. Expression levels of CASQ1, CASQ2, SLC25A4, and TNNT1 were significantly lower in all human OS cells, including MG-63, HOS, and Saos-2, relative to the human fetal osteoblast cell line (hFOB1.19) (Figure 6B). Furthermore, their contents in OS and normal tissues were also demonstrated with TCGA. The expression levels of CASQ1, CASQ2, and SLC25A4 were downregulated in OS, and these results corroborated with the RT-qPCR data (Supplementary Figure S4A). In addition to the expression levels of mRNA, the hub gene protein expressions in hFOB1.19, MG-63, HOS, and Saos-2 cells were further verified using Western blot. We revealed that the SLC25A4 protein expression in MG-63, HOS, and Saos-2 were lower relative to hFOB1.19, while the CASQ1 and CASQ2 protein contents did not decrease significantly (Figure 6C; Supplementary Figure S4B). Then, we validated the expression levels of SLC25A4 in four paired OS samples with adjacent normal tissues. Based on these findings, SLC25A4 was strongly diminished among the tumor group (Figure 6D).

SLC25A4 expression in pan-cancer was then explored, which indicated a decrease in the SLC25A4 expression in most tumors (Figure 7A). As the functions of proteins are often determined by their subcellular localization, the subcellular localization of SLC25A4 in OS was explored using the HPA database. The confocal images revealed that SLC25A4 was present in the mitochondria (Figure 7B). Whole genome pan-cancer analysis (ICGC/TCGA, Nature 2020) was performed to identify the genetic alterations and mutations in SLC25A4 via the cBioPortal database. The frequencies of SLC25A4 mutations in pan-cancer were 3%, and they were typically missense mutation, splice mutation, truncating mutation, amplification, and deep deletion (Figure 7C). Moreover, the frequencies of SLC25A4 mutations were significantly higher in OS patients, accounting for 5.71% (Figure 7D). The data also suggested that the overall SLC25A4-mutated patient survival was worse relative to those with no mutations, indicating that the genetic alterations in SLC25A4 genes was significantly associated with the prognosis of OS patients (Figure 7E). To further elucidate SLC25A4 significance in OS, we employed GSEA to screen the leading 10 networks that showed enrichment in the elevated- and reduced-SLC25A4 expression cohort (Figure 7F). GSEA showed that the high expression of SLC25A4 may potentially mediate certain inflammation, apoptosis and immune-related pathways, such as Oxidative Phosphorylation, TNFα Signaling Via NFκB, IL6/JAK2/STAT3 Signaling. On the other hand, the low expression of SLC25A4 was enriched in the G2M checkpoint, suggesting a possible association with cell proliferation. Therefore, GSEA revealed that SLC25A4 was associated with cell cycle, apoptosis, immune response, cell proliferation and inflammation. In order to further explore the functionality of this molecule, we conducted molecular basic experiments.

In order to investigate the effects of SLC25A4 knockdown on OS, the HOS and MG-63 cells were selected to establish the SLC25A4 knockdown cells. Small interfering (si) RNAs, targeting SLC25A4 and negative control siRNA (siNC) synthesized by RiboBio were transfected into HOS and MG-63 cells. RT-qPCR and Western blot assessments were performed to verify the effects of SLC25A4 knockdown. Based on our results, siSLC25A4-3 exhibited the optimal knockdown effects (Figures 8A,B). Therefore, siSLC25A4-3 was selected for subsequent experiments to validate the effects of SLC25A4 knockdown on OS.

First, we explored the SLC25A4 knockdown-mediated effect on OS cell proliferation. CCK-8 assay was used to evaluate OS cell proliferation, which showed that the SLC25A4 knockdown cell growth rate, including both HOS and MG-63 cells, was significantly increased (Figure 8C). Furthermore, consistent with the CCK-8 data, EdU assessment showed that the downregulation of SLC25A4 expression significantly increased the proliferation rates of HOS and MG-63 cells (Figure 8D). Together, these findings indicated that SLC25A4 deficiency could enhance the proliferative ability of HOS and MG-63 cells.

Given the importance of migratory and invasive features of tumor cells in tumor development and metastasis, wound healing and Transwell assessments were next conducted on SLC25A4 knockdown cells. The wound healing assay revealed that SLC25A4 knockdown enhanced the HOS and MG-63 cell migration (Figure 8E). Furthermore, the Transwell assay, which can be used to explore both cell migration and cell invasion, revealed that SLC25A4 knockdown predominantly promoted the migration and invasion abilities of HOS and MG-63 cells, which were related to the metastasis of OS (Figure 8F).

ANT1, as a component of mPTP, was closely associated with apoptosis; therefore, its effect on OS cell apoptosis was assessed via flow cytometry. Knockdown of SLC25A4 decreased the overall apoptosis rates of HOS and MG-63 cells (Figure 8G).

Since SLC25A4 knockdown played an inhibitory role in tumorigenesis and progression of OS, we further investigated whether SLC25A4 overexpression might have the opposite effect. The SLC25A4-overexpressing plasmid was used to upregulate SLC25A4 expression in HOS and MG-63 cells. SLC25A4 overexpression efficiency was verified using RT-qPCR and Western blot (Figures 9A,B).

SLC25A4 expression-mediated regulation of OS cell proliferation, migration and invasion was investigated by constructing SLC25A4-overexpressing HOS and MG-63 cells. After the overexpression of SLC25A4, HOS and MG-63 cells showed a dominant inhibition of proliferation ability, which was confirmed using CCK-8 and EdU assays (Figures 9C,D). From the perspective of cell migration, wound healing assay confirmed that elevated SLC25A4 expression levels impaired the migration ability of OS cells (Figure 9E). Additionally, the Transwell assay showed that SLC25A4 overexpression strongly suppressed OS cell migratory and invasive abilities (Figure 9F); moreover, flow cytometry revealed that the overexpression of SLC25A4 significantly accelerated HOS and MG-63 cell apoptosis, confirming the promise of SLC25A4 as a potential therapeutic target for OS (Figure 9G).