MMP14 Contributes to HDAC Inhibition-Induced Radiosensitization of Glioblastoma

Glioblastoma (GBM) is the most common and malignant primary brain tumor in adults. Radiotherapy has long been an important treatment method of GBM. However, the intrinsic radioresistance of GBM cells is a key reason of poor therapeutic efficiency. Recently, many studies have shown that using the histone deacetylase (HDAC) inhibitor suberoylanilide hydroxamic acid (SAHA) in radiotherapy may improve the prognosis of GBM patients, but the underlying molecular mechanisms remain unclear. In this study, Gene Expression Omnibus (GEO) datasets GSE153982 and GSE131956 were analyzed to evaluate radiation-induced changes of gene expression in GBM without or with SAHA treatment, respectively. Additionally, the survival-associated genes of GBM patients were screened using the Chinese Glioma Genome Atlas (CGGA) database. Taking the intersection of these three datasets, 11 survival-associated genes were discovered to be activated by irradiation and regulated by SAHA. The expressions of these genes were further verified in human GBM cell lines U251, T98G, and U251 homologous radioresistant cells (U251R) by reverse transcription-quantitative polymerase chain reaction (RT-qPCR). It was found that MMP14 mRNA was considerably highly expressed in the radioresistant cell lines and was reduced by SAHA treatment. Transfection of MMP14 siRNA (siMMP14) suppressed cell survivals of these GBM cells after irradiation. Taken together, our results reveal for the first time that the MMP14 gene contributed to SAHA-induced radiosensitization of GBM.


Introduction
Glioblastoma (GBM) is the most common and aggressive primary brain tumor in adults [1]. Despite considerable advances in fundamental research on GBM that have been achieved in recent years, the median survival time for GBM patients following conventional therapy is only approximately 15 months [2]. For a long time, radiotherapy has been a major treatment strategy for GBM. However, the intrinsic radioresistance of GBM cells is the main reason for unsatisfactory therapy, and the expression of inherent radioresistance genes may be adversely correlated with the survival of GBM patients. Due to the poor prognosis of GBM patients and the limited effectiveness of traditional radiotherapy, there has been a steady increase in novel radioresistance-related biomarkers and radiosensitizers over the last decade.
Histone deacetylase inhibitor (HDACi), a tiny compound that can increase histone acetylation levels by inhibiting the activity of histone deacetylase (HDAC) [3], has the potential to increase tumor radiosensitivity [4,5]. HDACs are enzyme families whose enzymatic activity controls the acetylation state of protein lysine residues [6]. Acetylation of histones influences gene expression via influencing chromatin structure [4,7]. Inhibition

Evaluation of Radiation-Induced Changes in Gene Expression
To assess the changes in radiation-induced gene expression in GBM, we analyzed the GEO dataset GSE152982, which contains RNA-seq data generated from patient-derived GBM cells treated with 4 Gy radiation. RNA samples with three replicates from this dataset were extracted 48 h after irradiation. A cluster analysis of differentially expressed genes (DEGs) between irradiated and non-irradiated cells was performed ( Figure 1A), and a volcano plot was displayed ( Figure 1B), which showed that the gene expression patterns of these two groups were significantly different. In total, 428 genes were upregulated, and 588 genes were downregulated in irradiated cells compared with non-irradiated cells ( Figure 1C). A Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs was conducted, and the top ten KEGG enrichment pathways are shown in Figure 1D. Notably, many DEGs were involved in the cell cycle and p53 signaling pathways. In addition, we performed gene set enrichment analysis (GSEA) to functionally annotate the gene expression differences between irradiated and non-irradiated cells and revealed that the transcriptional signature of DEGs was enriched in the pathways of microtubule cytoskeleton organization, mitotic cell cycle, mitotic sister chromatid segregation, nuclear division, and sister chromatid segregation ( Figure 1E). Because cell cycle arrest and DNA

Screening of Survival-Associated Genes
GBM is classified as a grade IV glioma according to World Health Organization (WHO) definitions [22]. To identify survival-associated genes in GBM, we analyzed 249 high-grade glioma samples with WHO grade IV in dataset "mRNAseq_693" from the CGGA. Subjects with an overall survival <30% (overall survival <8.5 months) were classified as having short-term survival, whereas those with an overall survival >70% (overall survival >20 months) were classified as having long-term survival. DEG analysis of the two groups was performed using the Bioconductor package limma. Totally, 143 DEGs were identified, including 140 upregulated genes and 3 downregulated genes (Figure 2A). According to the KEGG pathway enrichment analysis, the survival-associated DEGs were strongly connected to the pathways of focal adhesion, phosphatidylinositol-3 kinase (PI3K)-Akt signaling pathway, extracellular matrix (ECM)-receptor interaction, and regulation of actin cytoskeleton pathways ( Figure 2B). A chord plot was created to provide an overview of the relationships between the genes linked to multiple processes ( Figure 2C). DEGs annotated to specific Gene Ontology (GO) terms could be identified using this representation. For example, MMP14, a significantly upregulated gene, was involved in the pathways of extracellular matrix organization, extracellular structure organization, and cell-substrate adhesion processes. A total of 143 DEGs were assigned to 50 sub-categories of GO terms from three main categories: cellular component (CC), biological process (BP), and molecular function (MF) ( Figure 2D). The top two sub-groups that contained most genes in the BP group were extracellular matrix organization (GO:0030198) and extracellular structure organization (GO:0043062). The top sub-group that contained most genes in the CC group was the collagen-containing extracellular matrix (GO:0062023). Moreover, the terms of extracellular matrix structural constituent (GO:0005201) and growth factor binding (GO:0019838) were associated with the greatest numbers of genes. The z-score was automatically computed using the Bioconductor package GOplot, which indicated whether a GO term was reduced (negative value) or increased (positive value). These analyses effectively limited the scope of the subsequent gene screening.

Analysis of Histone Acetylation-Controlled Genes
Histone acetylation plays an essential role in gene transcriptional regulation, and histone H3 lysine 27 acetylation (H3K27ac) is an epigenetic mark for active regulatory regions [23,24]. Therefore, we hypothesized that the expression of genes with H3K27ac gain at their promoters may be dramatically changed after SAHA treatment. To identify the potential genes controlled by histone acetylation, we downloaded the raw data of GEO dataset GSE131956, which contains H3K27ac chromatin immunoprecipitation-sequencing (ChIP-seq) of GBM [25]. The input control signal was subtracted from the ChIP signal before plotting. Annotating the characteristic peaks with ChIPSeeker [26], we found that 28.80% of the peaks were in the promoter regions ( Figure 3A). that the cell-substrate junction (GO:0030055) had the greatest numbers of genes in the CC group. Genes whose expression may be altered as a result of SAHA treatment were identified based on these analyses.    In addition, an enrichment of H3K27ac ChIP peaks around the transcription start sites (TSSs) was observed ( Figure 3B,D). After the annotation of peaks, 5474 genes with H3K27ac peaks at their promoters were extracted as candidate genes, and the top ten pathways of KEGG enrichment analysis are shown in Figure 3C. These candidate genes were closely related to Alzheimer's disease, amyotrophic lateral sclerosis, and Huntington's disease. Significantly enriched GO terms were also analyzed. Figure 3E illustrates that the transcription coregulator activity (GO:0003712) had the greatest number of genes in the MF group ribonucleoprotein complex biogenesis (GO:0022613), that the ncRNA metabolic process (GO:0034660) had the highest number of genes in the BP group, and that the cell-substrate junction (GO:0030055) had the greatest numbers of genes in the CC group.
Genes whose expression may be altered as a result of SAHA treatment were identified based on these analyses.

MMP14 Was the Potential Target of SAHA
To learn more about the major potential biomarker genes contributing to SAHAinduced radiosensitization of GBM, the intersection of the above three datasets was analyzed and 11 genes were obtained ( Figure 4A). "Radio_genes" indicates the DEGs after irradiation, "CGGA_genes" represents the survival-associated genes, and "ChIP_genes" represents the candidate genes that may change after SAHA treatment. H3K27ac was found to be highly enriched in the promoters of these 11 genes ( Figure 5). We hypothesized that these 11 genes may contribute to SAHA-induced radiosensitization in GBM.

MMP14 Was the Potential Target of SAHA
To learn more about the major potential biomarker genes contributing to SAHA-induced radiosensitization of GBM, the intersection of the above three datasets was analyzed and 11 genes were obtained ( Figure 4A). "Radio_genes" indicates the DEGs after irradiation, "CGGA_genes" represents the survival-associated genes, and "ChIP_genes" represents the candidate genes that may change after SAHA treatment. H3K27ac was found to be highly enriched in the promoters of these 11 genes ( Figure 5). We hypothesized that these 11 genes may contribute to SAHA-induced radiosensitization in GBM.  U251, U251R, and T98G cells after SAHA treatment. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001. Abbreviations: CGGA, Chinese Glioma Genome Atlas; ChIP, chromatin immunoprecipitation; H3K27ac, histone H3 lysine 27 acetylation; TSS, transcription start site; SAHA, suberoylanilide hydroxamic acid; IR, irradiated; Ctrl, control. To verify the bioinformatic analysis results, the in vitro experiments were performed using human GBM cell lines (U251 and T98G) and U251 homologous radioresistant cells (U251R). U251 and T98G cell lines with various radiosensitivities are widely used in the study of human GBM [27]. The U251R cell line was previously established from U251 us- To verify the bioinformatic analysis results, the in vitro experiments were performed using human GBM cell lines (U251 and T98G) and U251 homologous radioresistant cells (U251R). U251 and T98G cell lines with various radiosensitivities are widely used in the study of human GBM [27]. The U251R cell line was previously established from U251 using fractionated irradiation of X-rays in our laboratory, which showed an increased survival after irradiation in comparison with U251 [28]. These data suggested that the highly radioresistant cells might represent a poor clinical prognosis. It was found that T98G cells had the highest radioresistance among these cell lines, followed by U251R and U251 cells ( Figure 4B). To choose the best drug concentration of SAHA for further research, we evaluated the survival of these cell lines using the Cell Counting Kit-8 (CCK-8) assay at various concentrations and times ( Figure 4C) and chose the drug treatment conditions of 5 µM for 24 h. Western blot analysis revealed that after this treatment of SAHA, H3K27ac protein expression was lower in untreated cells than in SAHA-treated cells ( Figure 4D).
To further assess the potential involvement of the key candidate genes in the radiosensitivity of GBM cells, we examined the RNA expression levels of the key candidate genes in different GBM cell lines. Analysis of the CGGA datasets showed that RNA expression levels of 11 candidate genes were higher in the low-survival group (Table 1). Therefore, we hypothesized that the RNA expression levels of candidate genes would be higher in radioresistant cell lines. Consequently, the RNA expression levels of F3, IGFBP3, SOCS3, MMP14, and SPRY1 increased in U251, U251R, and T98G cells in an ordered manner ( Figure 6A), which was positively correlated with the radioresistance of the three cells ( Figure 4B). Furthermore, only GAL3ST4, SERPINE1, and PAMR1 expression was increased in T98G cells. However, there were no discernible changes in the RNA expression levels of PLAT, PBX3, PLAU, and SPRY1 in the three GBM cells. After 24 h of 5 µM SAHA treatment, there was a substantial increase in the RNA expression levels of F3, GAL3ST4, IGFBP3, PLAT, and SERPINE1, and there was a significant decrease in the RNA expression levels of MMP14, PBX3, and PLAU in the three cell lines ( Figure 6A). Moreover, there were no significant variations in the RNA expression levels of PAMR1 and SPRY1 in GBM cells after SAHA treatment ( Figure 6B). Furthermore, we verified the changes in the RNA expression levels of these candidate genes after irradiation ( Figure 6C). Because the RNA samples in GEO dataset GSE153982 were extracted 48 h after 4 Gy irradiation [29], we treated GBM cell lines with the same dose in this study. After irradiation for 48 h, the RNA expression levels of F3 and MMP14 were increased in the three cell lines, which was consistent with the results of our bioinformatics analysis of "Radio_genes" ( Table 1). The RNA expression levels of IGFBP3, SERPINE1, and PAMR1 were enhanced in the U251 and U251R cell lines. In U251R cells, the RNA expression level of GAL3ST4 was increased. PLAT levels decreased in U251R cells, whereas SOCS3 levels decreased in T98G cells. Furthermore, no significant changes in RNA expression of PBX3, PLAU, PAMR1, and SPRY1 were detected in the three cell lines.
According to the reverse transcription-quantitative polymerase chain reaction (RT-qPCR) data, MMP14 RNA expression was considerably higher in radioresistant cell lines and was reduced following SAHA treatment. In addition, the change in RNA expression level of MMP14 was compatible with the bioinformatics analysis of "Radio_genes". As a result, we hypothesized that the decrease in MMP14 expression was involved in the radiation sensitizing effect of SAHA in GBM. According to the reverse transcription-quantitative polymerase chain reaction (RT-qPCR) data, MMP14 RNA expression was considerably higher in radioresistant cell lines and was reduced following SAHA treatment. In addition, the change in RNA expression level of MMP14 was compatible with the bioinformatics analysis of "Radio_genes". As a result, we hypothesized that the decrease in MMP14 expression was involved in the radiation sensitizing effect of SAHA in GBM.

The Role of MMP14 in Radiosensitivity
MMP14 is highly expressed in numerous types of cancer, where it promotes angiogenesis, inflammation, cancer cell invasion, and metastasis [30]. Kaplan-Meier disease-free survival analysis showed that patients with high MMP14 expression in GBM had a significantly poorer prognosis than those with low MMP14 expression ( Figure 7A) (data from GEPIA database). Therefore, we assumed that MMP14 was linked to radiosensitivity in GBM.  SAHA had an obvious radiosensitizing effect on GBM cells ( Figure 7B-D). After SAHA treatment, the RNA expression level of MMP14 was significantly reduced ( Figure 6B). Furthermore, the protein expression level of MMP14 was reduced in three GBM cells after SAHA treatment, as determined by Western blot analysis ( Figure 7E). These results indicate that the application of SAHA and downregulation of MMP14 were largely related.
To further investigate the function of MMP14 in GBM, a small inferring RNA (siRNA) targeting MMP14 (siMMP14) and a scrambled small inferring negative control (siNC) were transfected into GBM cells. The transfection of siMMP14 significantly reduced the MMP14 levels in GBM cells compared with siNC ( Figure 7F). When MMP14 was knocked down by MMP14 siRNA, the survival of U251, U251R, and T98G cells was significantly decreased after irradiation ( Figure 7G). These results demonstrate that the decrease in the RNA expression level of MMP14 might be a reason for the radiosensitization of GBM cells with SAHA.

Discussion
GBM is a very aggressive tumor. Because of the limitations of standard therapy, a synergistic combination of therapies is required to counteract this aggressive brain tumor [31]. The drug development of radiosensitizer is an important link in the synergistic approaches to overcome current limitations of GBM treatment. HDACi, as represented by SAHA, is the most used radiosensitizing medicine, and it has also been used to treat GBM in conjunction with chemotherapy and radiation therapy [14]. HDACi can inhibit DNA repair responses and induce apoptosis, which contributes to enhanced sensitivity of tumor cells to chemotherapy and radiotherapy [32]. SAHA, for example, has been demonstrated to limit GBM development by reducing PRC2 function by decreasing EZH2 expression [33]. Furthermore, the application of SAHA suppressed c-Myc protein levels and extended animal survival in a patient-derived xenograft model [34]. However, the pharmacological targets of SAHA in the radiation sensitizing effect are poorly understood in GBM.
To identify the potential targets of SAHA for radiosensitization, we performed a series of bioinformatics analyses of GBM. Three databases were used to screen for potential targets. Containing RNA-seq data of the irradiated GBM cell line with three replicates, GEO dataset GSE153982 was analyzed to identify radiation-induced changes in gene expression. We believe that these "Radio_genes" are involved in the modulation of radiosensitivity in GBM. The CGGA database was used to screen for survival-associated genes to reduce the scope of the investigation. The intersection of the two groups was believed to be a possible radioresistant biomarker of GBM.
To further investigate the alteration of these candidate genes after SAHA treatment, H3K27ac ChIP-seq data (GEO dataset GSE131956) were analyzed. Depending on the gene target, histone acetylation marks can be present in the promoter regions where they can function by recruiting non-histone proteins to further influence gene expression [35]. Therefore, we hypothesized that genes with abundant H3K27ac marks in promoter regions would show notable alterations after SAHA treatment. Finally, the intersection of the three datasets yielded 11 potential genes. These 11 genes met three conditions, that is, their expression level changed after irradiation, their high expression was related to low survival of GBM, and their promoter regions had H3K27ac enrichment. Several previous studies have shown that IGFBP3, SOCS3, and MMP14 are related to radiosensitivity in certain types of cancers [36][37][38], which is consistent with our results that these three genes were highly expressed in radioresistant cells. Therefore, we assumed that these 11 genes were involved in SAHA-induced radiosensitization in GBM.
To verify the above conclusion, we examined the mRNA expression levels of 11 genes in three GBM cell lines. The results of RT-qPCR data confirmed that the RNA expression levels of these genes, except for PAMR1 and SPRY1, changed after SAHA treatment. Apart from PLAT, PBX3, and PLAU, the RNA expression levels of candidate genes were high in radioresistant GBM cell lines. Furthermore, the alterations in the RNA expression levels of F3, IGFBP3, MMP14, PAMR1, and SERPINE1 after irradiation were consistent with our bioinformatic study.
Since MMP14 is highly expressed in radioresistant cells and is decreased after SAHA treatment and since its expression trend after irradiation is compatible with the results of our bioinformatics study of GSE153982, we speculate that MMP14 is a possible target of SAHA. Based on GEPIA database analysis, a high MMP14 expression was associated with a poor prognosis in GBM patients. Some recent studies showed that the inhibition of MMP14 improves the efficiency of chemotherapy and radiotherapy by reducing the migration and invasion of cancer [20,36], which is consistent with the results in our study. Although MMP14 is an attractive target for cancer therapy, the development of drugs targeting MMP14 in vivo remains a difficult task.
Our study shows that the reduction in MMP14 can increase radiosensitivity in GBM, while the radiosensitizer SAHA can directly reduce the expression of MMP14. This indicates that the inhibition of MMP14 contributes to SAHA-induced radiosensitization in GBM. There exist several lines of evidence suggesting that SAHA and other HDACi have good application prospects as a radiosensitizer for GBM treatment [10,11]. The current study finds a novel function of SAHA in inhibiting the expression of MMP14 and provides basis for the clinical application of HDACi in the treatment of radioresistant GBM. However, the specific mechanism through which SAHA influences MMP14 expression requires further study. Although histone acetylation is typically associated with gene expression enhancement, there are instances where histone acetylation contributes to the suppression of gene expression. It has been shown that SAHA can block at or upstream of Akt in disrupting the PI3K/Akt pathway [39]. Additionally, the PI3K/Akt signaling pathway has been shown to activate MMP14 expression and enhance cell migration and invasion [40]. Therefore, we speculate that SAHA may repress the expression of MMP14 by blocking the PI3K/Akt pathway. Interestingly, histone acetylation, as it turns out, can also directly inhibit gene expression [41,42]. Some studies suggest that SAHA may restore the expression of a DNA damage repair gene called Ogg1 gene by modifying histone deacetylation and remodeling DNA methylation [43]. In other words, some epigenetic modifications caused by SAHA may inhibit MMP14 expression. Collectively, this study demonstrated that MMP14 was a key element in explaining the mechanism of SAHA treatment and provided an important direction for subsequent research.

Cell Culture and Irradiation
Human GBM cell lines U251 and T98G were purchased from the Cell Bank of the Chinese Academy of Sciences. The U251R cell line was previously developed from its parental cell line U251 by exposure to 2 Gy radiation/day for 30 fractions (5 fractions/week in general) with a total dose of 60 Gy [28]. U251, U251R, and T98G cells were cultured in DMEM (Gibco, Thermo Fisher Scientific, Waltham, MA, USA). All media were supplemented with 10% fetal bovine serum, 100 units/mL penicillin, and 100 mg/mL streptomycin (Gibco). Cells were incubated at 37 • C in 5% CO2 and subcultured every 3 days.
Cells in the log phase were irradiated with a dose rate of 0.883 Gy/min X-ray (X-RAD 320, PXI Inc., North Branford, CT, USA; 12 mA, 2 mm aluminum filtration) at room temperature.

SAHA Treatment
SAHA (SML0061, Sigma-Aldrich Inc., Shanghai, China) was freshly stocked in dimethyl sulfoxide (DMSO) at 5 mM. It was diluted into the medium of U251, U251R, and T98G cultures in the proper concentration, and the cultures were incubated without change of culture medium for proper hours. As the control of SAHA, U251, U251R, and T98G cells were incubated with 0.1% DMSO but without SAHA.

Colony Formation Assay
The radiosensitivity of U251, U251R, and T98G cells was assessed using a cell colony formation assay. Briefly, cells were plated in the six-well plates at a density of 150, 250, 450, 1000, and 2000 cells/well. After full attachment, they were exposed to 0, 2, 4, 6, and 8 Gy radiation. Approximately 2 weeks after radiation, cell colonies were fixed with methanol for 20 min and stained with 0.1% crystal violet for 30 min to count them. The number of colonies with more than 50 cells was determined. The colony formation assay was repeated three times for each cell sample. The cell survival curve was fitted using a single-hit multitarget model using GraphPad Prism 8 software (GraphPad Software, LLC, San Diego, USA). Sensitization enhancement ratio (SER) was expressed as an enhancement ratio determined at a survival fraction (SF) of 37%.

RNA Extraction and RT-qPCR Assay
The expression levels of common differentially expressed genes were measured by RT-qPCR in vitro. Total cellular RNA was extracted using TRIzol reagent (Invitrogen, San Diego, CA, USA). Total RNA (1 µg) was reverse transcribed into cDNA using a FastKing RT Kit (with gDNase) (Tiangen Biotechnology, Beijing, China). RT-qPCR was performed using SuperReal PreMix Plus (SYBR Green) (Tiangen Biotechnology, Beijing, China) in 20 µL reaction reagents using the MX3000P platform according to the manufacturer's protocol. The primers used in the RT-qPCR assays are listed in Supplementary Table S1.

Cell Proliferation Assay
Cell proliferation was measured using the CCK-8 assay (Dojindo Laboratories, Kumamoto, Japan) according to the manufacturer's instructions. Briefly, cells were seeded in 96-well plates at the appropriate concentration and cultured at 37 • C in an incubator for 4 h. When cells had adhered, 10 µL of CCK-8 working buffer was added to the 96-well plates and incubated at 37 • C for 2 h. The optical density at 450 nm was measured using a microplate reader (Tecan Infinite M200 Pro, Männedorf, Switzerland) to evaluate the cell proliferation rate. The CCK-8 assay for cell proliferation was repeated three times with three replicates for each cell sample.

Gene Expression Analysis of Clinical Database
The CGGA database (http://www.cgga.org.cn, accessed 30 March 2021) was used for this study. Data from mRNAseq_693 were used for survival-associated gene analyses. "Clinical Data" and "Expression Data" were downloaded. High-grade gliomas with WHO grade IV in the "mRNAseq_693" dataset were studied. Subjects with an overall survival <30% (overall survival <8.5 months) were classified as having short-term survival, whereas those with an overall survival >70% (overall survival >20 months) were classified as having long-term survival. Differential expression analysis of the two groups was performed using the Bioconductor package limma.

RNA Interference and Drug Treatment
U251, U251R, and T98G cells were transfected with MMP14 siRNA (siMMP14) using riboFECT CP Transfection Agent (RiboBio, Guangzhou, China) according to the manufacturer's protocol. The sequences of siRNAs are listed as follows: siMMP14 (5 GGT CTC AAA TGG CAA CAT A 3 ). The negative control siRNA had a random sequence. The transduction efficiency was consistently between 90% and 95%. U251, U251R, and T98G cells were treated with 5 mM SAHA for 24 h.

Statistical Analysis
All experiments were repeated at least three times, and the data are presented as mean ± standard deviation (SD). The differences between indicated groups were evaluated by Student's t-test or one-way analysis of variance (ANOVA) using GraphPad Prism 8 (GraphPad software, San Diego, CA, USA). RNA-seq analysis was performed using the Bioconductor package limma v3.36.0. Statistical significance was set at p < 0.05.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author upon reasonable request.