In silico Analysis of Transcription Factors Associated to Differentially Expressed Genes in Irradiated Glioblastoma Cell Lines

1. Introduction

Glioblastoma multiforme (GBM) is one of the most frequent tumors in the central nervous system and the most malignant tumor among gliomas. In the past two decades, cytogenetic and molecular genetic studies have identified a number of recurrent chromosomal abnormalities and genetic alterations in malignant gliomas, particularly in GBM [1]. It was already described that GBM harbors combinations of the following genetic alterations: loss of heterozygozity of 10q, EGFR amplification, TP53 mutations, p16^INK4a deletion and PTEN mutations [2]. New integrative genomics studies provided a comprehensive view of the complicated genomic landscape of GBM, revealing a set of core signaling pathways commonly activated in GBM involving TP53, RB, and RTK (receptor tyrosine kinase) pathways [3, 4]. The majority of GBM tumors present genetic alterations in all three pathways, which helps to stimulate cell proliferation and enhance cell survival while allowing tumor cells to escaping from cell-cycle checkpoints, senescence, and apoptosis. This approach also identified previously unknown genetic alterations in IDH1/2, NF1, ERBB2, and NFKBIA genes [1].

The current GBM treatment involves aggressive management including surgery, adjuvant temozolomide-based chemotherapy, and radiotherapy [5], but GBM patients still present a dismal prognosis, and the median survival is 14.6 months from diagnosis [6]. Although radiotherapy has been found to significantly prolong survival rates for GBM patients, radioresistance is a typical characteristic of this tumor [7].

Current genome-wide studies and the molecular characterization of GBM have allowed the identification of potential new targets, development of novel therapeutic small molecules and monoclonal antibodies and initiation of clinical trials with these targets [6, 8-10]. However, there is a wide molecular diversity and heterogeneity associated with the aberrantly GBM signaling pathways, culminating in the relative lack of success of these new approaches [10]. Recently, an alternative strategy involves the selective targeting of GBM stem cells, which are resistant to chemo- and radiotherapy. But still, almost all small-molecule inhibitors designed to target these cells failed to demonstrate the effectiveness of this strategy, compared with the conventional therapy [11].

Considering that most of the treatment protocols are still ineffective, novel approaches are needed towards killing of GBM cells. Transcription machinery, as well as its regulatory elements is also a feasible new target for the application of molecular therapies. Transcription of DNA is dependent on the spatially and temporally coordinated interaction between transcriptional machinery involving RNA polymerase II, transcription factors (TFs)) and transcriptional regulatory components (promoter elements, enhancers, silencers and locus control regions) [12, 13]. The low level of transcription, directed by the general transcription factors associated to RNA polymerase core enzyme, is known as basal transcription [14]. However, there is a rapidly expanding number of ‘context-dependent’ transcription factors that bind DNA and these TFs are capable of positively or negatively regulating the transcription process depending on the context of their binding sites, the complement of protein interactions and other environmental influences [15].

Postgenomic analyses of major transcription factor families, in both malignant and nonmalignant cell types, have opened new discussions about TF function. The mechanisms by which TFs act in cancer cell systems appear to exhibit a restricted repertoire of skills and plasticity displayed by normal cell systems [16]. The evolution of a restricted malignant transcriptome can be seen clearly in the nuclear receptor superfamily, but is also apparent in the MYC and AP-1 networks [17]. Oncogenic transcriptional rigidity reflects the simultaneous deregulation of target loci such that proliferative and survival signals are enhanced and antimitotic inputs are either limited or lost. Co-repressor proteins significantly contribute with the disruption of these processes [16]. Therefore, understanding mechanisms involved in gene regulation and transcriptional network may lead to a better knowledge about the crucial functions of TFs, providing information to explore possibilities of their application as molecular targets in cancer therapy [18].

A valuable tool to study the transcription machinery is the DNA microarray technology [19], which measures the transcript expression of thousands of genes to identify changes in expression profiles at different biological conditions [20-24], thus allowing to compare different cell types under diverse treatment conditions. The influence of TP53 status on transcriptional profiles was previously described in tumor cell lines [25, 26]. Expression signatures of irradiated GBM cells were already performed for cell lines that are proficient and deficient for TP53 [27, 28].

Recently, information on the regulation of gene expression can also be used within the context of functional enrichment tests, and different databases containing TFs binding sites and other regulatory motifs are available, allowing to scan promoter regions of genes to detect the presence of target motifs [29]. This information allow to determining whether a set of pre-selected genes is under control of TFs. FatiGO + [30] is a web-based tool capable of associating TFs that are common to a gene set used as parameters. This TF prediction method was already applied to a GBM dataset obtained from public repositories of microarray experiments, and the up-regulation of two predicted TFs, E2F1 and E2F4, was validated for several GBM cell lines [31], demonstrating the suitability of this method.

In the current study, we aimed to identify TFs that could be predicted from significant differentially expressed genes (previously obtained in microarray experiments in irradiated GBM cells) using an in silico analysis.

We found few predicted TFs that were common between GBM cell lines, while several exclusive TFs were found for each cell line, indicating that the transcriptional response to ionizing radiation is very particular to each cell line examined in our microarray study, a fact that can be due to the genetic heterogeneity inherent to GBM cells. In spite of this, there was a convergence of biological functions among cell lines; the most relevant processes were related to apoptosis, cell proliferation, cell cycle, DNA repair, oxidative stress, among others. Furthermore, the present results also showed several TFs that were already reported as associated to cancer and stress responses.

2. Materials and methods

2.1. Briefly characterization of the experiment that provided the statistically modulated genes used for TF prediction

2.1.1. Cell culture and irradiation

Human GBM T98G and U87MG cell lines were supplied by the American Type Culture Collection (ATCC) (Rockville, Maryland, USA) and gently donated by Dr. Mari C. Sogayar (Universidade de São Paulo, Brazil). U343MG-a (U343), a cell line established from a primary malignant astrocytoma in an adult [32], was kindly donated by Dr. James T. Rutka (The Arthur and Sonia Labatt Brain Tumour Research Center, Canada); U251MG cell lines was also purchased from the ATCC (Rockville, MD, USA) and gently donated by Dr. Guido Lenz (Universidade Federal do Rio Grande do Sul, Brazil) [33]. All cell lines grown in the presence of DMEM + HAM F10 medium (Sigma-Aldrich, St. Louis, USA) plus 10% fetal calf serum (Cultilab, Campinas, Brazil), and kept at 37°C and 5% CO2, until they reach semiconfluency. Cells were sub-cultured and 1x10⁶ cells were seeded in 25 cm² flasks, being incubated at 37°C for 48 h, and irradiated with 8 Gy of gamma-rays (⁶⁰Co source, dose rate of 2.0 Gy / min., Unit Gammatron S-80, Siemens, 1.25 MeV, HC-FMRP/USP).

2.1.2. cDNA microarrays method and analysis

Two experiments with irradiated and sham-irradiated GBM cells were carried out using a glass slide microarrays containing ~4300 clones of cDNA probe (in replicates) from the human IMAGE Consortium cDNA library [34]; kindly provided by Dr. Catherine Nguyen (INSERM-CNRS, Marseille, France)], and prepared according to the protocol described by Hegde et al [35]. Microarrays were spotted onto glass slides (Corning, Lowell, MA, USA) by using a Generation III Array Spotter (Amersham Molecular Dynamics, Sunnyvale, USA) according to the manufacturer’s instructions.

Total RNA extraction was performed for all cell lines, 30 min. and 6 h after irradiation, using the Trizol reagent (Invitrogen, Carlsbad, USA) according to manufacturer’s instructions. Each cDNA sample was spotted twice in the slide (duplicate spots). The cDNA complex probes were prepared using the CyScribe Post Labeling Kit (Amersham Biosciences, Buckinghamshire, UK) as previously described [23]. Hybridizations were carried out using an automatic system (Automatic Slide Processor, Amersham Biosciences, UK) and signals were immediately captured after the final wash procedure, using a Generation III laser scanner (Amersham Biosciences, UK). This array platform was already used in several studies [22-26, 36].

2.1.3. Data acquisition and gene expression analysis

The provided microarray data was filtered and normalized [25, 36]. Following the normalization procedure, microarray data was exported to tab-delimited tables in MEV format and analyzed in MEV (v. 3.1) software [37].

The gene set submitted to SAM (Significance Analysis of Microarray [20]) were previously obtained by a t-test (α=5%) comparing irradiated (8 Gy) versus unirradiated (controls) T98G, U251MG, U343MG-a and U87MG cell lines, separately, considering two time points (30 min. and 6 h). The overall results are displayed in Table 1. The complete gene lists are available at http//www.rge.fmrp.usp.br/passos/genesgbm01/

Condition	Number of genes		Fold Change variation
up-regulated	down-regulated
U343MG-a (30 min.)	7	116	+1.53 to -2.42
U343MG-a (6 h)	3	11	+1.83 to -1.39
U87MG (30 min.)	56	73	+1.88 to - 2.95
U87MG (6 h)	86	54	+1.68 to -1.95
T98G (30 min)	32	0	+2.26 to +1.13
T98G (6 h)	16	7	+2.70 to -1.63
U251MG (30 min.)	12	69	+1.85 to -1.40
U251MG (6 h)	17	20	+2.28 to -2.18

Table 1.

Overall quantitative results on significant differentially expressed genes obtained by the DNA microarray method, and analysis performed by SAM – Significance Analysis of Microarray (FDR < 5 %), for the comparison irradiated versus un-irradiated cells.RNA samples from U87, U343, T98 and U251 cells were collected at 30 min. and 6 h following irradiation with 8 Gy of gamma-rays. Fold-change (+) or (-) means up- and down-regulation in transcript expression, respectively.

The list of significantly modulated genes was obtained for a FDR < 5%. U343 cells showed 123 and 14 significantly differentially expressed genes at 30 min and 6 h after irradiation, respectively, whereas U87 showed 129 genes at 30 min and 140 genes at 6 h; T98G cell line displayed 32 and 23 significantly up-regulated genes at 30 min. and 6 h, respectively, whereas U251 showed 81 genes at 30 min. and 37 genes at 6 h (Table 2).

2.2. Transcriptional factor analysis

The analysis of TFs related to the significant differentially expressed genes (SAM) was performed by applying the FatiGO + [30]. This program uses the TRANSFAC [38], and CisRed [39] transcription factors database, including their respective binding sites and regulated genes.

FatiGO + analyzes if the pre –selected set of genes (provided after SAM analysis), are under control of the same TF, and search for significant enrichments to each TF that is associated to the gene list compared to the complete reference list, containing ~4300 clones that were spotted onto the microarray slide [29].

The p-values obtained in the analysis of regulatory elements have been established by the program using the Fisher's exact test for multiple comparisons (unadjusted p-value). The Enrichment Index (EI) calculated for each TF corresponds to the increment obtained regarding the number of genes (%) statistically modulated (SAM) that are associated to a specific TF (List #1) divided by the total number of genes (%) in the array set that were predicted as targets for the same TF (List #2):

EI = % gene List #1/ % gene list #2

The TFs were selected according to unadjusted p-values < 0.05. The genes were submitted to FatiGO + v3.2, using the Gene symbol identifier and the selected gene distance of 10 kb. After selecting the TFs associated to modulated genes (SAM), a search was conducted in PubMed (http://www.ncbi.nlm. Nih.gov / sites / entrez /) looking for biological functions of those TFs.

2.3. Quantitative real-time PCR (qPCR)

We analyzed the transcript expression of HEB, a predict TF that was found associated to 57.7% of up-regulated genes in U87 cells, 30 min. after IR. The reverse transcription step was carried out in the remaining RNA samples from microarray experiments, with the Superscript III Reverse Transcriptase kit (Invitrogen, USA), according to manufacturer’s instructions. The integrity of cDNA samples was validated by the amplification of the endogenous B2M gene and visualization in agarose gel electrophoresis. qPCR was carried out using SYBR green master mix (Applied Biosystems, Foster City, USA) and the expression levels were estimated by the Relative Expression Software Tool (REST) [49], using 10000 interactions as setup parameter. All primers (Integrated DNA Technologies, Coralville, USA) were designed in Primer3 software [50] and are displayed on Table 2. The reactions were carried out in the Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, USA) equipment, using primer sets with an annealing temperature near 60°C and an amplicon of 100–120 bp. The PCR cycle was the following: pre-heating at 50°C for 2 min., 10 min. at 95°C (denaturation step), followed by 40 cycles at 95°C for 15 sec., and at 60°C for 60 sec. The dissociation curves were set up as following: 95°C for 15 sec., 60°C for 20 sec. and 95°C for 15 sec.

Primer	Sequence	PCR product size (pb)
B2M – forward	5’- AGGCTATCCAGCGTACTCCA - 3’	112
B2M – reverse	5’ - TCAATGTCGGATGGATGAAA - 3’
HEB – forward	5’ - CCGCTTGAGTTATCCTCCAC - 3’	116
HEB – reverse	5’ - GTGAGGCAGCAACGTAAGGT - 3’

Table 2.

Primer sequences used in Real Time qPCR; the housekeeping B2M gene was used as internal control.

2.4. Western Blot (WB)

Protein extraction was performed with the Trizol reagent (Invitrogen, Carlsbad-USA) according to the manufacturer's instructions, using the same samples for RNA extraction. These samples were obtained from U87 cells collected at 30 min. post-irradiation. The expression of HEB was analyzed by Western blot, using ACTB as internal control. Samples were prepared with 30 μg of total protein. After electrophoresis, proteins were transferred from the gel to the membrane Invitrolon PVDF using the XCell IITM Blot Module system (Invitrogen, Carlsbad - USA). The immunodetection and protein visualization were conducted with the WesternBreeze Chromogenic kit (Invitrogen, Carlsbad - USA). The antibodies used in this study were anti-HEB (Santa Cruz, Santa Cruz, USA), and anti-ACTB (Cell Signaling, Danvers, USA), dilution of 1:1000.

We performed densitometric analysis of WB bands using the GelPro Analyzer (MediaCybernetics, Rockville, USA) 4.0, and the relative expression of HEB was calculated relatively to ACTB.

3. Results

In the FatiGO + analysis, the lists of statistically modulated genes (SAM) were up-loaded in order to find TFs that were significantly associated with up-regulated and down-regulated genes for non-adjusted p-values < 0.05 (Table 3).

A Venn diagram was constructed based on the numbers of predicted TFs from data set previously obtained for each cell line (microarray experiments) (Fig. 1). TFs predicted for 30 min and 6 h were pooled together. Each cell line showed a number of exclusive TFs, but we also observed common TFs between cell lines. Out of 18 exclusive TFs found for U87MG cell line, PEBP (p = 0.008), Bach2 (p = 0.007), Freac-4 (p = 0.003), HLV (p = 0.006), Evi-1 (p = 0.009) displayed the lowest p-values, while PPARG and SEF-1 displayed the highest EI (31.3). U343 presented 9 exclusive TFs; High values of EI were found for MAF (33.1), E2F:DP-1 (22.0), PR (45.5) and STAT3 (38.5), and ARP-1 was the TF presenting the lowest p-value (0.009). T98G cells displayed only 6 exclusive TFs: EBF, Pax, Pbx1b, C/EBP, Poly A downstream element and Pax-9; two of them, EBF and PolyA showed low p-values, 0.005 and 0.007, respectively. Regarding U251MG cells, 13 TFs were predicted, and only APOLYA presented a high EI (27.8) (Table 3).

Interestingly, STAT3 was the common TF found for TP53 wild-type cells; however, this TF was associated with up-regulated genes in U87, and with down-regulated genes in U343. Only one TF (VBP) was common among three cell lines (U87, U343 and U251), being associated with down-regulated genes. Among the TP53 mutant cell lines, ATF4 was common between T98 and U251, associated with up-regulated genes (30 min. and 6 h). Two TFs were found common between U343 and U251, TEF (associated to down-regulated genes, 30 min.) and MAF (associated to up-regulated, 6 h). Finally, C/BPGamma was commonly predicted for up-regulated genes in U87 (30 min.) and U251 (6 h) cell lines (Fig. 1).

Therefore, our results showed that most of the predicted TFs were exclusive to each cell line and few TFs were common among the GBM cell lines; these results indicate that the transcriptional response to ionizing radiation is very particular to each cell line, and most probably this can be due to the genetic heterogeneity of GBM cells.

By using the real time qPCR method, we confirmed the expression of HEB to validate the in silico prediction for this TF. By using the REST 2009 software, we found that HEB was statistically up-regulated (+2.6) when comparing irradiated and sham-irradiated U87 cell lines (30 min.) (Fig. 2A). Primer efficiency was also determined for B2M (0.9615) and HEB (0.9652).

We also look for HEB protein expression by Western Blot; both ACTB and HEB antibodies were used for irradiated and sham-irrradiated U87, 30 min after irradiation (Fig.2B). The relative expression values calculated by densitometric analysis showed that HEB expression was 1.7 higher in irradiated (8 Gy) cells, relatively to the control value (Fig. 2C).

U343MG-a
Collection time	Transcription factor	% of genes (List 1)	% of genes (List 2)	EI	p-value
30 min. (↓)	ARP-1	14.8	6.7	2.2	0.009
TEF	37.5	25.3	1.5	0.013
VBP	19.3	10.4	1.9	0.013
Imperfect Hogness/Goldberg BOX	2.3	0.2	14.2	0.016
Muscle initiator sequence-20	20.5	12.2	1.7	0.031
Elk-1	44.3	33.1	1.3	0.038
Sox-5	4.6	1.5	3.1	0.048
ACAAT	12.5	6.7	1.9	0.049
6 h (↑)	MAF	50.0	1.5	33.1	0.031
E2F-4:DP-1	50.0	2.3	22.0	0.046
6 h (↓)	ICSBP	40.0	3.7	10.9	0.013
PR	20.0	0.4	45.5	0.024
STAT3	20.0	0.5	38.5	0.028
ARP-1	40.0	6.7	6.0	0.040
U87 MG
Collection time	Transcription factor	% of genes (List 1)	% of genes (List 2)	EI	p-value
30 min. (↑)	C/EBPgamma	85.0	64.4	1.3	0.007
AP-1	80.0	61.1	1.3	0.014
HEB	57.5	38.5	1.5	0.021
SREBP-1	97.5	84.5	1.2	0.024
FOXP3	87.5	72.0	1.2	0.032
PPARG	2.5	0.1	31.3	0.046
SEF-1	2.5	0.1	31.3	0.046
30 min. (↓)	Bach2	63.5	43.9	1.4	0.007
PEBP	28.9	14.4	2.0	0.008
COUP-TF:HNF-4	11.5	3.7	3.1	0.014
MEF-3	5.8	1.2	5.0	0.026
FOX	78.9	64.8	1.2	0.039
6 h (↑)	DEC	45.6	31.9	1.4	0.032
STAT3	3.5	0.5	6.8	0.042
6 h (↓)	Freac-4	10.5	1.6	6.8	0.003
HLF	50.0	28.2	1.8	0.006
Evi-1	97.4	81.6	1.2	0.009
VBP	23.7	10.4	2.3	0.015
TCF-4	55.3	37.2	1.5	0.028
AP-1	79.0	61.1	1.3	0.028
Gfi-1	10.5	3.1	3.4	0.031
CRE-BP1	15.8	6.6	2.4	0.039
HNF-4alpha	29.0	16.4	1.8	0.047
T98G
Collection time	Transcription factor	% of genes (List 1)	% of genes (List 2)	EI	p-value
30 min. (↑)	EBF	29.6	10.5	2.8	0.005
ATF4	22.2	8.4	2.6	0.023
Pax	92.6	75.3	1.2	0.041
6 h (↑)	Pbx1b	14.3	1.8	8.2	0.026
6 h (↓)	C/EBP	100.0	20.0	5.0	0.002
Poly A downstream element	75.0	12.6	6.0	0.007
Pax-9	75.0	23.8	3.2	0.045
U251MG
Collection time	Transcription factor	% of genes (List 1)	% of genes (List 2)	EI	p-value
30 min. (↑)	SMAD-4	40.0	10.4	3.9	0.015
PTF1-beta	20.0	2.7	7.4	0.030
APOLYA	10.0	0.4	27.8	0.039
30 min. (↓)	VBP	21.7	10.4	2.1	0.025
HNF-6	6.5	1.4	4.7	0.030
E2F	52.2	36.0	1.4	0.030
TEF	39.1	25.3	1.5	0.040
TTF1	23.9	13.3	1.8	0.047
CDP CR1	54.4	39.7	1.4	0.049
POU1F1	54.4	39.7	1.4	0.049
6 h (↑)	MAF	15.4	1.5	10.2	0.017
CREB	46.2	17.8	2.6	0.018
ATF4	30.8	8.4	3.7	0.020
MEIS1B:HOXA9	15.4	1.9	8.2	0.025
C/EBPgamma	92.3	64.4	1.4	0.041
6 h (↓)	HES1	45.5	15.6	2.9	0.019
Lmo2 complex	27.3	6.0	4.5	0.026
ATATA	18.2	2.9	6.2	0.041

Table 3.

Transcription factors associated with statistically modulated genes (SAM, FDR ≤ 5 %), as predicted by the FATIGO + v3.2., analysis performed for U343MG-a, U87MG, T98G and U251MG cell lines (30 min. and 6 h post-irradiation). We used gene lists that showed patterns of repression (↓) and induction (↑) in irradiated cells compared with mock-irradiated. The Enrichment Index (EI) calculated for each TF corresponds to the increment regarding the number of genes (%) statistically modulated (SAM) that are associated to a specific TF (List #1) divided by the total number of genes (%) in the array set that were predicted as targets for the same TF (List #2). The gene distance for the analysis of the TFs was 10 kb.

Figure 1.

Venn diagram showing predicted TFs associated with significant differentially expressed genes (from microarray experiments) selected for four GBM cell lines, comparing irradiated versus sham-irradiated cells, collected at 30 min. and 6 h following irradiation. TF prediction was carried out using FatiGO + v3.2.

Figure 2.

HEB expression. A) HEB expression levels obtained by the qPCR method. This TF was found associated with up-regulated genes in U87 cells, 30 min. after irradiation. Boxes represent the interquartile range; the dotted line represents the median value; whiskers represent the minimum and maximum observations. B) Protein expression analyzed by Western Blot using antibodies for HEB (Santa Cruz) and ACTB (Cell Signalling) as endogenous control. C) Densitometric analysis of Western Blot bands using the Gel Pro Analyzer 4.0 software displayed for HEB expression relatively to ACTB.

4. Discussion

Recently, genome wide technologies, such as DNA microarrays, provide a huge amount of information about gene expression, but require additional bioinformatics analyses for data interpretation. In order to reduce complex signatures to a small number of activated transcriptional elements, new bioinformatics tools have been developed. To date, genome-wide TF-binding regions and sites were identified using a variety of indirect methods and data sets, revealing abundant binding sites for different TFs in mammalian cells [40-43]. Using lists of differentially expressed genes that were generated by microarray experiments, it is possible to predict TFs that can target common binding sites to a gene set. In the current study, we performed an in silico analysis (FatiGO + v3.2.) to identify TFs from a list of significant differentially expressed genes selected for irradiated GBM cell lines in microarray experiments. Only few predicted TFs were common to GBM cell lines, while several TFs were exclusive to each cell line, indicating that the transcriptional response to ionizing radiation is very peculiar to each cell line examined in our microarray study. The most relevant predicted TFs are discussed below. While few predicted TFs were shared between different cell lines, several TFs were found exclusive to each cell line, except U251.

4.1. Commonly predicted TFs for two or three GBM cell lines

We found few TFs that were predicted for more than one cell line: MAF, TEF, ATF4, STAT3, VBP and C/EBPGamma. Most of these TFs (STAT3, TEF and VBP) are related to apoptosis, while other biological classes were also found, such as oxidative stress (ATF4), differentiation (MAF) and nucleotide excision repair (C/EBPGamma).

STAT3, signal transducer and activator of transcription 3 (acute-phase response factor), is part of the STAT family of cytoplasmic latent transcription factors, and was predicted for the TP53 wild type cells, U87 (up-regulated genes, 6 h) and U343 (down-regulated genes, 6 h). Phosphorylated STAT3 leads to transcriptional activation of downstream genes involved in processes such as cell proliferation, suppression of apoptosis, and angiogenesis [44, 45]. It was demonstrated that STAT3 is constitutively activated and overexpressed in human gliomas; STAT3 activation correlates with malignancy [46, 47], while STAT3 inhibition reduces the lethality of GBM tumors in vivo [48], and its inhibition have been tested in phase 0 trial in head and neck cancers [49].

Thyrotroph embryonic factor (TEF) and human hepatic leukemia factor (HLF) are members of the PAR (proline and acidic amino acid-rich) subfamily of basic region/leucine zipper (bZIP) transcription factors. The chicken vitellogenin gene-binding protein (VBP) is also a bZIP TF member and is considered as the chicken homologue of TEF. TEF was predicted from down-regulated genes in U343 and U251 (30 min.), and its homolog, VBP, from down-regulated genes in U87 (6 h), U343 (30 min.) and U251 (30 min.). PAR bZIP proteins have recently been shown to be involved in amino acid and neurotransmitter metabolism in both liver and brain [50]. PAR bZIP proteins are also able to transactivate the promoter of bcl-gS which is directly involved in apoptosis induction. Consistently, transfection of TEF induces the expression of endogenous bcl-gS in cancer cells, independently on TP53 [51].

Activating transcription factor 4 (tax-responsive enhancer element B67); activating transcription factor 4C (ATF4) belongs to the large ATF/CREB family of transcription factors [52] and was predicted from up-regulated genes in T98 (30 min) and U251 (6 h). Up-regulation of ATF4 is directly involved in the endoplasmic reticulum (ER) stress through induction of CHOP in GBM treated with Nelfinavir (protease inhibitor class of drugs) [53] or in concert with PERK, GADD34 and EIF2alpha in Hela cells submitted to hypoxia [54]. Therefore, activation of ATF4 was already reported in GBM treated cells.

MAF, the v-maf musculoaponeurotic fibrosarcoma oncogene homolog (avian) is a unique subclass of bZIP proteins and was predicted from up-regulated genes in U343 and U251 at 6 h. Depending on the binding site and binding partner, the encoded protein can be a transcriptional activator or repressor. Members of the MAF family appear to play important roles in the regulation of differentiation [55]. MAF was found up-regulated in various cancers, such as colon cancer (but only in tumors that presented high levels of COX-2 expression) [56], a small subset of myelomas, hairy cell leukemia, T- and NK-cell neoplasms and small cell lymphomas [57].

C/EBPGamma, a member of the CCAAT/enhancer-binding protein (C/EBP) family of transcription factors was predicted from up-regulated genes in U87 (30 min.) and U251 (6 h). This TF regulates the expression of ERCC5 [58], and is a participant of DNA repair [59], particularly in the nucleotide excision repair [60].

The predicted TFs represent the overall GBM response to irradiation, since they were selected for more than one cell line, and as mentioned above, their functions are directly associated with stress responses involving apoptosis, DNA repair and ER stress. Moreover, as an example of STAT3, which is in clinical trial [49], predicted TFs may constitute potential targets to be investigated and validated in cancer treatment.

4.2. Exclusively predicted TFs for irradiated U343 cell line

Few predicted TFs (PR, E2F4:DP-1, and ARP-1) associated with statistically significant modulated genes in irradiated U343 cells were found exclusive to this cell line. The functions of these TFs are mainly associated with cell cycle/ tumor growth, being involved in various types of cancer. The overexpression of E2F4 and its binding partner DP-1 revealed a dual function of E2F4, which acts as an activator as well as a repressor, being implicated in positive regulation of the cell cycle [61]. In a previous work, the up-regulation of E2F4 was confirmed for several GBM cell lines [31], demonstrating the potential of this TF as molecular target in cancer therapy.

Progesterone receptor (PR), a nuclear receptor transcription factor was associated with down-regulated genes in U343 cells (6 h). Steroid hormones participate in several physiological and pathological processes in the brain, including the regulation of tumor cell growth [62-64]. Progesterone exerts many of its effects by the interaction with specific intracellular receptors [62, 65].

ARP-1, also known as orphan nuclear receptor chicken ovalbumin upstream promoter transcription factor II (COUP-TFII) is a member of the steroid/thyroid nuclear hormone receptor superfamily [66] and can act as transcriptional repressor or activator. ARP-1 plays critical roles in organogenesis [67-70], and is a major angiogenesis regulator within the tumor microenvironment during pancreatic tumor progression and metastasis [71]. Besides, ARP-1 was associated with therapy response in oligodendroglial tumors with 1p/19q loss [72]. In the present study, this TF was predicted from down-regulated genes (30 min. and 6 h), indicating its possible involvement in radiation responses restricted to U343 cell line.

4.3. Exclusively predicted TFs for irradiated U87 cell line

A high number of TFs with low p-value (Bach2, PEBP, Freac-4, HLF and Evi-1) and high EI (PPARG and SEF-1) was associated with significant expressed genes in U87 cell line.

PPARG is a member of the peroxisome proliferator activated receptor (PPAR) family, a subfamily of the nuclear receptor superfamily [73]. The protein level of this receptor has been recently identified as a significant prognostic marker [74]. Interestingly, recent studies have shown that PPARG is expressed in normal and malignant human brain, and the treatment with PPARG agonists induces growth arrest and apoptosis in brain tumor cells in vitro and in animal models in vivo [75-77]. Recent findings show that PPARG agonists regulate growth and expansion of brain tumor stem cells [78] and also altered the expression of stemness genes [79]. Unfortunately, clinical trials also failed to demonstrate the effectiveness of such agonists as a monotherapy for cancer treatment, a fact which stimulates the search for combination treatments to enhance their effects [80].

Basic leucine zipper transcription factor 2 (Bach2) is an evolutionarily related member of the BTB-basic region leucine zipper transcription factor family. Bach2 can function as transcriptional activator and repressor [81]. This TF down-regulates cell proliferation of the neuroblastoma cell line N1E-115 and negatively affects their potential to differentiate, being considered as gatekeeper of the differentiated status [82]. Bach2 presents high frequency of loss of heterozygosity of the Bach2 gene in human B cell lymphomas [83]. Consistent with its putative role as a tumor suppressor, Bach2 enhances apoptosis in response to oxidative stress [84, 85].

The transcription Factor PEBP, also called Raf kinase inhibitor protein (RKIP) is a member of the phosphatidylethanolamine-binding protein (PEBP) family. RKIP plays a pivotal modulatory role in several protein kinase signaling cascades. RKIP regulates the activity of the Raf/MEK/ERK, which is responsible for proliferation and differentiation of diverse cell types [86]. It has been reported that RKIP was poorly expressed in primary tumors, being absent in various metastatic cancers; its induction sensitize resistant tumor cells to apoptosis by various chemo- and immunotherapeutic drugs, as well as inhibitors of metastasis [87].The absence of RKIP is also associated with highly malignant behavior and poor survival of patients [88].

The forkhead domain is a monomeric DNA binding motif that defines a rapidly growing family of eukaryotic transcriptional regulators. We found Freac-4, also known as Forkhead Box D1 (FOXD1) associated with down-regulated genes in U87 cells (6 h). This gene was found repressed in chemoresistant tumors, as analyzed by microarrays [89]. However, FOXD1 and FOXD2 were highly expressed in prostate cancer and lymph node metastases, among various cancer types [90]. In another study, using kidney-derived cell lines, it was suggested that FOXD1 may be regulated by TP53, WTAR (a mutated form of WT1) and WT1 [91].

As already mentioned, HLF is a member of PAR bZIP transcription factors, and was associated with down-regulated genes, 6 h after irradiation in U87. PAR bZIP proteins are also involved in apoptosis induction [51].The fused gene E2A-HLF was responsible for the development of lymphoid malignancies in 60 % of the transgenic mice [92].

Activator protein one (AP1) transcription factors are a family of jun and fos proteins, whose subunits present diverse pro/anti-cancer effects, like inhibition or increase in proliferation, inhibition of apoptosis and angiogenesis [93, 94]. AP-1 is one of the genes early activated after radiation in primary human B cells [95]. The inhibition of AP-1 blocks the proliferation of breast tumor cells by suppressing the growth factor signaling [96]. The modulation of AP-1 activity may be a new attempt to reduce the malignant transformation. However, only the function involved with malignancy should be targeted [97], since AP-1 presents oncogenic and anti-oncogenic properties. This TF was associated to U87 cells (up-regulated genes, 30 min.) and (down-regulated genes, 6 h).

The EVI1 gene encodes a zinc finger transcription factor with important roles in normal development and leukemogenesis. Reports in animal model and findings in in vitro studies. showed that EVI1 affected cellular proliferation, differentiation, and apoptosis [98]. EVI-1 was also found overexpressed in infratentorial ependymomas, it can promote proliferation of ependymal tumor cells, and its expression indicates an unfavorable prognosis [99].

U87 cell line presented several predicted TFs with significant p values and higher enrichment index than other cell lines. Most of the predicted TFs are related to apoptosis (PPARG, Batch2, PEBP, HLF, AP1 and EVI1), but they were associated with up- or down-regulated genes. Overall, the biological functions of these TFs are related to cell proliferation, differentiation and tumor growth, indicating the relevance of their deregulation in cancer development and malignancy, and possibly, in tumor responses to anti-cancer therapies.

4.4. Exclusively predicted TFs for irradiated T98G cell line

Only two TFs were predicted for T98G cells: EBF and C/EBP. The early B-cell factors (EBF) are a family of four highly conserved DNA-binding transcription factors with an atypical zinc-finger and helix-loop-helix motif. Zardo and colleagues found that the EBF3 locus on the human chromosome 10q is deleted or methylated in brain tumors [100]. Functional studies revealed that EBF3 activates genes involved in cell cycle arrest and apoptosis, while in opposite, it represses genes involved in cell survival and proliferation [101]. Therefore, EBFs represent a novel tumor suppressor whose inactivation blocks normal development and contributes to tumorigenesis of diverse types of human cancer [102].

CEBP is also known as basic leucine zipper transcription factor, CCAAT/enhancer binding protein alpha (CEBPA), which directly interacts with CDK2 and CDK4 and arrests cell proliferation by inhibition of these kinases [103]. CEBPA is crucial for normal granulopoiesis, and dominant-negative mutations of CEBPA gene were found in patients with myeloblastic subtypes (M1 and M2) of acute myeloid leukemia [104]. CEBPA also plays a role in DNA damage response dependent on TP53, as observed in keratinocytes [105]. C/EBPA was found silenced in human squamous cell carcinoma (SCC), and loss of C/EBPA confers susceptibility to UVB-induced skin SCCs involving defective cell cycle arrest in response to UVB [106]. Interestingly, these findings indicate the role of CEBP in DNA damage responses, and possibly, the potential of this TF to be explored as therapeutic molecular target.

4.5. Validation of TF prediction

As a predicted TF associated with up-regulated genes, HEB (p-value = 0.021 and EI = 1.5) was chosen to be studied in terms of expression levels, aiming to validate the in silico analysis, although for a single TF. Interestingly, we showed that HEB transcript expression was up-regulated (+2.6) in irradiated U87 cell line, 30 min. after irradiation, while HEB protein expression analyzed by Western blot was 1.7 higher in irradiated (8 Gy) cells, relatively to un-exposed controls.

HEB is a member of the class A basic helix-loop-helix (bHLH) family that participates in the nervous system development [107, 108]. According to O'Neil et al. [109], the repression of E47/HEB has been associated with the induction of leukemia in mice. In another study, it was demonstrated the induction of HEB in gliomas compared with non-neoplastic brain tissue [110]. Moreover, HEB seems to be involved in cell proliferation control of neural stem cells and also progenitor cells, being important to sustain their undifferentiated state during embryonic and adult neurogenesis [108]. Although HEB expression has not yet been correlated with radiation responses in GBM cells, in the present study, we found its association with significant differentially expressed genes at 30 min. following irradiation in U87 cell line; interestingly, we also showed that HEB transcript and protein expression was induced in irradiated U87 cell line, 30 min. after irradiation. This finding, although restrict to one TF, indicates the validity of the TF prediction by in silico analysis.

5. Conclusions: Lessons from TF prediction in irradiated GBM cells

The present findings about prediction of TFs associated to differentially expressed genes in GBM cell lines showed that few TFs were shared among different GBM cell lines, while several TFs were found exclusive to each cell line, indicating that the transcriptional response to ionizing radiation is very particular to each cell line, probably due to the genetic heterogeneity, which is characteristic of GBM cell lines. In spite of this observation, several biological functions were similar among cell lines, such as apoptosis, cell proliferation, cell cycle control, DNA repair, ER stress, and differentiation. Furthermore, most of the predicted TFs were already reported as differentially expressed, deleted or mutated in cancer, including GBM. However, apart the similarity of biological functions, different pathways seems to be associated to the predicted TFs. Interestingly, we could not find TP53 as a TF associated to the data set (List #1) analyzed in the present study, even for the GBM cells that were wild-type for TP53 gene, and even considering the presence of TP53 cDNA clone in the microarray slide. It is possible that the TP53 protein could not be activated in GBM cells, impairing its action as transcription factor, as previously suggested by other authors [111].

The most intriguing finding refers to apoptotic related TFs. Probably, predicted TFs related to apoptosis control, and found associated with expressed genes at early time (30 minutes) following irradiation, are related to survival in GBM cells; this is supported by reports showing that in general, these cells are very resistant to undergo apoptosis, even under conditions of drug treatment or radiation exposure [112-114]. In fact, GBM cells seem to be capable of activating several pathways to escaping from cell killing by anticancer therapies.

Even considering the relevance of our findings, some methodological limitations should be mentioned regarding in silico prediction of TFs. Despite the great advancement in terms of DNA binding sites detection, it is hard to determine which sites are functional regulatory elements that influence transcription. It is possible that a considerable fraction of these binding sites are nonfunctional and may constitute biological noise [115]. Other choices, such as ChIP experiments, may overcome this concern by detecting indirect TF-DNA interactions through protein/protein interaction [116].

In spite of the limitation mentioned above, in a previous study, we validated the expression of E2F [31] and the HEB expression was confirmed in the present study, both of them in GBM cell lines. In addition, we selected predicted TFs that were associated with stress response genes, and importantly, the TFs were reported as deregulated or mutated in different cancer types, thus indicating the relevance of further studies to better exploring the role of TFs in the context of therapeutic strategies based on molecular targets.