Identification of radiation responsive genes and transcriptome profiling via complete RNA sequencing in a stable radioresistant U87 glioblastoma model

The absence of major progress in the treatment of glioblastoma (GBM) is partly attributable to our poor understanding of both GBM tumor biology and the acquirement of treatment resistance in recurrent GBMs. Recurrent GBMs are characterized by their resistance to radiation. In this study, we used an established stable U87 radioresistant GBM model and total RNA sequencing to shed light on global mRNA expression changes following irradiation. We identified many genes, the expressions of which were altered in our radioresistant GBM model, that have never before been reported to be associated with the development of radioresistant GBM and should be concertedly further investigated to understand their roles in radioresistance. These genes were enriched in various biological processes such as inflammatory response, cell migration, positive regulation of epithelial to mesenchymal transition, angiogenesis, apoptosis, positive regulation of T-cell migration, positive regulation of macrophage chemotaxis, T-cell antigen processing and presentation, and microglial cell activation involved in immune response genes. These findings furnish crucial information for elucidating the molecular mechanisms associated with radioresistance in GBM. Therapeutically, with the global alterations of multiple biological pathways observed in irradiated GBM cells, an effective GBM therapy may require a cocktail carrying multiple agents targeting multiple implicated pathways in order to have a chance at making a substantial impact on improving the overall GBM survival.


INTRODUCTION
Glioblastoma (GBM), with an overall survival of less than 1.5 years and a 5-year survival rate of 5%, is the most common, malignant primary cancer of the central nervous system in adults, with an estimated 12,120 new diagnoses in the United States alone in 2016 [1][2][3]. The current standard treatment regimen for GBM consists of maximal safe surgical resection, followed by the Stupp regimen, which consists of radiation therapy combined with concomitant and adjuvant temozolomide [4,5]. Despite this, recurrencecharacterized by radioresistance-is inevitable [6,7]. To elucidate the underlying mechanism of radioresistance, we recently established a stable, radioresistant GBM model and Research Paper also presented evidence partly attributing the radioresistance of GBM cells to radiation-induced upregulation of tumor promoters acid ceramidase (ASAH1) and sphingosine-1phosphate (S1P) [8,9]. ASAH1 is a lysosomal cysteine amidase that converts ceramides, which trigger senescence and apoptosis, into sphingosine and free fatty acids [10][11][12][13][14][15][16][17]. Subsequently, sphingosine is phosphorylated by sphingosine kinase 1 (SPHK1) or 2 (SPHK2) to form a tumor promoter, S1P [10][11][12][13][14][15]. However, the radiation effects on global gene expression at the messenger ribonucleic acid (mRNA) level in a stable radioresistant GBM model have never been investigated.
RNA sequencing (RNA-seq) has become a powerful technique for transcriptome profiling of cell lines of interest to explain phenotypic variations [18,19]. This study used RNA-seq to reveal many crucial radiationresponsive genes that may enable GBM cells to acquire resistance to radiation, providing the essential basis for further investigations of the role of these differential mRNA expressions in acquiring radioresistance.

Irradiation of GBM cells induced differential expression of 1094 radiation-responsive genes
We previously generated and described a stable radioresitant GBM model. Briefly, U87 cells received a total radiation of 10 Gy to generate radioresistant U87-10gy cells. Over the course of weeks, most cells died and less than ~1% of cells survived the irradiation. These radioresitant U87-10 gy cells were allowed to grow to confluence and were perpetuated for experiments. Total RNAs from U87 and U87-10 gy cells were harvested and subjected to further analysis. To screen for global mRNA changes following irradiation, we profiled transcriptomes of the control U87 cell line and its derivatives irradiated U87-10gy cells via RNA sequencing (Figure 1 and Supplementary Table 1). Using genes with twofold changes or greater with statistically significant values (P < 0.05) as criteria, we identified 1094 radiation responsive genes that were upregulated or downregulated in U87-10gy vs U87 cells. Among these 1094 radiation responsive genes, 427 were upregulated and 667 were downregulated (Supplementary Tables 2 and 3).

Upregulation of genes promoting tumor aggressiveness and invasion following irradiation
In an effort to categorize genes based on their functions, we performed a gene ontology analysis. This analysis revealed that upregulated genes were enriched in inflammatory response, cell migration, positive regulation of epithelial to mesenchymal transition, angiogenesis, cell proliferation, cell growth, positive regulation of the canonical Wnt signaling pathway, response to hypoxia, activation of mitogen-activated protein kinase activity, metalloendopeptidase activity, cellular response to fibroblast growth factor stimulus, cellular response to transforming growth factor beta stimulus, cellular response to tumor necrosis factor, and ribonuclease A activity genes (Tables 1 and 2, and Supplementary Tables 2 and 4).

Tumor suppressor, immune response, P-53dependent apoptosis, and cell adhesion genes are enriched in the pool of downregulated genes
The gene ontology analysis was conducted to analyze differentially expressed genes. Compared with control U87 cells, downregulated genes in irradiated U87-10gy cells were enriched in signal transduction, transcription DNAtemplated, metabolic pathway, apoptotic process, cell adhesion, intrinsic apoptotic signaling pathway by p53 class mediator, positive regulation of T-cell migration, positive regulation of macrophage chemotaxis, T-cell antigen processing and presentation, and microglial cell activation involved in immune response genes (Tables 1 and 3, and  Supplementary Tables 3 and 5).
Downregulation of tumor suppressor genes after irradiation has been described as a major mechanism for GBM cells to enhance their survival [41]. Highlighting the negative association with survival, RNA-seq revealed downregulation of the following representative well-known genes involved in the negative regulation of cell survival: S1PR1, PARP15, HOXA11, and ADGRG1 [42][43][44][45]. Improved prognosis in patients with GBMs has been associated with the high expression of S1PR1, suggesting its function as a tumor suppressor [42]. GBMs become less susceptible to radiation and chemotherapy when HOXA11 expression is suppressed [44]. The loss of ADGRG1's function promotes radioresistance of glioma-initiating cells [45].

DISCUSSION
Investigators have shown that irradiation of GBM cells can transform them into a more malignant form [6,7]. In their transcriptome analysis of glioma, Ma et al. describe that radioresistance of gliomas within hours following irradiation may be due to the inactivation of early proapoptotic molecules and late activation of to RNA sequencing. Experiments were performed in triplicates, and only genes with statically significant changes (~2000 genes) were utilized for the heat map. Green: relatively high expression, Red: relatively low expression. (P < 0.05) antiapoptotic genes [41]. It is important to note that in their study, the transcriptome analysis was performed within hours (6-48 hours) following radiation, while in our study, the analysis was done in a stable radioresistant GBM model that we recently developed [8]. Since their study performed analysis within hours following radiation, their findings likely include cells that would not survive radiation long-term. On the other hand, our stable radioresistant GBM model selected out these cells, and our final RNA-seq involved only truly radioresistant cells. Consistent with the Ma et al. study, the aberrant gene expressions observed in irradiated U87-10gy cells were enriched in genes involved in enhancing tumor malignancy and invasion. Similar to the Ma et al. study, we also found upregulation of antiapoptotic genes BNIP3 and SOD2 in irradiated U87-10gy cells. Exclusively enriched in upregulated genes are positive regulation of epithelial to mesenchymal transition, metalloendopeptidase activity genes, and response to hypoxia genes. Upregulated metalloproteases mRNA expressions include MME, MMP2, MMP3, MMP7, MMP12, ADAM9, and ADAM12.
Metalloproteases have been strongly implicated in encouraging tumor invasion and metastasis of many cancers by degrading the extracellular protein matrix [21,22]. Epithelial to mesenchymal transition, a process that overlaps with the acquirement of stem cell properties characterized by increased cell motility and resistance to chemo-and radiotherapy, is an important inducer of cancer stem-like phenotypes and is associated with an aggressive phenotype in glioma [46,47]. Epithelial to mesenchymal transition is typically induced by TGFB3, which was also upregulated in irradiated U87-10gy cells [47,48]. Hypoxia, which is frequent in GBM, is a major inducer of epithelial to mesenchymal transition as well and is also a main promoter of GBM invasion [47,49]. In GBM, hypoxia-inducible factor 1-alpha (HIF-1α) and carbonic anhydrase 9 expressions are induced by hypoxia and promote angiogenesis, migration, cell survival, proliferation, epithelial to mesenchymal transition, and radio-and chemoresistance [49,50]. In line with this, we found HIF-1α and carbonic anhydrase 9 expressions upregulated in irradiated U87-10gy cells ( Table 2).  Downregulated genes were enriched in tumor suppressor, positive regulation of immune response, P-53 dependent apoptosis, and cell adhesion genes (Tables 1  and 3). On the basis of gene ontology, we uncovered 16 genes associated with apoptosis, 10 genes associated with positive regulation of immune response, and 29 genes associated with cell adhesion that were downregulated (Tables 1 and 3). Ma et al. suggest that suppressing apoptotic potential through gene expression regulation of irradiated cells helps explain the radioresistant nature of irradiated GBM cells [41]. Consistent with our findings, multiple studies have identified downregulation of some of the apoptotic genes discovered in this study to play major roles in attenuating GBM apoptosis, especially for BBC3, DCC, BEX2, CASP1, IL1B, and SFRP2 [51-56]. For example, in their study of GBM cell migration and proliferation, investigators identified BBC3 as a potent inhibitor of cell migration and proliferation in GBM [51]. Downregulation of DCC expression was suggested as an important marker in tumor malignancy and recurrence in astrocytic tumors [52]. To further enhance their survival, GBM cells produce an immunosuppressive microenvironment to escape immune surveillance [57]. One mechanism to achieve this objective is through the secretion of transforming growth factor B (TGF-β) to block T-cell activation and proliferation [57]. In addition to TGF-β, we uncovered many other downregulated genes involved in the activation of the immune system, especially genes mediating T-cell antigen processing and presentation, that may help with immune evasion of radioresistant GBM cells (Tables 2  and 3).
A previous study of our recently described radioresistant GBM model revealed that radiation-induced accumulation of the ASAH1 protein level, as measured by immunoblotting, may enable GBM cells to survive radiation [8,9]. However, RNA-seq data demonstrated no changes in the mRNA expression of ASAH1 between U87 and U87-10gy cells. This discrepancy can be explained by the reported absence of the high degree of, or even negative correlation (~40% on average), between mRNA and protein expressions [19,[58][59][60][61][62]. Similar mRNA levels can produce varied levels of protein of interest depending on the regulations between the transcripts and protein product such as the ability of the cells to stabilize the mRNA, the rates of translation, the rates of protein degradation, etc. [19,59,60]. As seen in this study, irradiation of U87 cells resulted in significant gene expression changes, which may alter post-transcriptional regulation and ultimately affect the resultant protein expression level.
Therapeutically, with the global alterations of multiple biological pathways observed in irradiated GBM cells, an effective GBM therapy may require a cocktail carrying multiple agents targeting multiple implicated pathways in order to have a chance at making a substantial impact on improving the overall GBM survival. These findings of aberrantly expressed mRNAs following irradiation provide a crucial comprehensive starting point to understand the complex mechanism of radioresistance in GBM, and should be combined with immunoblotting or other techniques of direct measurement of protein levels to supply a more accurate picture of how cells can be altered to be radioresistant. Many mRNAs, whose expressions were altered in our radioresistant GBM model, have never before been reported to be associated with the development of radioresistant GBM and should be further investigated to understand their roles in radioresistance. The absence of a major progress in the treatment of GBM is partly attributable to our poor understanding of both GBM tumor biology and the acquirement of treatment resistance in recurrent GBMs. Recurrent GBMs are characterized by their resistance to radiation. In this study, we used an established stable radioresistant GBM model to shed light on global mRNA expression changes after irradiation. Many mRNAs, the expressions of which were altered in our radioresistant U87 GBM model, have never before been reported to be associated with the development of radioresistant GBMs and should be concertedly further investigated to understand their roles in radioresistance.

Reagents and cells
The U87 and U87-10gy glioblastoma cell lines were cultured in Eagle's minimum essential medium containing 10% (v/v) fetal bovine serum. Culture medium materials were obtained from Life Technologies, Inc. (Grand Island, NY, USA).

RNA library preparation and sequencing
RNA sequencing libraries were prepared using the TruSeq Stranded mRNA Library Prep Kit (Illumina, Inc., San Diego, CA, USA) according to the manufacturer`s protocol. The RNA concentration was measured with a NanoDrop 2000c spectrophotometer (Thermo Scientific Inc., Waltham, MA, USA). Integrity was assessed using an Agilent 2200 TapeStation instrument (Agilent Technologies, Santa Clara, CA, USA). Briefly, polyA mRNA from an input of 500 ng high-quality total RNA (RINe>8) was purified and fragmented. First strand complementary deoxyribonucleic acid (cDNA) syntheses were performed at 25° C for 10 minutes, 42° C for 15 minutes, and 70° C for 15 minutes, using random hexameres and ProtoScript II Reverse Transcriptase (New England BioLabs Inc., Ipswich, MA, USA). In a second strand cDNA synthesis, the RNA templates were removed and a second replacement strand was generated by incorporating deoxyuridine triphosphate (in place of deoxythymidine triphosphate, to keep strand information) to generate ds cDNA. The blunt-ended cDNA was cleaned up from the second strand reaction mix with beads. Next, the 3′ ends of the cDNA were adenylated and then indexing adaptors were ligated. Polymerase chain reactions (15 cycles of 98° C for 10 seconds, 60° C for 30 seconds, and 72° C for 30 seconds) were used to selectively enrich those DNA fragments that have adapter molecules on both ends and to amplify the amount of DNA in the library.
The libraries were quantified using the Promega QuantiFluor dsDNA System on a Quantus Fluorometer (Promega, Madison, WI). The size and purity of the libraries were analyzed using the High Sensitivity D1000 Screen Tape on an Agilent 2200 TapeStation instrument. The libraries were normalized, pooled, and subjected to cluster, and pair read sequencing was performed for 150 cycles on a HiSeq4000 instrument (Illumina, Inc., San Diego, CA, USA), according to the manufacturer's instructions.

Generation of the stable radioresistant GBM model
The method and detail to generate the stable radioresistant GBM model was previously described by us [8]. U87 cell lines were grown to confluence and then irradiated with a Pantak HF320 X-ray machine (Agfa NDT Ltd., Reading, UK) operating at 300 kV at a dosage of 2.09 Gy/min to a total radiation dose of 10 Gy, to generate the U87-10gy cell lines. Following radiation, these irradiated cells were allowed to recover and to grow to confluence. The U87-10gy cell line was stable for continued passages for further studies.