Identification of a novel glycolysis-related gene signature for predicting the prognosis of osteosarcoma patients

Filtering genes through gene set enrichment analysis (GSEA)

The mRNA expression data was downloaded from the genomic data commons data portal (https://portal.gdc.cancer.gov/). The clinical features and prognostic information of the samples were downloaded from the TARGET database (https://ocg.cancer.gov/programs/target). Hallmark gene sets were downloaded from the Molecular Signatures Database (MSigDB). GSEA was performed to analyze specific gene sets among the hallmark gene sets, which differed significantly between metastatic and non-metastatic OS patients. Seven gene sets, namely, UV response, MYC targets V1, MYC targets V2, hypoxia, E2F targets, glycolysis, and unfolded protein response, were significantly enriched (P < 0.05) (Figure 1 and Table 1). As glycolysis has been clearly shown to be related to metastasis [19], and glycolysis in OS is poorly researched, we chose the glycolysis gene set (P = 0.035), which contained 198 genes, to further analyze how the function of glycolysis-related genes differs between patients with metastatic and non-metastatic OS.

Figure 1. Enrichment plots of seven gene sets which were importantly differentiated between in metastasis and non-metastasis tissues. The horizontal axis represents genes in gene sets ranked by decreasing risk score. The vertical axis represents enrichment scores. The enrichment score increased with the number of enriched genes.

Identification and construction of a three-mRNA signature to predict the prognosis of patients with OS

According to the glycolysis gene set in MSigDB, 198 glycolysis-related genes were selected. First, we compared the expression profiles of these genes between metastatic and non-metastatic samples in the TARGET database. A heatmap was drawn to illustrate the expression profiles (Figure 2A). We found 23 glycolysis-related genes, including TXN, ABCB6, SLC16A3, CAPN5, B3GAT3, ENO1, PGAM1, ANKZF1, P4HA1, and FAM162A, that were altered in metastatic tissues. Then, we performed Kaplan-Meier (KM) analysis and found that 10 of the 23 genes showed a high correlation with the overall survival rates of patients with OS (Supplementary Figure 1).

Figure 2. Identification and construction of a three-mRNA signature to predict the prognosis of patients with OS. (A) Heatmap was constructed showing 23 differentially expressed glycolysis-related genes in metastasis OS tissues compared with non-metastasis tissues. (B) Forest plot of univariate Cox regression analysis of the survival-related 10 differentially expressed genes in OS. (C) Multivariate cox regression analysis identified a risk signature include 3 glycolysis-related genes among the 10 differentially expressed genes. (D) Distribution of risk score in the high-risk group and the low-risk group. (E) Survival status between the high-risk group and the low-risk group. (F) Heatmap of the expression profile of the included glycolysis-related genes.

Subsequently, we performed a univariate Cox regression analysis of these 10 genes to further search for genes associated with patients’ overall survival and prognosis. Six genes were found to correlate with the prognosis of patients (P < 0.05) (Figure 2B). Finally, these mRNAs were performed by multivariate Cox regression analysis in order to establish a risk signature. Accordingly, three genes (ABCB6, P4HA1, and STC2) were identified as “risk genes” (Figure 2C and Table 2). We calculated the risk scores of all samples in the dataset by applying this three-mRNA signature as follows: expression of ABCB6 × 0.946 + expression of P4HA1 × 0.413 + expression of STC2 × 0.435. After ranking the risk scores, we considered the median score as the cut-off and classified the samples into low- and high-risk groups (Figure 2D). Thereafter, the survival rates of patients in the two groups were compared; low death rates and high survival rates were sorted in the low-risk group (Figure 2E). We also drew a heatmap to illustrate the expression profiles of the three genes in both groups (Figure 2F). The expression of all three genes was low in the low-risk group. Thus, our constructed risk signature could be a good predictor for the prognosis of OS patients.

Predicting survival rate and metastasis in OS patients based on the glycolysis-related risk signature

We used Kaplan-Meier analysis to study the differential survival rate between the two groups and observed a significantly higher survival rate was shown in the low-risk group (P < 0.001) (Figure 3A). Next, we performed a receiver operating characteristic (ROC) curve analysis to assess whether the risk signature was efficacious in forecasting outcomes for OS patients. The areas under the curve (AUC) in this risk signature were 0.884 and 0.790 for predicting 1- and 3-year survival, respectively (Figure 3B, 3C). We also used the ROC curve analysis to compare the predictive efficiencies of the three-gene signature and the single genes alone. The results showed that the three-mRNA signature indeed had a higher AUC value for predicting the survival rate of OS patients compared to the individual genes, for both 1 year and 3 years. (Figure 3B, 3C).

Figure 3. The glycolysis-related risk signature could predict survival rate and metastasis of OS patients. (A) Kaplan-Meier survival curve of overall survival rate among OS patients from the low-risk group and the high-risk group. Patients in high-risk group had the poorer prognosis. (B) Survival prediction between risk score and single gene was assessed by time-dependent receiver operating characteristic (ROC) curve for 1 years. (C) Survival prediction between risk score and single gene was assessed by time-dependent receiver operating characteristic (ROC) curve for 3 years. (D) Risk scores among metastasis and non-metastasis groups. (E, F) Kaplan-Meier survival curve of overall survival rate for prognostic value of risk signature for the patients divided by clinical feature of metastasis.

The primary cause of death in OS patients is lung metastasis. We found that the metastatic group in our study showed higher risk scores than the non-metastatic group (P < 0.001) (Figure 3D). We stratified OS patients into two subgroups based on their metastasis status and found that the risk signature was more efficient in predicting the overall survival rates of metastatic patients (P < 0.01) compared to non-metastatic patients (P = 0.091) (Figure 3F), indicating that the risk signature was associated with metastasis in OS.

Additionally, we explored whether the risk signature could represent glycolysis activity. Firstly, we performed GSEA and found that four glycolysis-related gene sets (GO_GLYCOLYTIC_PROCESS, HALLMARK_GLYCOLYSIS, KEGG_GLYCOLYSIS_GLUCONEOGENESIS, and REACTOME_GLYCOLYSIS), were significantly enriched in the high-risk group (Supplementary Figure 2A and Table 3). Secondly, core genes in the glycolytic pathway, including HK2, PGAM1, ENO1, and LDHA, showed an elevated expression in the high-risk group (P < 0.001) (Supplementary Figure 2B). These results suggested that the risk signature could simultaneously reflect the activity of glycolysis and predict the survival and metastasis of patients with OS.

Independence of the three-mRNA signature in predicting the prognosis of patients with osteosarcoma

To verify whether the risk scores were independent of other clinical characteristics, we performed univariate and multivariate Cox regression analyses. The age of 15 was the median age among the 85 OS patients, with 48 male and 37 female patients were included in the TARGET database. This group included 33 patients who presented with metastasis, and 52 without. Through both univariate and multivariate Cox analyses, the risk scores and the clinicopathological parameters of metastasis were demonstrated to be independent prognostic factors in OS to predict the prognosis of OS patients (P < 0.001) (Figure 4A, 4B).

Figure 4. Risk score based on the glycolysis-related risk signature could be an independent prognostic factor. (A) A forest plot of univariate cox regression analysis of risk score and different clinical feature in OS. (B) A forest plot of multivariate cox regression analysis of risk score and different clinical feature in OS.

Validation of the prognostic risk signature

A GEO data set (GSE21257) was chosen by us as a validation cohort to ensure the reliability of the established glycolysis-related risk signature. This cohort comprised 53 cases of which 25 were in the high-risk group and 28 in the low-risk group, based on the risk signature. Low survival rates of patients were consistently seen in the high-risk group (Figure 5A). The risk signature’s prognostic efficiency was assessed using the ROC curve analysis in the validation group; the AUC was 0.740 at 1 year and 0.759 at 3 years (Figure 5B). Moreover, in this validation cohort patients with metastasis had a higher risk score than those without metastasis (P<0.01) (Figure 5C). Higher mortality and lower patient survival rates were observed in the high-risk group compared with the low-risk group (Figure 5D, 5E). A heatmap was drawn according to the expression profiles of the risk genes between both groups (Figure 5F). We found that the risk genes showed a higher expression in the high-risk group compared to the low-risk group. We further verified the prognostic potential of the three glycolysis-related risk genes in the R2 database (http://r2.amc.nl). The high-P4HA1 expression group had a lower overall survival probability and metastasis-free survival probability than the low-P4HA1 expression group (Supplementary Figure 3A). The high-STC2 expression group had a lower overall survival probability than the low-STC2 expression group (Supplementary Figure 3B). The high-ABCB6 expression group also showed a lower metastasis-free survival probability than the low-ABCB6 expression group (Supplementary Figure 3C). We also collected five metastatic human OS tissues and five non-metastatic human OS tissues. The results of western blotting showed that the risk genes (P4HA1, ABCB6, and STC2) had higher expressions in metastatic human OS tissues than non-metastatic human OS tissues (Supplementary Figure 4). These results demonstrated that the risk model constructed from the three risk genes was a good predictive instrument for prognosis in the validation cohort.

Figure 5. Validation of the prognosis risk signature in GSE21257. (A) Kaplan-Meier survival curve of overall survival rate among OS patients from the low-risk group and the high-risk group in GSE21257. Patients in high-risk group had the poorer prognosis. (B) Survival prediction of the risk signature in GSE21257 was assessed by time-dependent receiver operating characteristic (ROC) curve for 1 and 3 years. (C) Risk scores among metastasis and non-metastasis groups in GSE21257. (D) Distribution of risk score in the high-risk group and the low-risk group in GSE21257. (E) Survival status between the high-risk group and the low-risk group in GSE21257. (F) Heatmap of the expression profile of the included glycolysis-related genes in GSE21257.

Changes in the tumor microenvironment of OS patients based on the glycolysis-related risk signature

Recent research has emphasized the pivotal importance of whole organism nourishment and metabolism in the function of immune cells [31], but the relationship between immune response and glycolysis in OS is unknown. To further explore the potential connection between immune response and glycolysis in OS, we calculated the immune score of the patients in the training and validation cohorts, using the ESTIMATE algorithm. In the training cohort, we found that the patients in the low-risk group had immune scores higher than those in the high-risk group (P = 0.059, almost statistically significant) (Figure 6A), and the similar results were seen in the validation cohort (P < 0.001) (Figure 6B).

Figure 6. Changes of tumor microenvironment of osteosarcoma patients based on glycolysis-related risk signature. Immune score was calculated between high-risk and low-risk group in training cohort (A) and validation cohort (B). Kaplan-Meier survival curve of overall survival rate among high-risk and low-risk group in training cohort (C) and validation cohort (D).

After dividing the patients into high- and low-immune score groups according to the median value of the immune scores in the training and validation cohorts, we performed Kaplan-Meier survival analysis of the data. The results showed that a higher immune score was related to a good prognosis in both the training cohort (P < 0.01) (Figure 6C) and validation cohorts (P = 0.062, almost statistically significant) (Figure 6D). In addition, we performed single-sample gene set enrichment analysis (ssGSEA) to calculate the relative abundance of 29 immune-related items of each sample in both cohorts. The 29 immune gene sets included immune cell types and functions, tumor-infiltrating lymphocytes (TILs), proinflammatory, para-inflammation (PI), cytokine and cytokine receptor (CCR), human leukocyte antigen (HLA), regulatory T (Treg) cells, and immune checkpoint. We found that immune infiltration levels were higher in the low-risk group in both the training cohort (Figure 7A) and validation cohorts (Figure 7B). In particular, the “APC_co_inhibition”, “Neutrophils”, “Th2_cells”, “CD8+_T_cells” and “HLA” gene sets showed differential expression between low- and high-risk groups in both cohorts. Meanwhile, we analyzed the correlation among immune cell types, immune pathways, and the risk scores, and found that the proportions of CCR, CD8+_T_cells, neutrophils, TILs, dendritic cells (DCs), HLA, inflammatory promoting cells and immune checkpoints showed negative correlation with the risk scores in both training cohort (Figure 8A) and validation cohort (Figure 8B). We also verified the differential expression of CD8+ T cell between metastatic and non-metastatic human OS tissues in clinical samples. The result suggested that patients with metastasis had lower contents of CD8+ T cell than those without metastasis (Supplementary Figure 5). Thus, we demonstrated that the risk signature correlated with the immune response and immune infiltration levels in OS samples. Furthermore, these results also showed that high immune infiltration levels could suggest good prognosis in OS patients.

Figure 7. Analysis of 29 immune gene sets between high- and low-risk groups of OS. The heatmap was used to visualize the proportions of these gene sets between high- and low-risk groups of OS in training cohort (A) and validation cohort (B).

Figure 8. Analysis of the correlation between risk scores and the immune gene sets. Spearman rank analysis was used to determine the association between risk score and the immune gene sets in training cohort (A) and validation cohort (B).

Validation of the function of P4HA1 in vitro and in vivo

To further examine the role of the three genes in the metastasis of OS, P4HA1 was analyzed in more detail. We verified the expression of P4HA1 in OS cell lines and the hFOB1.19 human osteoblast. The results of western blotting indicated that P4HA1 showed a higher expression in human OS cells than in the osteoblast cell line hFOB1.19 (Figure 9A). To further verify the function of P4HA1 in OS progression, we transfected si-P4HA1 into HOS and 143B OS cell lines. Knockdown efficiency was verified by western blotting (Figure 9B) and si-P4HA1-1 selected for the subsequent experiments. According to the results of the wound-healing assay, knockdown of P4HA1 could restrain the migration of OS in the HOS and 143B cell lines (Figure 9C). Furthermore, the transwell invasion assay results corroborated with those of the wound-healing assay in HOS and 143B cell lines (Figure 9D), indicating that P4HA1 indeed played a vital role in the metastasis of OS in vitro. Next, we further verified the role of P4HA1 in OS metastasis in vivo. The numbers of lung metastatic nodules in P4HA1-knockdown orthotopic OS-implanted mice were fewer than those in the control group (Figure 9E). These results further confirmed our previous results (Figure 2A–2C).

Figure 9. P4HA1 promoted the metastasis of OS in vitro. (A) Western blotting analysis of the expression of P4HA1 in human osteoblastic cell line and osteosarcoma cell lines. (B) Western blotting analysis of the expression of P4HA1 in HOS and 143B cell lines after transfection of si-control, si-P4HA1-1 and si-P4HA1-2. (C) Wound healing assay analysis in HOS and 143B cell lines after transfection of si-control and si-P4HA1-1. Representative images of migration were shown in the right panel. The degrees to which the wounds healed was shown in the histogram. (D) Transwell invasion assay analysis in HOS and 143B cell lines after transfection of si-control and si-P4HA1-1. Representative images of invasion were shown in the right panel. The proportions of invading cells were shown in the histogram. The bars indicate the mean±s.d. Statistically significant differences (t-test), **P<0.01, ***P<0.001. (E) The number of lung metastasis nodules formed by P4HA1-knockdown 143B cells and control cells in orthotopic osteosarcoma implanted mice was shown in the right panel. Statistically significant differences (t-test), **P<0.01. Representative images of lung morphology and lung metastatic nodules formed by P4HA1-knockdown 143B cells and control cells were shown in the left panel. Corresponding HE staining was shown in the middle panel. Scale bars = 200 μm.

Data gathering

We downloaded the mRNA expression data of patients with OS from the genomic data commons data portal (https://portal.gdc.cancer.gov/). Complete clinical information of the above patients were downloaded from the TARGET database (https://ocg.cancer.gov/programs/target). The mRNA expression data and clinical characteristics of OS samples in GSE21257 were downloaded from the NCBI gene expression omnibus (https://www.ncbi.nlm.nih.gov/geo/, GSE21257). All of the clinical information in the two datasets were shown in Supplementary Table 1.

Gene set enrichment analysis (GSEA)

Gene sets which were significantly different between non-metastasis and metastasis groups in OS samples were screened out by performing GSEA (http://www.broadinstitute.org/gsea/index.jsp). The glycolysis-related genes were downloaded from the Molecular Signatures Database (MSigDB) in GSEA.

Establishment and validation of a glycolysis-related risk signature

We performed univariate cox analysis and Kaplan-Meier (KM) survival analysis to assess the role of glycolysis-related genes played in predicting the prognosis of OS patients. Statistically different genes acquired in both univariate variable Cox analysis and KM survival analysis were selected for further analysis in multivariable Cox analysis. Then, we obtained the prognostic risk score formula through the multivariate cox regression analysis according to a linear combination of expression levels weighted with the regression coefficients. The formula was as follows:

$\text{Risk score=}{{{\displaystyle \sum }}^{\text{}}}_{i=1}^n({\text{coef}}_{\text{i}}*{\text{Expr}}_{\text{i}})$

where Expr_i is the expression levels of gene i in a particular patient and coef_i is the coefficients of the gene i in the multivariate cox regression analysis. The median of risk scores were took by us as the cut-off, through which the patients were divided into high- and low-risk groups. Next, the KM survival analysis was performed to assess the differential survival rate between different groups and the ROC curves were used to evaluate the predictive efficiency of this risk signature. At last, we performed univariate and multivariate cox regression analyses to verify whether the risk score calculated by risk signature could be an independent prognostic factors of OS patients in the TARGET dataset. To further verify the accuracy of the risk signature, we chose the GSE21257 dataset as the validation cohort containing 53 patients with OS, which had intact clinical information. All of data in the two cohorts were adjusted in R software via the “sva” package. The TARGET dataset was recognized as the training cohort. According the risk signature, patients in GSE21257 were divided into low-and high-risk groups. Lastly, we performed KM survival analysis and ROC curves to assess the efficiency and accuracy of the risk signature in GSE21257 dataset.

Immune scores in OS samples

Firstly, we used the “limma” package in R software to normalize the matrix data of gene expression [44]. Then, immune score of each sample was computed by ESTIMATE algorithm according to the normalized gene expression [45].

SsGSEA in OS samples

To quantify the proportions of immune signatures in the OS samples based on the ssGSEA score, we used the 29 immune signatures, including cell types, functions, and pathways [46].

Immunohistochemistry (IHC)

IHC assays were conducted as reported previously [47]. Images were obtained with Leica microscope (Leica, DM4000b).

Osteosarcoma cell lines

The OS cell lines (143B, U2OS, Saos2, HOS and SJSA-1) and osteoblast cell line (hFOB1.19) were acquired from ATCC (Manassas, VA, USA). DMEM culture medium mixed with 10% FBS and 1% penicillin/streptomycin was used to culture the cell lines. All cell lines were cultivated in the incubator at 37° C with 5% CO2.

Western blotting

Firstly, we used the RIPA buffer (Beyotime, Shanghai, China) mixed with a protease and a phosphatase inhibitor cocktail (Sigma-Aldrich, USA) in advance to lyse the cells for 30 minutes on ice. After centrifugation, supernatant of lysate was mixed with loading buffer (Beyotime, Shanghai, China). After degeneration of lysates, they were separated by SDS-PAGE. The Bio-Rad (Hercules, CA, USA) was used shift the protein from SDS-PAGE to polyvinylidene difluoride (PVDF) membranes. Then, 5% non-fat milk diluted by TBST was incubated with the PVDF membranes for one hour at room temperature, after which the membranes were incubated with primary antibodies at 4° C overnight. The next day, TBST was used to wash the membranes for three times after which incubated with secondary antibody (Sigma-Aldrich) for 1h at room temperature. Lastly, after washing the membranes for three times by TBST, ECL detection reagents were used to detect the signals. Primary antibodies: P4HA1 (12658-1-AP, Proteintech).

RNA interference

P4HA1 siRNA duplexes and corresponding si-Control were purchased from RiboBio (Guangzhou, China). Lipofectamine 3000 reagent (Invitrogen, Carlsbad, CA, USA) was used to transfect cells based on the manufacturer's protocol. The sequences targeting P4HA1 were as follows: 5’-CCTGTGGTGTCTCGAATTAATdTdT-3’ (si-P4HA1-1); 5’- GGAAUUACAGGUAGCAAAU-3’(si-P4HA1-2).

Transwell invasion assay and wound-healing migration assay

Cell invasion capabilities were detected by using Transwell chambers precoated with the Matrigel matrix as previously described [47]. We seeded 5*10⁴ cells into the upper compartment with 100μl serum-free culture medium. Meanwhile, we added 600μl DMEM with 10% FBS into the lower compartment as a chemoattractant. Then, the plates were incubated for 12h at 37° C, after which we used a cotton bud to clear the non-invaded cells in the upper chamber and used 4% paraformaldehyde to fix the invaded cells on the bottom surface. Fifteen minutes later, 0.1% crystal violet was used to stain the invaded cells. We used an inverted microscope (Olympus) to acquire the images, and the corresponding cells were counted by ImageJ software (National Institutes of Health, MD). Then, we performed the wound-healing migration assay. Six-well plate was used to seed the cells at the density of 80% in advance. Then, we used a sterile 100 μL pipette tip to create a “wound” when the cells reached full confluence. Afterwards, we observed the breadth of “wound” at different time points with a fluorescence microscope (Olympus). The distances were measured by ImageJ software (National Institutes of Health, MD).

Cell transfection

The P4HA1 lentiviral shRNA plasmid was supplied by Genomeditech (Shanghai, China). The sequences targeting P4HA1 was as follows: shP4HA1, CCGGCCTGTGGTGTCTCGAATTAATCTCGAGATTAATTCGAGACACCACAGGTTTTTG. The lentivirus was attained by transferring the plasmid into HEK-293T cell lines with Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA) in according with the manufacturer’s instructions.

In vivo animal experiments

The animal experiments were approved by the Animal Research Committee of the Shanghai General Hospital (Ethical code: 2019AW004) and were performed in accordance with established guidelines. For the growth and metastasis assays in vivo, four-week-old male nude mice were orthotopically inoculated into the right tibia with a micro syringe. A total of 1 × 10⁶ corresponding cells suspended in 20 μL of PBS were injected into each nude mouse. One month later, the mice were sacrificed, the tumors and individual lung tissues were excised and fixed with 4% paraformaldehyde. Metastatic lung tissues were analyzed by H&E staining.

Statistical analysis

Statistical analyses were performed mostly based on the R 3.5.1 (https://www.r-project.org/) and GraphPad Prism 5 (GraphPad Software Inc, La Jolla, CA). The univariate cox regression analysis was used to assess the correlation between gene expressions and prognosis of patients. The multivariate cox regression analysis was used to establish the risk signature based on the genes related to prognosis of patients with OS. Differential analysis of survival rates between high- and low-risk groups generated by the KM analysis was defined by log-rank tests. The package of “survival ROC” in R software was utilized to generate ROC curve to assess the predictive ability of the established risk signature. P < 0.05 was considered be significant in all statistical tests.