Expression pattern of HDAC family in osteosarcoma and normal tissues
HDAC family members were grouped into four distinct classes based on homology and structure. Class I HDACs (HDAC1/2/3/8) were highly homologous to Rpd3, class II HDACs (HDAC4/5/6/7/9/10) were homologous to Hda1, and class IV HDAC (HDAC11) was similar to class I/II zinc-dependent HDACs. Class III HDACs (Sirt 1-7) were silent information regulator-2 related enzymes, which require NAD+ instead of zinc as a cofactor. The expression of HDAC family members were evaluated in osteosarcoma and normal tissues (Fig S1). Comparing with the normal tissues, the osteosarcoma tissues showed higher expression level of HDAC2 in GSE16088 (P<0.001, Fig 1a), GSE36001 (P<0.001, Fig 1b) and GSE42352 (P=0.004, Fig 1c). These results indicate that elevated expression of HDAC2 might play an oncogenetic role in osteosarcoma.
Correlation between HDAC2 expression and clinicopathological characteristics
The relationship between HDAC2 expression and clinicopathological characteristics was evaluated in GSE21257 and TCGA cohort. High HDAC2 mRNA expression was positively correlated with tumor metastasis (P=0.007, Fig 2a and b), but not with age, gender, Huvos grade, tumor location and histological subtypes (all P>0.05, Fig 2a). Similarly, the higher expression of HDAC2 mRNA was remarkably correlated with better chemotherapy efficacy (P=0.048, Fig S2), but not with age, gender and metastasis (all P>0.05). Overall, elevated HDAC2 expression was significantly associated with advanced clinicopathological characteristics.
Prognostic value of HDAC2 in osteosarcoma patients
Kaplan-Meier analysis was applied to evaluate the prognostic value of HDAC2 mRNA expression in osteosarcoma patients. Patients with high HDAC2 expression had remarkably worse overall survival than those with low HDAC2 expression in GSE21257 (HR=2.82, 95%CI (1.16-6.90), P=0.0227, Fig 3a) and TCGA (HR=1.87, 95%CI (1.24-2.81), P=0.003, Fig S3a) cohorts. Patients with high HDAC2 expression possessed high mortality and shorter survival time (Fig 3b and Fig S3b). In GSE21257 cohort, time-dependent ROC curve for 3- and 5-year overall survival predictions were generated, and the AUC were 0.715 and 0.662, respectively (Fig 3c). The C-index of HDAC2 was 0.626. The calibration curve for validation of the 3-year survival prediction ability of HDAC2 revealed that the predicted overall survival accorded with the observed overall survival (Fig 3d). Similar results were obtained in TCGA cohort (Fig S3b-e). Generally, HDAC2 performed well at predicting overall survival of osteosarcoma.
HDAC2 served as independent prognostic indicator in osteosarcoma
Prognostic factors of overall survival for osteosarcoma were identified by using univariate and multivariate Cox regression analyses. In GSE21257, age (>17 years vs ≤17 years), gender (male vs female), Huvos grade (grade 4, 3, 2 vs grade 1), location1 (upper limbs vs lower limbs), location2 (distal limbs vs proximal limbs), histology (Fibroblastic, chondroblastic and others vs osteoblastic), and HDAC2 expression (high vs low). The univariate Cox regression model revealed that Huvos grade ((HR=0.496, 95%CI (0.284-0.867), P=0.014) and HDAC2 ((HR=2.824, 95%CI (1.156-6.898), P=0.023) were significantly associated with overall survival of osteosarcoma (Table 2). Multivariate Cox regression model revealed that Huvos grade and HDAC2 were independent prognostic indicator of overall survival (both P<0.05, Table 2). Consistent results were obtained in TCGA cohort (Table S1).
Construction and validation of a predictive nomogram
A prognostic nomogram was constructed to predict 3- and 5-year overall survival based on the multivariate Cox regression model (Fig 4a and Fig S4a). The C-index of the nomogram were 0.748 and 0.891 in GSE21257 and TCGA cohorts, respectively, which was superior to HDAC2 in terms of predicting overall survival in osteosarcoma. The calibration curve for validation of the 3-year survival prediction ability of nomogram revealed that the predicted overall survival accorded with the observed overall survival (Fig 4b and Fig S4b). Decision curve analysis showed better prediction ability of nomogram than Huvos grade or HDAC2 alone (Fig 4c and Fig S4c).
Protein-protein interaction network
Protein-protein interaction network analysis using STRING database revealed that 10 genes including MTA1, MTA2, CHD3, CHD4, RBBP4, RBBP7, SIN3A, EZH2, MBD2 and SMARCA5 were interacted with HDAC2 (Fig 5a). The network was consisted of 11 nodes and 54 edges. The nodes were represented the proteins and the lines were represented the interaction between them.
Functional annotations and signaling pathways enrichment
Functional enrichment analyses of the 11 involved genes were performed and visualized in bubble chart in GSE21257 cohort. The top ten enriched biological process (BP) of the 11 genes were involved in histone deacetylation, covalent chromatin modification, chromatin organization, histone modification, regulation of gene expression, epigenetic, negative regulation of gene expression, epigenetic, regulation of transcription, DNA-Templated, transcription, DNA-templated, Regulation of transcription by RNA polymerase II, and regulation of RNA metabolic process (Fig 5b). The top ten enriched molecular function (MF) of the 11 genes were involved in protein deacetylase activity, histone deacetylase activity, RNA polymerase II repressing transcription factor binding, chromatin binding, nucleic acid binding, RNA polymerase II transcription factor binding, DNA binding, catalytic activity, acting on a protein, hydrolase activity, and histone deacetylase binding (Fig 5c). The top ten enriched cellular component (CC) of the 11 genes were involved in histone deacetylase complex, nuclear transcriptional repressor complex, transcriptional repressor complex, nuclear chromatin, NuRD complex, SWI/SNF superfamily-type complex, nucleoplasm part, catalytic complex, ESC/E(Z) complex, and sin3 complex (Fig 5d). The Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathways of the 11 genes were involved in thyroid hormone signaling pathway, Huntingtons disease, transcriptional misregulation in cancer, viral carcinogenesis, and human papillomavirus infection (Fig 5e). Gene set enrichment analysis was conducted between high HDAC2 group and low HDAC2 group to identify the biological signaling pathways. Multiple ontology-associated hallmark pathways including regulation of nuclear division, nuclear receptor transcription, regulation of cAMP mediated signaling, and stem cell (Fig 6).