Histone Deacetylase-2 Serves as an Independent Prognostic Indicator in Osteosarcoma

Background Osteosarcoma is the most common primary malignant tumor of skeleton in adolescence. Histone deacetylase 2 (HDAC2), a member of class I histone deacetylase, is putatively involved in tumorigenesis of human malignancies. This study aimed to evaluate the expression pattern and prognostic value of HDAC2 in osteosarcoma. Methods Four datasets were obtained from the gene expression omnibus (GEO) database to explore the expression and prognostic value of HDAC2. Level 3 mRNA expression proles and clinical data were obtained in The Cancer Genome Atlas (TCGA) for validation. Expression pattern of HDAC2 were illustrated in GSE16088, GSE36001 and GSE42352. The prognostic value of HDAC2 was evaluated and validated by Kaplan-Meier analyses, receiver operating characteristic (ROC), concordance index (C-index) and calibration curve in GSE21257 and TCGA. Multivariate Cox regression analysis, nomogram, and decision curve analysis (DCA) were performed to assess the prognosis predictive capability. Protein-protein interaction (PPI) and gene set enrichment analysis (GSEA) were applied to further understand the molecular network and regulatory mechanisms. increased in osteosarcoma tissues. High HDAC2 expression was associated with tumor metastasis and chemotherapy ecacy. Kaplan-Meier analysis demonstrated that high HDAC2 predicted worse overall survival. The ROC curve showed good performance in survival prediction. Cox regression demonstrated that HDAC2 could be an independent prognostic indicator. GSEA revealed patients with high HDAC2 expression were enriched with multiple ontological signatures. Elevated expression of may identify an aggressive subgroup in osteosarcoma and serve as an independent prognostic indicator in these patients.


Backgrounds
Osteosarcoma, derived from primitive mesenchymal bone cells and characterized by high-metastasis potential, is the most common primary malignant tumor of skeleton in adolescence. 1 2 The percentage of osteosarcoma cases accounts for 2-3% of all childhood cancers, and the 5-year survival ranges from 65.7-69.6%, which is far more worse than other type of childhood cancers. 3 As osteosarcoma is a radioresistant tumor, surgical resection and systemic chemotherapy are the standard managements for osteosarcoma. However, due to the unknown molecular mechanisms of tumorigenesis, effective and precise therapeutic strategies still remain a challenge. 4 5 Therefore, the identi cation of novel bio-marks predicting the prognosis and treatment response is urgently needed.
Histone deacetylases are key enzymes that responsible for the deacetylation of lysine residues at the Nterminal regions of core histones (H2A, H2B, H3 and H4). Histone deacetylase 2 (HDAC2) belongs to class I histone deacetylase family, which is homology with yeast RPD3. Previous studies showed that HDAC2 was involved in transcriptional regulation, cell cycle progression and developmental events by Page 3/12 forming transcriptional repressor complexes with MDA, SIN3 and YY1. 6 7 HDAC2 also plays an important role in embryonic development, memory formation, synaptic plasticity and in ammatory diseases. [8][9][10][11] Elevated expression of HDAC2 was found in various cancerous diseases. Overexpression of HDAC2 in non-small cell lung cancer increased bronectin level and promoted cell invasion and metastasis. 12 In ovarian cancer cells, HDAC2 expression was increased after exposure to chemotherapy, and in higher order chromatin changes and cellular DNA damage responses. 13 However, the signi cance of HDAC2 in osteosarcoma remains unknown. In the present study, we investigated the expression levels, prognostic value, protein interaction network and enriched signaling pathways of HDAC2 via bioinformatics data mining. Our results showed that HDAC2 was overexpressed in osteosarcoma, and associated with tumor metastasis, chemotherapy e cacy. HDAC2 functioned well for estimating prognosis, and could serve as an independent prognostic indicator for the patients. A prognostic nomogram incorporating HDAC2 and chemotherapy response may accurately predict overall survival as well. HDAC2 interacted with multiple proteins, including MTA1, MTA2, CHD3, CHD4, RBBP4, RBBP7, SINA3A, EZH2, MBD2 and SMARCA5, and was involved with multiple tumor-associated pathways.
GSE21257 included a total of 53 osteosarcoma tissues with complete clinical information, and were used for correlation analysis and survival analysis, see Table 1 for details.

TCGA osteosarcoma cohort
The normalized gene expression data of TCGA osteosarcoma cohort, shown as log 2 (x+1) transformed RSEM normalized counts, were downloaded from the TCGA hub in USCS Xena (https://tcga.xenahubs.net). In total, 263 samples with survival information were enrolled in the study for survival analysis and gene set enrichment analysis.

Integration of protein-protein interaction (PPI) network construction
Protein-protein interaction analysis for HDAC2 were performed by STRING database (https://stringdb.org/) 18 . Based on the subsistent data, the protein interaction network was visualized by the Cytoscape software (http://cytoscape.org/) 19 . A node in the network denotes protein, and the edge denotes the interactions. Subsequently, the gene oncology (GO) enrichment and KEGG pathways data were obtained from STRING database and visualized in R software (version 3.6.1).
Gene set enrichment analysis (GSEA) GSEA was conducted to explore the high-HDAC2 associated signaling pathways on GSEA software (https://www.broadinstitute.org/gesa) 20 . The gene expression data in GSE21257 osteosarcoma cohort together with the gene sets in Hallmark MSigDB datasets were applied in GSEA.

Statistical analysis
Statistical analysis was performed by using the R software and GraphPad Prism (version 8.0). All group data are presented as mean ± SD. Two-group independent sample comparisons were preformed using Mann-Whitney U test. Pearson chi-square test was applied to determine the relationship between HDAC2 expression and the clinicopathological features. For survival analysis, Kaplan-Meier curve was generated with log-rank test. The cut-off value was optimized by using X-tile software (version 3.6.1) 21 . Timedependent ROC analysis was performed to assess the predictive accuracy of the predictors. Calibration curves were generated to explore the performance of survival prediction. The nomogram model was constructed based on the multivariate Cox regression. DCA curves was generated to assess the clinical practicability of the nomogram. P<0.05 was considered statistically signi cant.

Results
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 Fig 2a). Similarly, the higher expression of HDAC2 mRNA was remarkably correlated with better chemotherapy e cacy (P=0.048, Fig S2), but not with age, gender and metastasis (all P>0.05). Overall, elevated HDAC2 expression was signi cantly 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, timedependent 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.  (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 modi cation, chromatin organization, histone modi cation, 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).

Discussion
Novel therapeutic strategies such as target therapy and immunotherapy, have been developed for various tumor management, due to the deepen understanding of molecular and tumor microenvironment involved in oncogenesis. However, there is no breakthrough progress in the treatment of osteosarcoma in decades, and patients' survival remain poor. In the present study, we reported that HDAC2 is aberrantly expressed in osteosarcoma and associated with overall survival. Moreover, multiple proteins and mitosisassociated pathways are found to be correlated with HDAC2.
Histone deacetylases (HDACs) are known for promoting transcriptional repression and silencing. 22 The class I HDACs are the most frequently overexpressed in human malignancies, and alter cellular epigenetic programming to promote cell proliferation, metastasis and apoptosis. 23 24 The p300/YY1/miR-500a-5p/HDAC2 signaling axis regulates cell proliferation in colorectal cancer. 25 Dan Sun et al demonstrated that p53-HDAC2 axis plays a vital role in the regulation of the DNA damage response and contributes to chemo-sensitivity of osteosarcoma cells. 26 Jing Li et al reported that HDAC2 could upregulate the expression of IL-6 and trigger the migration of osteosarcoma cells, implying that targeted inhibition of HDAC2/IL-6 might be a potential approach for osteosarcoma therapy. 27 In this study, by combination analysis of multiple gene expression data of osteosarcoma from GEO, TCGA and Kaplan-Meier Plotter cohorts, we demonstrated that HDAC2 was overexpressed in osteosarcoma tissues, and acted as an unfavorable marker for patients' survival.
Through the protein-protein interaction network analysis, we found 10 genes involved with HDAC2.
Metastasis associated 1 (MTA1) and metastasis associated 1 family member 2 (MTA2) are proteins identi ed in metastatic cells and involved in cancer in ltration and metastasis. 28 Chromodomain Helicase DNA Binding Protein 3 (CHD3) and Chromodomain Helicase DNA Binding Protein 4 (CHD4) belongs to the CHD family which represents the main component of the nucleosome remodeling and deacetylase complex and plays an important role in epigenetic transcriptional repression. 29 RB Binding Protein 4 (RBBP4) and RB Binding Protein 7 (RBBP7) encode ubiquitously expressed nuclear proteins which belong to a highly conserved subfamily of WD-repeat proteins, and are the common subunit of histone methyltransferase PRC2 and histone deacetylase NuRD and SIN3A. 30 SIN3A is a subunit of histone deacetylase complex and related to tumor progression in breast cancer. 31 32 Translational dysregulation of enhancer of zeste homolog 2 (EZH2) is frequently observed in various cancer types.
Speci c strategies for the managing of EZH2 developed as precise cancer medicines. 33 Methyl-CpGbinding protein 2 (MBD2) is demonstrated to play a critical role in "rewriting" the cancer methylome by the maintenance and spread of DNA methylation at CpG islands and shores. 34 SMARCA5 partners with RSF-1 to compose the RSF complex, which belongs to the ISWI family of chromatin remodelers, is involved in cancer cell proliferation, invasion and migration. 35 36 These genes are associated with epigenetics regulation, and found to be important in tumor cell cycle control, DNA damage and repair. Thus, multiple ontology-associated hallmark pathways including regulation of nuclear division, nuclear receptor transcription, regulation of cAMP mediated signaling, and stem cell are enriched in osteosarcoma patients with high HDAC2 expression.

Conclusions
In conclusion, we demonstrate that HDAC2 expression levels are signi cantly elevated in osteosarcoma patients. High expression of HDAC2 correlates with tumor progression, and is an independent prognostic indicator in patients. Multiple epigenetic-associated genes and ontology-associated pathways exist crosstalk with HDAC2. Thus, elevated expression of HDAC2 may identify an aggressive subgroup in osteosarcoma and serve as an independent prognostic indicator in these patients.