An HDAC9-associated immune-related signature predicts bladder cancer prognosis

Background The close relationship between histone deacetylase 9 (HDAC9) and immunity has attracted attention. We constructed an immune signature for HDAC9, a vital epigenetic modification, to predict the survival status and treatment benefits in bladder cancer (BC). Methods An exhaustive analysis of HDAC9 and immunology via the tumor and immune system interaction database (TISIDB) was performed, and an immune prognostic risk signature was developed based on genes enriched in the top five immune-related pathways under high HDAC9 status. Comprehensive analysis of survival curves and Cox regression were used to estimate the effectiveness of the risk signature. The relationship between immunological characteristics and the risk score was evaluated, and the mechanisms were also explored. Results In the TISIDB, HDAC9 was closely related to various immunological characteristics. The risk signature was obtained based on genes related to prognosis enriched in the top five immune-related pathways under high HDAC9 status. The survival rate of the high-risk BC patients was poor. The risk score was closely related to multiple immunological characteristics, drug sensitivity, immunotherapy benefits and biofunctions. Conclusion An immune-related prognostic signature established for HDAC9 expression status could independently predict the prognosis of BC patients. The use of this signature could help clinicians make personalized treatment decisions.


Background
As the most common urinary tumor in the whole world, the morbidity of bladder cancer (BC) is increasing each year in China [1]. Once BC progresses to the advanced stage of metastasis, conventional treatments are no longer effective [2,3]. Furthermore, the clinical outcome of BC in the past few years has been poor because of high recurrence and drug resistance rates [4]. Thus, identifying new markers for prognosis and treatment for BC is necessary.
Histone deacetylation is considered a vital epigenetic modification in vivo. Histone deacetylases (HDACs) restrain related gene expression by catalyzing the deacetylation of histones or nonhistone proteins [5][6][7]. Histone deacetylase 9 (HDAC9) is a IIa HDAC subtype that has shown dual roles in tumors in multiple studies. In pancreatic ductal adenocarcinoma, high HDAC9 expression was significantly related to poor prognosis [8], and overexpression of HDAC9 promoted progression of the malignant phenotype of oral squamous cell carcinoma [9]. Additionally, elevated HDAC9 expression activated angiogenesis and invasion in triplenegative breast cancer [10]. However, HDAC9 also exhibited inhibitory effects in tumors. HDAC9 stimulated the expression of the ATDC target gene P53 by inhibiting ATDC expression, leading to inhibition of tumorigenesis [11], and it suppressed the proliferation of gastric cancer cells and enhanced sensitivity to cisplatin [12]. Our previous study suggested that low HDAC9 expression might facilitate clear cell renal cell carcinoma (ccRCC) cell growth, and we further found that HDAC9 was closely related to immunity and increased the infiltration levels of a variety of immune cells [13]. Other studies have revealed that HDAC9 deficiency significantly suppresses the immune response and inflammatory response [14][15][16]. Currently, immunotherapy has great potential in BC, and several immune-related drugs are currently in development [17]. Therefore, we sought to explore the pathogenesis of HDAC9 in BC, especially the relationship between HDAC9 and immunological characteristics.
In this study, a comprehensive analysis of HDAC9 was carried out to explore the connection between HDAC9 and the immune phenotype in BC. Importantly, we constructed an immune signature according to the expression status of HDAC9, a vital epigenetic modification, to predict the survival status and treatment benefits in BC.

Data acquisition from multiple public databases
The mRNA and clinically related factor data were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/), GSE32548 from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and the IMvigor210 trial [18]. TCGA data comprised the training cohort, GEO data were used as the validation cohort, and IMvigor210 trial data were used to evaluate the effect of immunotherapy. Immunohistochemical images of HDAC9 in BC samples and normal samples were retrieved from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/). The HDAC9 antibody product number used on the HPA website was HPA028926 (Atlas Inc., Stockholm, Sweden). Then, the connections between HDAC9 and the immune phenotypes were studied via the tumor and immune system interaction database (TISIDB) (http://cis.hku.hk/TISIDB/index.php) [19].
value (FDR Q value) < 0.25 were considered to indicate significant enrichment, and genes in the top five enriched immune-related biological pathways were extracted.

Risk score calculation
Genes related to survival were distinguished via univariable Cox regression. The P value was corrected for multiple comparisons by the FDR method, and a P adjusted value (P. adj) < 0.05 was considered statistically significant. Then, through least absolute shrinkage and selection operator (LASSO) Cox regression (the glmnet package of R), a risk score was calculated. The BC samples were split into two risk groups (low-risk vs. high-risk). A receiver operating characteristic (ROC) curve was utilized to analyze the risk signature accuracy for predicting prognosis via the R timeROC package, and the overall survival (OS) of the two groups was compared via a survival curve using the R survival package. To assess the ability of the risk score to independently predict prognosis, Cox regressions were then performed. The performance of prognosis prediction was also assessed through a nomogram and calibrations via the rms R package. To compare the predictive ability of our risk signature to other published signatures, we included risk signatures based only on prognostic-related genes [20][21][22][23] and those based on other risk factors [24][25][26][27][28] in the comparison of the AUCs of the ROC curves via the timeROC R package.

Coexpression analysis
Coexpression analysis was performed to further understand connections between HDAC9 and the risk score genes via the psych package of R. The relevant results were visualized using the ggplot2 R package. The P value was corrected for multiple comparisons by the FDR method, and P. adj < 0.05 was considered statistically significant.

The infiltration of immune cells and the risk score
To appraise the immune infiltration levels between two different risk groups, cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) was carried out. In a previous study, CIBERSORT could calculate the relative abundance of 22 immune cells using transcriptome data [29].

Tumor microenvironment (TME) and the risk score
The infiltration of stromal cells (also called stromal score), infiltration of immune cells (also called immune score) and tumor purity are considered important indicators of the TME, and these aspects were assessed via Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) [30]. Relationships between the ESTIMATE results and the risk score were then assessed. The scores of 13 specific stromal cells in BC patients from TCGA were also calculated by the xCell package of R software [31] and the relationships between the xCell results and the risk score were also evaluated. In addition, the connections between genes included in the risk signature and TME related parameters (including stromal cells, immune cells, stromal score, immune score and tumor purity) were also explored via the psych R package (the P value was corrected for multiple comparisons by the FDR method, and P. adj < 0.05 was considered statistically significant.). The above relevant results were visualized using the ggplot2 R package.

The drug sensitivity and risk score
The connections between the half-maximal inhibitory concentration (IC50) of six common chemotherapy drugs (cisplatin, docetaxel, methotrexate, gemcitabine, paclitaxel and doxorubicin) and the risk score were explored using the pRRophetic R package [32]. As responses, including stable disease (SD), progressive disease (PD), complete response (CR) and partial response (PR), of each BC case to anti-PD-1 therapy were included in the IMvigor210 trial data, the effects of immunotherapy in the two risk groups were also analyzed.

Functional assessment of the risk score
To explore the risk score-related biological functions, GSEA was performed. The "c2.cp.kegg. v7.4.symbols.gmt" gene set and "c5.all.v7.4.symbols.gmt" gene set of the MSigDB database were downloaded. NOM P value < 0.05 and FDR Q value < 0.25 were considered significant.

Statistical analysis
The Wilcoxon rank sum test was used to compare the significant differences between the two groups. Log-rank analysis (survival curve) was utilized for survival analysis, correlation analysis was performed with the Spearman test, and the chi-square test was used to assess differences in the responses of the two groups of patients to immunotherapy. The non-parametric test was performed to compare differences between AUCs of the ROC curves also via the time-ROC R package [33]. All statistical calculations were performed via R software (version 4.1.1).

Preliminary exploration of HDAC9 in BC
We included 411 tumor samples and 19 normal samples from TCGA, of which clinical information was available for 409 samples (Table 1). Compared with normal bladder tissues, HDAC9 was found to be lower in tumor tissues according to TCGA ( Fig 1A) and HPA ( Fig 1B  and 1C) data. It was worth noting that the significance level of the relationship between HDAC9 expression and immune subtypes in the BC ranked second among all tumors (Fig 1D  and 1E) [34]. HDAC9 was also closely related to various immunological characteristics [including tumor-infiltrating lymphocytes (TILs) (S1A-S1F Fig

Establishment of the risk signature
According to the top five immune-related pathways in the GSEA results, HDAC9 played both immunosuppressive and immune activation effects in BC (Table 2). In total, 593 genes in the five pathways were extracted, and genes related to survival (P. adj < 0.05) were identified via univariable Cox regression (S1 Table). Then, through LASSO Cox regression, a risk score was calculated. (Fig 2A and 2B) (S2 Table). The survival rate of the high-risk BC patients was lower than that of the low-risk patients (P < 0.001) (Fig 2C). The areas under the curve (AUCs) were 0.812 (1 year), 0.809 (3 years) and 0.813 (5 years), and the C index was 0.780 ( Fig 2D). The GSE32548 cohort, including 130 BC tissues and corresponding clinical data, was used to validate the signature (Fig 2E and 2F) ( Table 3). The survival rate of the high-risk BC patients was also lower (P = 0.009) (Fig 2G), with AUCs of 0.858 (1 year), 0.781 (3 years) and 0.786 (5 years) and a C index of 0.895 ( Fig 2H).
In the TCGA cohort, univariate Cox regression suggested that clinical stage, T stage, N stage and the risk score were significantly related to poor prognosis ( Fig 3A) (all P values < 0.05), and the risk score could be considered a prognostic factor based on multivariate Cox regression ( Fig 3B). In the GSE32548 cohort, univariate Cox regression revealed that T stage, grade and risk score were closely related to poor prognosis in BC (all P values < 0.05) (Fig 3C), and the risk score was also identified as a prognostic factor based on multivariate Cox regression ( Fig 3D). Additionally, the nomogram and calibration curves appropriately predicted 1-year, 3-year and 5-year OS for both TCGA (Fig 4A-4D) and GSE32548 (Fig 5A-5D). Upon comparison with other signatures constructed through only prognosis-related genes (Fig 6A-6C) or other risk factors (Fig 6D-6F), our risk signature was found to be superior (Tables 4 and 5).
The relationships between HDAC9 and selected genes were assessed by coexpression analysis (S5 Fig). CDK6 was identified as the most positively correlated with HDAC9 expression, whereas NR1H2 was the most negatively correlated with HDAC9 expression.

The risk score was related to infiltration of various immune cells
The results of CIBERSORT revealed enrichment for signatures of M0 ( Fig 7A) and M2 macrophages ( Fig 7B) and neutrophils ( Fig 7C) in high-risk BC patients (all P values < 0.05). Conversely, signatures of CD8 T cells (Fig 7D), activated memory CD4 T cells (Fig 7E), activated . The survival rate of the high-risk BC patients was worse than that of the low-risk patients (C). AUCs of the ROC curves were shown (D). In the verification cohort of the GEO, BC patients were distributed into different risk groups (E). Survival status of BC patients in different risk groups (F). The survival rate of the high-risk BC patients was worse than that of the low-risk patients (G). AUCs of the ROC curves were shown (H) TCGA, the Cancer Genome Atlas; GEO, Gene dendritic cells (Fig 7F), follicular helper T cells ( Fig 7G) and plasma cells ( Fig 7H) were enriched in low-risk BC patients (all P values < 0.05).

The risk score could affect the change of TME
The relationship between the TME (including stromal score, immune score and tumor purity) and the risk score was also evaluated. The results revealed that the stromal score was positively correlated with the risk score ( Fig 8A), but tumor purity was negatively correlated (Fig 8B) (all P values < 0.05). Nevertheless, there was no significant difference between the immune score and risk score ( Fig 8C). In addition, signatures of adipocytes (

Risk score and drug sensitivity
We predicted the IC50 of six common chemotherapy drugs in the different groups. Cisplatin ( Fig 10A) and docetaxel ( Fig 10B) exhibited higher IC50 values in low-risk patients; thus, the high-risk group was more sensitive to cisplatin and docetaxel (all P < 0.05). Methotrexate ( Fig  10C) and gemcitabine ( Fig 10D) exhibited a higher IC50 in high-risk BC patients; thus, lowrisk BC patients were more sensitive to methotrexate and gemcitabine (all P < 0.05). However, there were no significant differences between other chemotherapy drugs (paclitaxel and doxorubicin) and the risk score (Fig 10E and 10F). By analyzing data from the IMvigor210 trial (including 68 CR/PR patients and 230 SD/PD patients), we also found that the CR/PR rate in low-risk patients was higher than that in high-risk patients, which suggested that low-risk patients benefit significantly more from immunotherapy than high-risk patients ( Fig 10G).

Analysis of biological functions
The biological functions of the risk score were evaluated via GSEA. The most significant biofunctions enriched in the high-risk group based on Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses are listed in Tables 6 and 7, respectively. The most significant biofunctions enriched in the low-risk group based on KEGG analysis are listed in Table 8. Unfortunately, no significant pathways were found to be enriched in low-risk BC patients based on GO analyses.

Discussion
Currently, metabolic reprogramming, immune evasion and tumor-promoting inflammation are the three hallmarks of cancers [35]. As members of histone deacetylation metabolism, the role of HDACs in tumors has received widespread attention. As the HDAC IIa subtype, HDAC9 was found to play a dual role in tumors, which might be associated with the content of target proteins, alternative splicing, transcription factors that bind to the HDAC9 N-terminal region, and phosphorylation differences in multiple tumors [11,[36][37][38][39]. HDAC9 was initially found to be closely related to a variety of immunological parameters. And it was differentially expressed in five immune subtypes (previous research suggested that six immune subtypes could be identified in tumors, and there were significant differences in   immune infiltration and sensitivity to immunotherapy among subtypes, five immune subtypes could be identified in BC) [34], implying that HDAC9 could be used to classify immune subtypes, which would help individualized treatment. However, further follow-up work is required to improve these results.

AUCs 1 year (P value) 3 years (P value) 5 years (P value)
HDAC9Sig  found to regulate glucose metabolism [68][69][70][71][72]; MSH6 and TP53BP1 affected nucleotide metabolism by regulating DNA repair [73,74]; and NR1H2 and PHB could simultaneously influence both glucose and lipid metabolism [75,76]. Although immune-related pathways were activated under high HDAC9 expression, most enriched genes were metabolically related. In other words, changes in HDAC9, a regulator of protein posttranslational The risk score was significantly related to the TME based on the ESTIMATE algorithm. The stromal score was positively correlated with the risk score (A), and tumor purity was negatively correlated with the risk score (B). There was no significant difference between the risk score and immune score (C). TME, tumor microenvironment.
https://doi.org/10.1371/journal.pone.0264527.g008 modifications, could potentially be connected to metabolic disorders. Metabolic disorders might lead to changes in HDAC9 expression and activity, or metabolism and HDAC9 could both be regulated by the same upstream signal; this will require rigorous biological experiments to verify. Therefore, immunity and metabolism were closely linked in BC. By observing the coexpression of HDAC9 and the genes included in the risk signature and the degree of difference, more potential factors could be selected for subsequent research. By comparing with other risk signatures constructed only through prognosis-related genes, our risk signature was found to be superior in performance. One reason for this could be that the signature constructed by only prognosis-related genes was mixed with numerous factors that affect the prognosis of tumors (such as hypoxia and inflammation). The risk signature in our research used a metabolic factor (HDAC9) for stratification, and most of the immunerelated genes included in the risk signature were related to various metabolic processes, which effectively reduced confounding factors. We focused on the influence of metabolic processes on BC. Interestingly, by comparing BC signatures constructed with other risk factors, the signature we built still had advantages, which indicated that the metabolic process in BC should attract attention.
The connection between the risk score and TME was further analyzed. The ESTIMATE results showed that the risk score was positively correlated with the stromal score and negatively correlated with tumor purity, and there was no significant difference between the risk and docetaxel (B) exhibited higher IC50 values in low-risk patients; thus, the high-risk group was more sensitive to cisplatin and docetaxel. Methotrexate (C) and gemcitabine (D) exhibited a higher IC50 in high-risk BC patients; thus, low-risk BC patients were more sensitive to methotrexate and gemcitabine. However, there were no significant differences between other chemotherapy drugs (paclitaxel and doxorubicin) and the risk score (E-F). The CR/PR rate in low-risk patients was higher than that in high-risk patients (G), indicating more significant immunotherapy benefits for low-risk patients. IC50, half-maximal inhibitory concentration; SD, stable disease; PD, progressive disease; CR, complete response; PR, partial response; BC, bladder cancer; � , P <0.05; ���� , P < 0.0001.
https://doi.org/10.1371/journal.pone.0264527.g010 score and the immune score. However, the evaluation of specific immune cells through the CIBERSORT algorithm suggested that M0 and M2 macrophages and neutrophils were more enriched in the high-risk group, whereas CD8 T cells, activated memory CD4 T cells, follicular helper T cells, activated dendritic cells and plasma cells were more enriched in the low-risk group. M0 and M2 macrophages were positively correlated with a poor prognosis of BC [77,78], while a higher infiltration level of CD8 T cells was correlated with a good prognosis [79,80]. Therefore, the constructed risk score could appropriately assess immune cell infiltration, and it had an impact on immune cell infiltration. The lack of an association between overall immune score and risk score based on the ESTIMATE results suggested that the accumulation of negative effector immune cells in the high-risk group was close to that of positive effector immune cells in the low-risk group, resulting in no significant differences between the two groups. Additionally, the stromal score, which was significantly positively correlated with the risk score, could indicate that elevated stromal cells play an important role in different risk groups. Through further analysis, enrichment of signatures of adipocytes, chondrocytes, endothelial cells, fibroblasts, ly endothelial cells, MSCs, myocytes, and smooth muscle was found in the high-risk BC patients. Studies have shown that adipocytes promote breast cancer metastasis [81] and that chondrocyte can promote the progression of tumors by secreting inflammatory factors [82]. By combining the results of ESTIMATE analysis, we hypothesize that the decrease in tumor purity due to the increased stromal cell infiltration might contribute to the poor prognosis of high-risk patients.
Chemotherapy is an important treatment for BC, and cisplatin is considered to be the firstline chemotherapy for advanced BC [83]. However, due to drug toxicity and heterogeneity among patients, 40-60% of patients with bladder cancer do not respond to cisplatin-based systemic chemotherapy [84,85]. Studies have pointed out that the drug sensitivity of cisplatin in BC is closely related to ERKs and androgen receptors [86,87]. Despite the existence of various cancer treatment strategies, chemotherapy drugs are still widely used. Due to drug resistance and side effects, the therapeutic effects of chemotherapeutic drugs are often very different and difficult to predict, and the risk score calculated by various related risk factors could be significantly related to drug sensitivity [88]. The risk signature we built could develop personalized treatments by stratifying the risks of BC patients. Adjusting the drug dosage and choosing anticancer treatment strategies according to individual specificities could improve the prognosis of patients and minimize side effects. The six chemotherapy drugs included in the study affected the metabolic process to varying degrees. Among the chemotherapeutic drugs with differences between the two risk groups, cisplatin could induce cell cycle G1/S arrest [89], docetaxel promoted the formation of tubulin structures [90], methotrexate prohibited folic acid metabolism [91] and gemcitabine inhibited DNA replication [92]. The GO results indicated G1/S cell cycle arrest (GOBP_NEGATIVE_REGULATION_OF_CELL_CYCLE_G1_S_PHA-SE_TRANSITION), tubulin structure synthesis (GO_TUBULIN_BINDING), folic acid metabolism process (GO_FOLIC_ACID_METABOLIC_PROCESS), and the DNA replication (GO_POSITIVE_REGULATION_OF_DNA_REPLICATION) related biological functions were enriched in the high-risk group, implying that chemotherapy effectiveness depended on whether the drug mechanism conflicted with the metabolic changes of the patients themselves.
If they conflicted (such as methotrexate and gemcitabine), the effectiveness would be reduced, and vice versa (such as cisplatin and docetaxel). Chemotherapeutic drugs with similar mechanisms may mainly act on different metabolic changes in patients selectively, which might eventually lead to the opposite results. These should be the content of subsequent in vivo and in vitro precise pharmacological experiments. GSEA was conducted to explore biological mechanisms associated with the risk score. Based on GO analysis, a large number of stromal cells, junctions between cells and metabolic process-related biofunctions were significantly enriched in the high-risk BC patients. The results of KEGG revealed extracellular matrix (ECM), metabolic process and common BCrelated pathways (WNT, TGF-β and MAPK pathways) to be enriched in the high-risk phenotype [93][94][95][96][97][98]. The ECM is crucial for maintaining tissue homeostasis, and its disruption can promote tumor occurrence, progression, and metastasis by inducing EMT [99][100][101][102][103]. Hence, the GO, KEGG, ESTIMATE and xCell results were consistent.
Finally, our research has some limitations. First, this research was retrospective in design, and due to the comparison of patients from different cohorts, heterogeneity existed. Therefore, it is necessary to verify the prospective cohort used in the future. Second, the results of GSEA stratified by HDAC9 status and TISIDB suggested that HDAC9 could play a dual role in immune function, and further analysis should be conducted to clarify the net effect of HDAC9 on immune function and related mechanisms in future studies. Third, the inherent molecular mechanism of the risk signature affecting the sensitivity of tumor chemotherapy drugs is worthy of in-depth exploration.
In summary, an immune-related prognostic signature based on HDAC9 expression that could independently predict BC patient prognosis was constructed. The genes included in this risk score may constitute new targets for the treatment of BC, and further study is warranted. Moreover, the use of this signature will help clinicians make personalized treatment decisions.