Abstract
High-dose radiation is the main component of glioblastoma therapy.
Unfortunately, radio-resistance is a common problem and a major
contributor to tumor relapse. Understanding the molecular mechanisms
driving response to radiation is critical for identifying regulatory
routes that could be targeted to improve treatment response. We
conducted an integrated analysis in the U251 and U343 glioblastoma cell
lines to map early alterations in the expression of genes at three
levels: transcription, splicing, and translation in response to
ionizing radiation. Changes at the transcriptional level were the most
prevalent response. Downregulated genes are strongly associated
with cell cycle and DNA replication and linked to a coordinated module
of expression. Alterations in this group are likely driven by decreased
expression of the transcription factor FOXM1 and members of the E2F
family. Genes involved in RNA regulatory mechanisms were affected at
the mRNA, splicing, and translation levels, highlighting their
importance in radiation-response. We identified a number of oncogenic
factors, with an increased expression upon radiation exposure,
including BCL6, RRM2B, IDO1, FTH1, APIP, and LRIG2 and lncRNAs NEAT1
and FTX. Several of these targets have been previously implicated
in radio-resistance. Therefore, antagonizing their effects
post-radiation could increase therapeutic efficacy. Our integrated
analysis provides a comprehensive view of early response to radiation
in glioblastoma. We identify new biological processes involved in
altered expression of various oncogenic factors and suggest new target
options to increase radiation sensitivity and prevent relapse.
Subject terms: Cancer, DNA damage and repair, Non-coding RNAs,
Transcriptomics, Translation, Cancer genomics, Functional genomics, RNA
splicing, Cancer, Genetics, Molecular biology, Oncology
Introduction
Glioblastoma is the most common intracranial malignant brain tumor with
an aggressive clinical course. Standard of care entails maximally safe
resection followed by radiotherapy with concomitant and adjuvant
temozolomide as described in the landmark European Organization for
Research and Treatment of Cancer (EORTC) Brain Tumor and Radiotherapy
Groups and the National Cancer Institute of Canada study^[41]1.
Adjuvant radiation therapy is commonly delivered 4–6 weeks after
maximally safe resection.. Nonetheless, the median overall survival
remains approximately 16 months^[42]1,[43]2, and the recent addition of
tumor-treating fields to the standard of care has only increased median
overall survival to 20.5 months^[44]1. Recurrence occurs in part
because glioblastoma uses sophisticated cellular mechanisms to repair
DNA damage from double-stranded breaks caused by ionizing radiation,
specifically homologous recombination and non-homologous end-joining.
Thus, the repair machinery confers a mechanism for resistance to
radiation therapy. Ionizing radiation can also cause base damage and
single-strand breaks, which are repaired by base excision and
single-strand break repair mechanisms, respectively^[45]3. A
comprehensive analysis of molecular mechanisms driving resistance to
chemotherapy and radiation is required to surpass major barriers and
advance treatments for glioblastoma.
The Cancer Genome Atlas (TCGA) was instrumental in improving the
classification and identification of tumor drivers^[46]4, but its
datasets provide limited opportunities to investigate radiation
response. Thus, studies using cell and murine models are still the best
alternatives to evaluate radiation response at the genomic level. The
list of biomarkers associated with radiation resistance in glioblastoma
is still relatively small. Among the most relevant are
FOXM1^[47]5,[48]6, STAT3^[49]6, L1CAM^[50]7, NOTCH1^[51]8, RAD51^[52]9,
EZH2^[53]10, CHK1/ATR^[54]11, COX-2^[55]12, and XIAP^[56]13. Dissecting
how gene expression is altered by ionizing radiation is critical to
identify possible genes and pathways that could increase
radio-sensitivity. A few genomic studies^[57]14–[58]16 have explored
this question, but these analyses were restricted to describing changes
in transcription.
Gene expression is regulated at multiple levels, and RNA-mediated
mechanisms such as splicing and translation are particularly relevant
in cancer biology. A growing number of inhibitors against regulators of
splicing and translation are being identified^[59]17. Splicing
alterations are a common feature across cancer types and affect all
hallmarks of cancer^[60]18. Numerous splicing regulators display
altered expression in glioblastoma (e.g. PTBP1, hnRNPH, and RBM14) and
function as oncogenic factors^[61]19. Importantly, a genome-wide study
using patient-derived models revealed that transformation-specific
depended on RNA splicing machinery. The SF3b-complex protein PHF5A was
required for glioblastoma cells to survive, but not neural stem cells
(NSCs). Moreover, genome-wide splicing alterations after PHF5A loss
appear only in glioblastoma cells^[62]20. Translation regulation also
plays a critical role in glioblastoma development. Many translation
regulators such as elF4E, eEF2, Musashi1, HuR, IGF2BP3, and CPEB1
promote oncogenic activation in glioblastoma, and pathways linked to
translation regulation (e.g., mTOR) promote cancer phenotypes^[63]21.
To elucidate expression responses to radiation, we conducted an
integrated study in U251 and U343 glioblastoma cell lines covering
transcription (mRNAs and lncRNAs), splicing, and translation. We
determined that the downregulation of FOXM1 and members of the E2F
family are likely the major drivers of observed alterations in cell
cycle and DNA replication genes upon radiation exposure. Genes involved
in RNA regulatory mechanisms were particularly affected at the
transcription, splicing, and translation levels. In addition, we
identified several oncogenic factors and genes associated with poor
survival in glioblastoma that displayed increased expression upon
radiation exposure. Importantly, many have been implicated in
radio-resistance, and therefore, their inhibition in combination with
radiation could increase therapy efficacy.
Results
We collected samples from U251 and U343 glioblastoma cells at 0 (T0),
1 hour (T1), and 24 hours (T24) post irradiation and identified all the
alterations at mRNA, splicing, and translation levels (Supplementary
Fig. [64]S1A). A Principal Component Analysis (PCA) performed on
RNA-Seq and Ribo-Seq read counts revealed that most variation can be
explained by the cell type along the first principal component, while
radiation time-related changes were captured along the second principal
component (Supplementary Fig. [65]S1B). The distribution of fragment
lengths for ribosome footprints was enriched in the 28–30 nucleotides
range (Supplementary Fig. [66]S2). while the Ribo-seq metagene profiles
exhibit high periodicity within the coding domain sequence (CDS) as
expected since ribosomes traverse three nucleotides at a time
(Supplementary Fig. [67]S3). All the replicates in each condition show
high level of similarity (Supplementary Fig. [68]S1B and [69]S4A).
Changes in global transcriptome profile in response to radiation
We first conducted an integrated analysis to evaluate the early impact
of radiation on the expression profile of U251 and U343 GBM lines.
Using this approach, we focus on genes undergoing up regulation or
downregulation post irraditaion after adjusting for cell line
differences (See Methods). A relatively small number of genes displayed
altered expression at T1 (Supplementary Fig. [70]S4B,C). Downregulated
genes are mainly involved in transcription regulation and include 18
zinc finger transcription factors displaying high expression
correlation in glioblastoma samples from TCGA (Supplementary
Table [71]S1). Upregulated sets contain genes implicated in cell cycle
arrest, apoptosis, and stress such as ZFP36, FBXW7, SMAD7, BTG2, and
PLK3 (Supplementary Table [72]S1).
Since many alterations were observed when comparing the T1 vs. T24 time
points (Supplementary Table [73]S1 and Supplementary Fig. [74]S4C), we
opted to focus on genes showing the most marked changes (log[2]
fold-change >1.0 or <−1.0 and adjusted p-value < 0.05) to identify
biological processes and pathways most affected at the T24 time point.
Top enriched GO terms and pathways among downregulated genes include
chromatin remodeling, cell cycle, DNA replication, and repair
(Fig. [75]1A). Additionally, we identified several GO terms associated
with mRNA metabolism, decay, translation, and ncRNA processing,
suggesting active participation of RNA-mediated processes in
radio-response (Fig. [76]1B). Network analysis indicated the set of
genes in these categories is highly interconnected (Fig. [77]1C and
Supplementary Table [78]S2).
Figure 1.
[79]Figure 1
[80]Open in a new tab
Characteristics of downregulated genes at 24 hours (T24) after
radiation exposure in glioblastoma cell lines. (A) Enriched gene
ontology related to cell cycle, DNA replication, and repair among
downregulated genes. (B) RNA-related Gene Ontology (GO) terms enriched
among downregulated genes summarized using REVIGO^[81]97. (C)
Protein-protein interaction network, according to STRING^[82]98 showing
downregulated genes associated with RNA-related functions. Gene
clusters based on the strength of connection and gene function are
identified by color. Lines colors indicate the type of association:
light green indicates an association based on literature findings; blue
indicates gene co-occurrence; magenta indicates experimental evidence.
To expand the expression analysis, we employed Weighted Gene
Co-expression Network Analysis (WGCNA)^[83]22 to identify gene modules
with significant co-expression variation in response to radiation. All
identified modules, along with the complete list of genes in each
module, are shown in Supplementary Fig. [84]S5A,B and Supplementary
Table [85]S3. Seven modules were identified
[MATH: (Zs
ummary>5)
:MATH]
as tightly regulated, independent of the cell line (Supplementary
Fig. [86]S5B). Among modules with the highest significant correlation
(0.8, p-value < 1e−7), module 2 contains genes downregulated in T24,
with many involved in cell cycle, metabolism mRNA metabolism,
processing, splicing, and transport (Supplementary Table [87]S3),
corroborating results described above.
Next, we investigated downregulated genes with the gene set enrichment
analysis (GSEA) tool Enrichr^[88]23 and conducted expression
correlation analysis with Gliovis^[89]24. Based on their genomic
binding profiles and effect of gene expression, FOXM1 and the E2F
family of transcription factors emerged as potential regulators of a
large group of cell cycle/DNA replication-related genes in the affected
set (Fig. [90]2A, Supplementary Table [91]S4). In agreement, E2F1,
E2F2, E2F8, and FOXM1 displayed a significant decrease upon radiation.
FOXM1 and E2F factors have been previously implicated in chromatin
remodeling, cell cycle regulation, DNA repair, and
radio-resistance^[92]25,[93]26. All four factors are highly expressed
in glioblastoma with respect to low-grade glioma. Importantly, they
display high expression correlation with a large set of downregulated
genes implicated in cell cycle and DNA replication and among themselves
in glioblastoma samples in TCGA (Fig. [94]2B,C). We corroborated the
changes in expression of these transcription factors and some of their
potential targets by qRT-PCR in U343 cells and also observed similar
changes in the glioblastoma line A172 and the glioma stem cell line
3565 (Supplementary Fig. [95]S6).
Figure 2.
[96]Figure 2
[97]Open in a new tab
E2Fs and FOXM1 in glioblastoma. (A) Correlation of E2F1, E2F2, E2F8,
and FOXM1 with target genes involved in cell cycle. (B) Expression
levels of E2F1, E2F2, E2F8, and FOXM1 in gliomas grades II, III, and IV
in TCGA samples. (C) E2F1, E2F2, E2F8, and FOXM1 expression correlation
in glioblastoma (TCGA samples) using Gliovis^[98]24. ***p-value <
0.0001.
Upregulated genes at T24 are preferentially associated with the
extracellular matrix receptor interaction pathway, extracellular matrix
organization, axonogenesis, and response to type I interferon
(Fig. [99]3A, and Supplementary Table [100]S2). With respect to the
extracellular matrix, we observed changes in the expression levels of
several collagens (types II, IV, V, and XI), glycoproteins of the
laminin family (subunits
[MATH: α,β,andγ :MATH]
), and also integrins (subunits
[MATH: α :MATH]
, and
[MATH: β :MATH]
) (Fig. [101]3B, and Supplementary Table [102]S1). Collagen type IV is
highly expressed in glioblastoma and implicated in tumor
progression^[103]27. In addition, it has been observed that the
activation of two integrins, ITGB3 and ITGB5, contributes to
radio-resistance^[104]28.
Figure 3.
[105]Figure 3
[106]Open in a new tab
Global view of upregulated genes at T24 post-radiation in glioblastoma
cells. (A) Gene ontology analysis of upregulated genes (B)
Protein-protein interaction networks according to STRING^[107]99
showing genes associated with extracellular matrix organization and
response to interferon. Gene clusters based on the strength of
connection and gene function are identified by color. Lines colors
indicate type of association: light green, association based on
literature findings; blue indicates gene co-occurrence; magenta
indicates experimental evidence.
Radiation treatment also induced the expression of genes involved in
neuronal differentiation and axonogenesis. Some key genes in these
categories include SRC, VEGFA, EPHA4, DLG4, MAPK3, BMP4, and several
semaphorins. These genes can have very different effects on
glioblastoma development, with some factors activating oncogenic
programs and others behaving as tumor suppressors. Similarly, type I
interferon’s effects on treatment are varied. For instance, interferon
inhibited proliferation of glioma stem cells and their sphere-forming
capacity and induced STAT3 activation^[108]29. On the other hand,
chronic activation of type I IFN signaling has been linked to adaptive
resistance to therapy in many tumor types^[109]30.
Activation of oncogenic signals post-radiation could counteract
treatment effects and later contribute to relapse. We searched the set
of highly up-regulated genes post-radiation for previously identified
radio-resistance genes in glioblastoma, oncogenic factors and genes
whose high expression is associated with poor prognosis (Supplementary
Table [110]S5). In Table [111]1, we list these genes according to their
molecular function. Since several of these genes have never been
characterized in the context of glioblastoma, our results open new
opportunities to prevent radio-resistance and increase treatment
efficiency. Importantly, there are inhibitors available against several
of these proteins (Supplementary Table [112]S4).
Table 1.
List of oncogenic factors, genes whose high expression is associated
with poor survival and genes previously associated with
radio-resistance in GBM that showed increased expression upon
radiation. Genes are listed according to molecular function.
Function Genes
Membrane protein AQP1, ARHGEF2, BAALC, CSF1R, CSPG4, EPS8L2, ERBB3,
FGFR4, FYN, GPM6A, ITGB3, JUP
Protein kinase ANKK1, CDKN1A, CSF1R, ERBB3, FAM20C, FGFR4, FYN, IKBKE,
MERTK, PDGFRB, SRC, TEC
Gene expression regulation ARID3A, ASAH1, BCL3, BCL6, CBX7, CEBPB,
ELF3, FAM20C, FEZF1, HOXA1, HOXB9, JUP, KDM5B, LMO1, LMO2, MACC1, MAF,
MSI1, MUC1, NKX2-1, PML, PRDM6, RORC, SATB1, SREBF1, TP53BP1, ZMYM2
Enzymatic activity ACSS2, AGAP2, APIP, ARHGEF2, C1R, CARD16, CD24,
CDKN1A, CEBPB, CSF1R, CSPG4, CTSZ, CUL7, CYTH4, EPS8L2, ERBB3, FAM20C,
FGFR4, FTH1, FUCA1, FYN, GHDC, IDO1, IKBKE, ITGB3, JUP, KDM5B, MCF2,
MERTK, MFNG, MRAS, NKX2-1,PDE6G,PDGFRB, PML, QPRT, RRM2B, SERPINA5,
SFN, SGSH, SRC, SREBF1, TEC, TGFB1, ZMYM2
Phosphotransferase CSF1R, ERBB3, FAM20C, FYN, IKBKE, MERTK, PDGFRB,
SREBF1, TEC
Cell surface receptor BMP7, CSF1R, ERBB3, FGFR4, FYN, ITGB3, ITGB5,
LRIG2, MCF2, MERTK, MFNG, PDGFRB, PRDM6, SRC, TEC, TRPM8
Metabolism regulation ACSS2, APIP, BCL3, BCL6, BTG2, CEBPB,CRTC1,
CSPG4, CTSZ, ELF3, FUCA1, IDO1, ITGB3, JUP, MAF, MFNG, NKX2-1, PARP3,
PRDM6, PTGES, QPRT, RRM2B, SGSH, TGFB1, TP53BP1, TRPM8, USP9X, VEGFA,
ZMYM2
[113]Open in a new tab
Changes in lncRNA profile in response to radiation
lncRNAs have been implicated in the progression of
glioblastoma^[114]31, but their role in response to ionizing radiation
is still poorly understood. We identified 161 lncRNAs with expression
alterations in T1 vs. T24 comparisons. Analysis of this set with
LnC2Cancer^[115]32 identified several lncRNAs aberrantly expressed in
cancer and with relevance to prognosis (Supplementary Table [116]S1).
We also detected significant downregulation of MIR155HG, whose high
expression is associated with glioma progression and poor
survival^[117]33. Another downregulated lncRNA with relevance to
prognosis is linc000152, whose increased expression has been observed
in multiple tumor types^[118]34. On the other hand, we observed a
significant upregulation of two “oncogenic” lncRNAs, NEAT1 and FTX.
NEAT1 is associated with tumor growth, grade, and recurrence rate in
gliomas^[119]35, while FTX promotes cell proliferation and invasion
through negatively regulating miR-342–3p^[120]36. Thus, if further
studies corroborate NEAT1 and FTX as players in radio-resistance,
targeting these lncRNAs should be considered to improve treatment
response.
Effect of radiation on splicing
Alternative splicing impacts genes implicated in all hallmarks of
cancer^[121]37 and is an important component of changes in expression
triggered by ionizing radiation^[122]38. All types of splicing events
(exon skipping, alternative donor, and acceptor splice sites, multiple
exclusive exons, and intron retention) were affected similarly upon
exposure to radiation (Supplementary Table [123]S6). At T24, we
observed that transcripts associated with RNA-related functions
(especially translation), showed the most splicing alterations.
Affected transcripts encode ribosomal proteins, translation initiation
factors, regulators of translation, and genes involved in tRNA
processing and endoplasmic reticulum. Other enriched GO terms include
mRNA and ncRNA processing, mRNA degradation, and modification.
Catabolism is another process associated with several enriched terms,
suggesting that splicing alterations in genes involved in catabolic
routes could ultimately contribute to apoptosis (Fig. [124]4A,B, and
Supplementary Table [125]S7). Changes in the splicing profile are
likely driven by an alteration in the expression of splicing
regulators. In Table [126]2, we show a list of splicing factors
displaying strong expression alterations. Among those previously
connected to glioblastoma development, LGALS3 is the most extensively
characterized. LGALS3 is a galactosidase-binding lectin and non-classic
RNA binding protein implicated in pre-mRNA splicing and regulation of
proliferation, adhesion, and apoptosis; LGALS3 also is a marker of the
early stage of glioma^[127]39.
Figure 4.
[128]Figure 4
[129]Open in a new tab
Impact of radiation on the splicing profile of glioblastoma cells. (A)
GO-enriched terms among genes showing changes in splicing profiles at
T24. GO-enriched terms are summarized using REVIGO^[130]97. (B)
Protein-protein interaction networks according to STRING^[131]98
showing genes associated with RNA-related functions whose splicing
profiles displayed alterations at T24. Gene clusters based on the
strength of connection and gene function are identified by color. Lines
color indicate type of association: light green, an association based
on literature findings; blue indicates gene co-occurrence; magenta
indicates experimental evidence.
Table 2.
Splicing regulators showing changes in expression 24 hours
post-radiation. Factors showing an increase in the expression are shown
in bold, while factors showing a decrease in the expression are
represented in italic.
Gene ID Gene name Function
AHNAK2 Protein AHNAK2 splicing regulation
ESRP1 Epithelial splicing regulatory protein 1 regulation of mRNA
splicing
LGALS3 Galectin-3 signaling receptor binding
NOVA2 RNA-binding protein Nova-2 alternative splicing regulation
SNRPN Small nuclear ribonucleoprotein N spliceosomal snRNP assembly
ALYREF THO complex subunit 4 RNA binding
DDX39A ATP-dependent RNA helicase DDX39A RNA helicase
GEMIN4 Gem-associated protein 4 rRNA processing
HNRNPL Heterogeneous nuclear ribonucleoprotein L alternative splicing
regulation
LSM2 U6 snRNA-associated Sm-like protein LSm2 U6 snRNA-associated
Sm-like protein
MAGOHB Mago nashi homolog 2 exon-exon junction complex
PPIH Peptidyl-prolyl cis-trans isomerase H ribonucleoprotein complex
binding
RBMX RNA-binding motif protein, X chromosome regulation of mRNA
splicing
SNRPD1 Small nuclear ribonucleoprotein Sm D1 spliceosomal snRNP
assembly
SNRPE Small nuclear ribonucleoprotein E spliceosomal snRNP assembly
SRSF2 Serine/arginine-rich splicing factor 2 regulation of mRNA
splicing
SRSF3 Serine/arginine-rich splicing factor 3 regulation of mRNA
splicing
TTF2 Transcription termination factor 2 transcription regulation
[132]Open in a new tab
Differential translational efficiency
We used Ribo-seq^[133]40 to identify changes in translation efficiency
triggered by radiation (Supplementary Table [134]S8). Translation,
protein localization, and metabolism appear as top enriched terms among
downregulated genes in T1 vs. T24 comparisons (Supplementary
Table [135]S9). In particular, several ribosomal proteins, along with
translation initiation factors and mTOR, showed a significant decrease
in translation efficiency (Fig. [136]5A,B). Overall, these results
indicate repression of the translation machinery post-radiation
exposure and its strong auto-regulation. Since changes in components of
the translation machinery are occurring at all levels (transcription,
splicing, and translation) at T24, we expect that major translational
alterations take place in later stages of post-radiation.
Figure 5.
[137]Figure 5
[138]Open in a new tab
Impact of radiation on the translation profile of glioblastoma cells.
(A) GO-enriched terms among genes showing changes in translation
efficiency at T24. GO-enriched terms are summarized using
REVIGO^[139]97. (B) Protein-protein interaction network, according to
STRING^[140]98 showing genes whose translation efficiency decreased at
T24. Gene clusters based on the strength of connection and gene
function are identified by color. Line colors indicate the type of
association: light green, an association based on literature findings;
blue indicates gene co-occurrence; magenta indicates experimental
evidence.
In the upregulated set, we highlight three genes FTH1, APIP, and LRIG2
that could potentially counteract the impact of radiation
(Supplementary Table [141]S10). FTH1 encodes the heavy subunit of
ferritin, an essential component of iron homeostasis^[142]41. Pang et
al.^[143]42, reported that H-ferritin plays an important role in
radio-resistance in glioblastoma by reducing oxidative stress and
activating DNA repair mechanisms. The depletion of ferritin causes
down-regulation of ATM, leading to increased DNA sensitivity towards
radiation. APIP is involved in the methionine salvage pathway and has a
key role in various cell death processes. It can inhibit
mitochondria-mediated apoptosis by directly binding to APAF-1^[144]43.
LRIG2 is a member of the leucine-rich and immunoglobulin-like domain
family^[145]44, and its expression levels are positively correlated
with the glioma grade and poor survival. LRIG2 promotes proliferation
and inhibits apoptosis of glioblastoma cells through activation of EGFR
and PI3K/Akt pathway^[146]45.
Crosstalk between regulatory processes
Parallel analyses of transcription, splicing, and translation
alterations in the early response to radiation provided an opportunity
to identify crosstalk between different regulatory processes. The
datasets showed little overlap, with just a few genes showing
alterations in two different regulatory processes. However, we
identified several shared GO terms when comparing the results of
alternative splicing, mRNA levels (transcription), and translation
efficiency (Supplementary Table [147]S10). These terms show two main
groups of biological processes. The first group indicates that the
expression of genes involved in DNA and RNA synthesis and metabolism is
particularly compromised. The second group is related to translation
initiation. Ribosomal proteins were particularly affected (Figs. [148]4
and [149]5). There is growing support for the concept of specialized
ribosomes. According to this model, variations in the composition of
the ribosome due to the presence or absence of certain ribosomal
proteins or alternative isoforms could ultimately dictate which mRNAs
get preferentially translated^[150]46. Therefore, these alterations
could later lead to translation changes of a specific set of genes.
Discussion
We performed the first integrated analysis to define global changes
associated with the early response to radiation in glioblastoma. Our
approach allowed the identification of “conserved” alterations at the
transcription, splicing and translation levels and defined possible
crosstalk between different regulatory processes. Alterations at the
level of transcription were dominant, but changes affecting genes
implicated in RNA mediated regulation were ubiquitous; they indicate
that these processes are important components in radio-response and
suggest that more robust changes in splicing and translation might take
place later.
E2F1, E2F2, E2F8, and FOXM1 as major drivers of transcriptional responses
upon radiation
We observed marked changes in the mRNA levels of genes implicated in
cell cycle, DNA replication, and repair 24 hours (T24) after radiation.
Downregulation of several transcription factors, most of them members
of the zinc finger family, was observed at one hour post-radiation.
This group displays high correlations in expression within glioblastoma
samples from TCGA, suggesting that they might work together to regulate
gene expression. Unfortunately, most are poorly characterized, and the
lack of information has prevented establishing further connections to
changes in the cell cycle and DNA replication that we observed at T24.
GSEA and expression correlation analysis suggested that the
downregulation of members of the E2F family is likely responsible for
several of the expression changes we observed at T24. E2Fs have been
defined as major transcriptional regulators of the cell cycle. The
family has eight members that could act as activators or repressors
depending on the context, and are known to regulate one another. They
are upregulated in many tumors due to overexpression of
cyclin-dependent kinases (CDKs), inactivation of CDK inhibitors, or RB
Transcriptional Corepressor 1 (RB1) and are linked to poor prognosis.
Alterations in E2F genes can induce cancer in mice^[151]47,[152]48.
Specifically, we found that three E2F members showed decreased
expression upon radiation: E2F1, E2F2, and E2F8, all of which have been
previously implicated in glioblastoma development.
E2F1 is probably the best-characterized member of the E2F family.
Besides its known effect on cell cycle regulation and DNA replication,
it is also a positive regulator of telomerase activity, binding the
TERT promoter^[153]49. Recent studies show that lncRNAs and miRNAs
function in an antagonistic fashion to regulate E2F1 expression,
ultimately affecting cell proliferation, glioblastoma growth, and
response to therapy^[154]50–[155]52. E2F2 has been linked to the
maintenance of glioma stem cell phenotypes and cell
transformation^[156]53,[157]54. Several tumor suppressor miRNAs (let7b,
miR-125b, miR-218, and miR-138) decrease the proliferation and growth
of glioblastoma cells by targeting E2F2^[158]53,[159]55–[160]57.
Although still poorly characterized in the context of glioblastoma,
E2F8 drives an oncogenic phenotype in glioblastoma. Its expression is
modulated by HOXD-AS1, which serves as a sponge and prevents the
binding of miR-130a to E2F8 transcripts^[161]58. FOXM1 is another
potential regulator of the group of cell cycle and DNA replication
genes affected by radiation. FOXM1 is established as an important
player in chemo- and radio-resistance and a contributor to glioma stem
cell phenotypes^[162]5,[163]6,[164]59–[165]64. FOXM1 and E2F protein
have a close relationship and share target genes^[166]65. Additionally,
FOXM1- and E2F2-mediated cell cycle transitions are implicated in the
malignant progression of IDH1 mutant glioma^[167]66.
E2F and FOXM1 targeting could be considered as an option to increase
radio-sensitivity. Since the development of transcription factor
inhibitors is very challenging, an alternative to be considered is the
use of BET (bromodomain and external) inhibitors. BET is a family of
proteins that function as readers for histone acetylation and modulates
the transcription of oncogenic programs^[168]67. Recent studies in
glioblastoma with a new BET inhibitor, dBET6, showed promising results
and established that its effect on cancer phenotypes comes via
disruption of the transcriptional program regulated by E2F1^[169]68.
RNA processing and regulation as novel categories in radio-response
Besides the expected changes in expression of cell cycle, DNA
replication and repair genes, radiation affected preferentially the
expression of genes implicated in RNA processing and regulation.
Additionally, we identified a co-expression module containing multiple
genes associated with translation initiation, rRNA and snoRNA
processing, RNA localization, and ribonucleoprotein complex biogenesis.
Many regulators of RNA processing are implicated in glioblastoma
development, and splicing alterations affect all hallmarks of
cancer^[170]69. Radiation-induced changes in the splicing patterns of
oncogenic factors and tumor suppressors such as CDH11, CHN1, CIC,
EIF4A2, FGFR1, HNRNPA2B1, MDM2, NCOA1, NUMA1, RPL22, SRSF3, TPM3, APC,
CBLB, FAS, PTCH1, and SETD2. We also observed changes in expression of
four RNA processing regulators previously identified in
genomic/functional screening for RNA binding proteins contributing to
glioblastoma phenotypes: MAGOH, PPIH, ALYREF, and SNRPE^[171]70.
Potential new targets to increase radio-sensitivity and prevent relapse
Activation of oncogenic signals is an undesirable effect of radiation
that could influence treatment response and contribute to relapse. We
observed increased expression or translation and splicing alterations
of a number of pro-oncogenic factors, genes whose high expression is
associated with poor survival and genes previously implicated in
radio-resistance.
Among genes with the most marked increase in expression upon radiation,
we identified members of the Notch pathway (HES2, NOTCH3, MFNG, and
JAG2). Notch activation has been linked to radio-resistance in
glioblastoma, and Notch targeting improves the results of radiation
treatment^[172]71,[173]72. We also identified several genes associated
with the PI3K-Akt, Ras, and Rap1 signaling pathways that increased
expression levels upon radiation exposure. Targeting these pathways has
been explored as a therapeutic option in glioblastoma^[174]71,[175]73.
Other oncogenic factors relevant to glioblastoma that had increased
expression after radiation exposure include SRC, MUC1, LMO2, PML, PDGFR
[MATH: β :MATH]
, BCL3, and BCL6.
Anti-apoptotic genes (BCL6, RRM2B, and IDO1) also showed increased
expression upon radiation. BCL6 is a member of the ZBTB family of
transcription factors, which functions as a p53 pathway repressor. The
blockage of the interaction between BCL6 and its cofactors has been
established as a novel therapeutic route to treat glioblastoma^[176]74.
RRM2B is an enzyme essential for DNA synthesis and participates in DNA
repair, cell cycle arrest, and mitochondrial homeostasis. The depletion
of RRM2B resulted in ADR-induced apoptosis, growth inhibition, and
enhanced sensitivity to chemo- and radiotherapy^[177]75. IDO1 is a
rate-limiting metabolic enzyme involved in tryptophan metabolism that
is highly expressed in numerous tumor types^[178]76. The combination of
radiation therapy and IDO1 inhibition enhanced therapeutic
response^[179]77.
Among genes whose high expression correlates with decreased survival in
glioblastoma, we identified several components of the “matrisome” and
associated factors (FAM20C, SEMA3F, ADAMTSL4, ADAMTS14, SERPINA5, and
CRELD1). The core of the “matrisome” contains ECM proteins, while
associated proteins include ECM-modifying enzymes and ECM-binding
growth factors. This complex of proteins assembles and modifies
extracellular matrices, contributing to cell survival, proliferation,
differentiation, morphology, and migration^[180]78. In addition,
several genes of the proteinase inhibitor SERPIN family (SERPINA3,
SERPINA12, SERPINA5, and SERPINI1) implicated in ECM regulation^[181]79
were among those with high levels of expression upon radiation.
In conclusion, our results generated a list of candidates for
combination therapy. Contracting the effect of oncogenic factors and
genes linked to poor survival could increase radio-sensitivity and
treatment efficiency. Importantly, there are known inhibitors against
several of these proteins (Supplementary Table [182]S5). Moreover, RNA
processing and translation were determined to be important components
of radio-response. These additional vulnerable points could be explored
in therapy, as many inhibitors against components of the RNA processing
and translation machinery have been identified^[183]80,[184]81.
Methods
Cell culture and radiation treatment
A172 was obtained from ATCC. U251 and U343 cells were obtained from the
University of Uppsala (Sweden) and maintained in Dulbecco’s Modified
Eagle Medium (DMEM, Hyclone) supplemented with 10% fetal bovine serum,
1% Penicillin/Streptomycin at 37 °C in 5% CO[2]-humidified incubators
and were sub-cultured twice a week. Cells were plated after appropriate
dilution, and ionizing radiation treatment was performed on the next
day at a dose of 5 Gray (Gy). 5 Gy was chosen since in current
hypofractionated regimen, radiation therapy is delivered over the
course of 5 fractions at 5 Gy per fraction. A cabinet X-ray system
(CP-160 Cabinet X-Radiator; Faxitron X-Ray Corp., Tucson, AZ) was used.
Glioma stem line 3565 was established in Dr. Jeremy Rich;s lab; cells
were cultured in Neurosphere Media consisting of Neurobasal-A
supplemented with B-27, 50 ng/ml EGF, and 50 ng/ml hFGF.
qRT-PCR analysis
After exposure to ionizing radiation, cells were cultured for 14, 24,
and 48 hours (hrs). Total RNA was extracted using the TRIzol reagent
(Invitrogen, Carlsbad, CA) according to the manufacturer’s
instructions. Reverse transcription was performed using High-Capacity
cDNA reverse transcription kit (Applied Biosystems). Quantitative PCR
was performed using TaqMan Universal PCR Master Mix (Applied
Biosystems) or PowerUp SYBR Green Master Mix (Applied Biosystems) and
reactions were performed on
[MATH:
ViiATM :MATH]
7 Real-Time PCR System (Applied Biosystems). Data were acquired using
ViiA 7 RUO software (Applied Biosystems) and analyzed using the 2-ΔΔ CT
method with GAPDH as an endogenous control. Probes and primers used in
qRT-PCR are listed in Supplementary Table [185]S11.
RNA preparation, RNA-seq and ribosome profiling (Ribo-seq)
RNA was purified using a GeneJet RNA kit from Thermo Scientific. The
TruSeq Ribo Profile (Mammalian) kit from Illumina was used to prepare
material for ribosome profiling (Ribo-seq). RNA-seq and Ribo-seq
samples were prepared according to Illumina protocols and sequenced at
UTHSCSA Genome Sequencing Facility.
Overall strategy to identify gene expression alterations upon radiation
To identify the most relevant expression alterations in the early
response to radiation, we analyzed samples from U251 and U343 cells
collected at 0 (T0), 1 (T1), and 24 (T24) hours post-radiation. The
time of one hour and 24 hours were chosen as this reflects the
immediacy of the DNA damage response and resolution of double stranded
DNA breaks. The production of double stranded DNA breaks occurs nearly
immediately with the subsequent resolution of the breaks nearly
immediately (within 15–60 minutes after exposure to ionizing radiation)
with near resolution of all breaks by 24 hours. To capture the
progressive dynamics of expression alterations, we compared T0 to T1
samples and T1 to T24 samples. Our strategy to identify the most
relevant alterations in expression with maximal statistical power was
to combine all samples and use a design matrix with cell type defined
as a covariate with time points (Supplementary Fig. [186]S1).
Sequence data pre-processing and mapping
The quality of raw sequences reads from RNA-Seq and Ribo-Seq datasets
were assessed using FastQC^[187]82. Adaptor sequences and low-quality
score (phred quality score <5) bases were trimmed from RNA-Seq and
Ribo-Seq datasets with TrimGalore (v0.4.3)^[188]83. The trimmed reads
were then aligned to the human reference genome sequence (Ensembl
GRCh38.p7) using STAR aligner (v.2.5.2b)^[189]84 with GENCODE^[190]85
v25 as a guided reference annotation, allowing a mismatch of at most
two positions. All the reads mapping to rRNA and tRNA sequences were
filtered out before downstream analysis. Most reads in the Ribo-seq
samples mapped to the CDS. The periodicity analysis was performed using
ribotricer^[191]86. The number of reads assigned to annotated genes
included in the reference genome was obtained by htseq-count^[192]87.
Differential gene expression analysis
For differential expression analysis, we performed counting over exons
for the RNA-seq samples. For translational efficiency analyses,
counting was restricted to the CDS. Differential gene expression
analysis was performed by employing the DESeq2 package^[193]88, with
read counts from both U251 and U343 cell samples as inputs. Both the
cell line and time were treated as covariates along with their
interaction term. To find differentially expressed genes that show
changes between two time points, we used time specific contrasts. The
p-values were adjusted for controlling for the false discovery rate
(adjusted p-value) using the Benjamini and Hochberg (BH)
procedure^[194]89. Differentially expressed genes were defined with an
adjusted p-value < 0.05.
Weighted gene co-expression network analysis
WGCNA^[195]22 uses pairwise correlations on expression values to
identify genes significantly co-expressed across samples. We used this
approach to identify gene modules with significant co-expression
variations as an effect of radiation. The entire set of expressed
genes, defined here as those with one or higher transcripts per million
higher (TPM), followed by variance stabilization) from U251 and U343
samples were clustered separately using the signed network strategy. We
used the
[MATH:
Zsummary
:MATH]
^[196]90 statistic as a measure of calculating the degree of module
preservation between U251 and U343 cells.
[MATH:
Zsummary
:MATH]
is a composite statistic defined as the average of the density and
connectivity based statistic. Thus, both density and connectivity are
considered for defining the preservation of a module. Modules with
[MATH:
Zsummary
>5 :MATH]
were considered as significantly preserved. The expression profile of
all genes in each co-expression module can be summarized as one
“eigengene”. We used the eigengene-based connectivity (kME) defined as
the correlation of a gene with the corresponding module eigengene to
assess the connectivity of genes in a module. The intramodular hub
genes were then defined as genes with the highest module membership
values (kME > 0.9). All analysis was performed using the R package
WGCNA. The protein-coding hub genes were then selected for gene
ontology enrichment analysis.
Translational efficiency analysis
We used Riborex^[197]91 to perform differential translational
efficiency analysis. The underlying engine selected was DESeq2^[198]88.
DESeq2 estimates a single dispersion parameter per gene. However,
RNA-Seq and Ribo-Seq libraries can have different dispersion parameters
owing to different protocols. We estimated the dispersion parameters
for RNA-Seq and Ribo-Seq samples separately and found them to be
significantly different (mean difference = 0.04, p-value
[MATH: <2.2e−16
:MATH]
). This leads to a skew in translational-efficiency p-value
distribution since the estimated null model variance for the Wald test
is underestimated. To address this issue, we performed a p-value
correction using fdrtool^[199]92 that re-estimates the variance using
an empirical bayes approach. Normalized read counts for both Ribo-seq
and RNA-seq samples are provided in Supplementary Table [200]S12.
Alternative splicing analysis
Alternative splicing analysis was performed using rMATS^[201]93. All
reads were trimmed using cutadapt^[202]94 with parameters (-u -13 -U
-13) to ensure trimmed reads had equal lengths (138 bp). rMATs was run
with default parameters in paired end mode (-t paired) and read length
set to 138 bp (-len 138) using GENCODE GTF (v25) and STAR index for
GRCh38.
Gene ontology (GO) and pathway enrichment analysis
To classify the functions of differentially enriched genes, we
performed GO enrichment, and the Reactome pathway^[203]95 analysis
using Panther^[204]96. For both analyses, we considered terms to be
significant if BH adjusted p-values weree <0.05, and fold enrichment is
>2.0. Further, we used REVIGO^[205]97 to reduce redundancy of the
enriched GO terms and visualize the semantic clustering of the
identified top-scoring terms. We used STRING database (v10)^[206]98 to
construct protein-protein interaction networks and determine
associations among genes in a given dataset. The interactions are based
on experimental evidence procured from high-throughput experiments
text-mining, and co-occurrence. Only high-confidence (0.70) nodes were
retained.
Expression correlation analysis
Gene expression correlation analysis was done using Gliovis^[207]24
using glioblastoma samples (RNAseq) from the TCGA. To select correlated
genes, we used Pearson correlation,
[MATH: R>0.3 :MATH]
, and p-value < 0.05. A list of genes affecting survival in
glioblastoma was downloaded from GEPIA^[208]99. A list of long
non-coding RNAs (lncRNAs) implicated in glioma development was obtained
from Lnc2Cancer^[209]32. Drug-gene interactions were identified using
the Drug-Gene Interaction Database^[210]100.
Accession codes
The processed data of read abundance matrices are available through GEO
accession [211]GSE141013. Code for differential expression analysis and
translational efficiency analysis is available at
[212]https://github.com/saketkc/2019_radiation_gbm.
Supplementary information
[213]Supplementary Table S1^ (1.5MB, pdf)
[214]Supplementary Table S2^ (449.6KB, xlsx)
[215]Supplementary Table S3^ (106.4KB, xlsx)
[216]Supplementary Table S4^ (551.7KB, xlsx)
[217]Supplementary Table S5^ (450.3KB, xlsx)
[218]Supplementary Table S6^ (110.1KB, xlsx)
[219]Supplementary Table S7^ (245.8KB, xlsx)
[220]Supplementary Table S8^ (49.2KB, xlsx)
[221]Supplementary Table S9^ (40.2KB, xlsx)
[222]Supplementary Table S10^ (80.2KB, xlsx)
[223]Supplementary Table S11^ (17.3KB, xlsx)
[224]Supplementary Table S12^ (51.9KB, xlsx)
[225]Supplementary information.^ (12.1MB, xlsx)
Acknowledgements