Abstract
Chemical carcinogens may cause a multitude of effects inside cells, thereby affecting transcript levels of genes by direct activation of transcription factors (TF) or indirectly through the formation of DNA damage. As the temporal profiles of these responses may be profoundly different, examining time-dependent changes may provide new insights in TF networks related to cellular responses to chemical carcinogens. Therefore, we investigated in human hepatoma cells gene expression changes caused by benzo[a]pyrene at 12 time points after exposure, in relation to DNA adduct and cell cycle. Temporal profiles for functional gene sets demonstrate both early and late effects in up- and downregulation of relevant gene sets involved in cell cycle, apoptosis, DNA repair, and metabolism of amino acids and lipids. Many significant transcription regulation networks appeared to be around TF that are proto-oncogenes or tumor suppressor genes. The time series analysis tool Short Time-series Expression Miner (STEM) was used to identify time-dependent correlation of pathways, gene sets, TF networks, and biological parameters. Most correlations are with DNA adduct levels, which is an early response, and less with the later responses on G1 and S phase cells. The majority of the modulated genes in the Reactome pathways can be regulated by several of these TF, e.g., 73% by nuclear factor-kappa B and 34–42% by c-MYC, SRF, AP1, and E2F1. All these TF can also regulate one or more of the others. Our data indicate that a complex network of a few TF is responsible for the majority of the transcriptional changes induced by BaP. This network hardly changes over time, despite that the transcriptional profiles clearly alter, suggesting that also other regulatory mechanisms are involved.
Chemicals can cause cancer through a variety of mechanisms. Genotoxic carcinogens induce damage to the genetic material, either directly by covalently binding to DNA or indirectly by interfering with the mitotic machinery, ultimately leading to mutations in genes or large aberrations in chromosomes (Hayashi, 1992). If tumor suppressor genes are inactivated or proto-oncogenes are activated by this process, this may lead to cancer formation. The mode of action of nongenotoxic carcinogens does not imply damage to DNA or chromosomes and can be very diverse, ranging from stimulation of cell proliferation, suppression of apoptosis, and induction of oxidative stress to activation of biotransformation enzymes, suppression of the immune system, and so on (Shaw and Jones, 1994).
Polycyclic aromatic hydrocarbons (PAHs) comprise a large group of structurally related compounds, including many carcinogens, to which humans are daily exposed through the environment and food. PAHs are well-known carcinogens, especially for lung following exposure by inhalation, but also for other tissues such as the liver because of oral exposure. Many carcinogenic PAHs, such as benzo[a]pyrene (BaP), share both genotoxic and nongenotoxic properties. The genotoxic properties are because of biotransformation of BaP to DNA reactive metabolites and the generation of reactive oxygen species (ROS) (Balu et al., 2004; Caino et al., 2007; Cavalieri and Rogan, 1995; Cheng et al., 1989; Dipple et al., 1999; Khan and Anderson, 2001; Penning et al., 1999).
The nongenotoxic properties of PAH are assumed to act through activating the transcription factor AhR, the aromatic hydrocarbon receptor, which thereupon dimerizes with the AhR nuclear translocator (ARNT) and induces the transcription of its target genes, including phase I enzymes like CYP1A1, CYP1A2, and CYP1B1 and several phase II enzymes (Ma, 2001; Nebert et al., 2004). In addition, PAH-mediated oxidative stress leads via the activation of transcription factor Nrf2 to the induction of multiple phase I and II biotransformation enzymes, including aldo-keto reductases, glutathione S-transferase-P, NAD(P)H:quinone oxireductases, and UDP-glucuronosyl transferase 1A6 (Wang et al., 2007).
In order to protect the cells from the deleterious effects of DNA damage, cells possess many protective mechanisms. DNA repair can be by nucleotide excision repair (NER) for the PAH-DNA adducts and base excision repair (BER) for the oxidized damages (Braithwaite et al., 1998; Izumi et al., 2003). In order to allow for more time for repair, cell proliferation can be arrested at several phases of the cell cycle. Instead of this cell cycle arrest, apoptosis can also be induced, especially at high DNA damage levels (Bartek and Lukas, 2007). In these DNA damage response processes, many stress response, proto-oncogenes, and tumor suppressor genes are involved, such as nuclear factor-kappa B (NF-κB), AP-1, MYC, and p53 (Bolotina et al., 2007; Croce, 2008; Meek, 2009; Sherr, 2004).
Exposure of cells to a full carcinogen such as BaP induces multiple molecular responses through the activities of multiple transcription factors (TF) and the temporal profiles of these responses may be profoundly different. So far, however, their time dependencies and possible functional interactions have not been examined and are thus not well understood. Examining transcriptional changes over many time points and several days may provide new insights in TF networks related to cellular responses to chemical carcinogens.
The objective of our study therefore was to investigate time-dependent changes in cellular responses—both phenotypic markers and gene expression profiles—upon exposure to a full carcinogen, thereby using BaP as a model compound, and to relate these to possible TF networks. This is achieved by extensive time series analyses of global gene expression changes in human hepatoma cells (HepG2) in relation to phenotypic events, such as DNA adduct formation and cell cycle progression. Analyses are focused on functional gene sets, TF networks, and the application of time series tools. HepG2 cells are chosen as in vitro model for several reasons. They represent the liver as target organ for PAH-induced hepatocarcinogenesis following oral exposure, are metabolically competent with respect to biotransformation of mutagens and carcinogens, carry no p53 mutations, are frequently used in toxicology and gene expression studies, and are known to fully respond to the mutagenic and epigenetic effects caused by BaP (Hsu et al., 1993; Knasmuller et al., 1998; Staal et al., 2006; Westerink and Schoonen, 2007a, 2007b; Westerink et al., 2010; Wilkening et al., 2003).
MATERIALS AND METHODS
Cell Culture and Treatment
HepG2 cells were cultured and treated with 3μM BaP (purity 97%, CAS no. 50-32-8; Sigma-Aldrich, Zwijndrecht, the Netherlands) or vehicle control (dimethyl sulfoxide, 0.1%) as described previously (Staal et al., 2007). This is the lowest BaP dose with a maximum accumulation of cells in S phase, higher doses do not increase that but only extend the period. No other doses were investigated, as based on our previous studies, we knew that time has a much larger effect on gene expression profiles than dose (Staal et al., 2006, 2007). The cells were exposed for 3, 6, 9, 12, 15, 18, 24, 30, 36, 48, 54, or 60 h (∼90% confluency after 24 h and ∼95% after 60 h), whereafter they were either fixed with cold methanol and stored at −20°C for FACS analysis or media was removed from the culture flasks and Trizol (Gibco/BRL, Breda, The Netherlands) was added for RNA and DNA isolation. HepG2 cells have a doubling time of about 30 h and grow in clumps because of which they hardly reach full confluency. Two independent experiments were conducted.
BaP in Culture Medium
BaP levels in culture medium were determined following extraction with hexane, evaporation of hexane, dissolving the residue in acetonitrile, and analyses on an HPLC gradient system (Supelcosil LC-PAH [column 25 cm × 4.6 mm, 5 μm]; solvent gradient at 1.5 ml/min: 3 min 64% acetonitril—a linear gradient to 100% acetonitril in 20 min—a 15-min hold at this) with fluorescence detection (ex/em at 280/> 350nm).
γH2AX Foci Analysis
HepG2 cells cultured and exposed on coverslips were analyzed for DNA double-strand breaks by immunofluorescence staining of γH2AX foci as described (Rothkamm and Lobrich, 2003). Following fixation, foci were stained with the anti-phospho-Histone H2A-X (ser139) Clone JBW 301 antibody (Upstate Biotechnology, Lake Placid, NY) and Alexa fluor 488 goat anti-mouse IgG antibody (Invitrogen, Breda, The Netherlands). Cells were scored manually using a damage score, ranging from 0 to 4, and sample scores (100 cells/sample) are given as percentage of the maximal damage.
Flow Cytometry Analyses for Cell Cycle and Apoptosis
Analyses of cell cycle profiles and apoptosis were performed as previously described using a FACSort (Becton Dickinson, Sunnyvale, CA) (Staal et al., 2007). For DNA content, cells were stained with propidium iodide; apoptotic cells were visualized with the primary antibody M30 CytoDeath (Roche, Penzberg, Germany) and fluorescein isothiocyanate-conjugated anti-mouse Ig as secondary antibody (DakoCytomation, Glostrup, Denmark). For each sample, 10,000 cells were analyzed. Cells in the G0-1, S or G2-M phase were expressed as a percentage of the total number of cells.
RNA and DNA Isolation
RNA was isolated from the Trizol solutions according to the producer’s manual and purified with the RNeasy mini kit (Qiagen Westburg bv., Leusden, The Netherlands). After isolation, the remaining phases were used for DNA isolation according to manufacturer’s protocol. RNA and DNA quantity was measured on a spectrophotometer, and RNA quality was determined on a BioAnalyzer (Agilent Technologies, Breda, The Netherlands). Only RNA samples that showed clear 18S and 28S peaks and with a RIN level higher than 8 were used.
DNA Adduct Analyses
DNA adduct levels were determined by 32P-postlabelling according to the procedure originally described by Reddy and Randerath (1986) with modifications described by Godschalk et al. (1998).
Gene Expression Analyses with DNA Microarrays
Labeling and hybridization.
Labeling and hybridization of RNA samples was performed according to Agilent’s manual for microarrays (Agilent Technologies). Samples from BaP- or vehicle-treated cells were labeled with cyanine 3 (Cy3) or cyanine 5 (Cy5). Complementary RNA of the time matched–treated and control samples was applied on the G4110B Agilent 22K Human Oligo Microarray, hybridized, and washed according to Agilent’s manual, and slides were scanned on a ScanArrayExpress (Packard Biochip Technologies, Perkin Elmer Life Sciences, Boston, MA) with fixed laser power (100%) and photo multiplier tube gain (55% for Cy5 and 60% for Cy3). For each biological experiment, two hybridizations per time point were conducted, with swapped Cy3 and Cy5 dyes. In total, this resulted in 48 hybridizations.
Image analysis and processing.
The images (10 μm resolution; 16-bit tiff) were processed with ImaGene 6.0 to quantify spot signals (BioDiscovery Inc., Los Angeles, CA). Irregular spots were automatically flagged and excluded from data analysis. Data from ImaGene were further processed with GeneSight 4.1.6 (BioDiscovery Inc.). For each spot, background was subtracted and flagged spots as well as spots with a net expression level below 20 in both channels were omitted. Data were log base 2-transformed and locally weighted scatterplot smoothing normalized. If for a gene more probes were on the array, replicates were averaged while omitting outliers (> 2 SDs). Raw data are available at ArrayExpress (accession numbers E-TOXM-22 and E-TOXM-23 for the first and second experiment, respectively) (www.ebi.ac.uk/arrayexpress/).
Gene Expression Data Analyses and Data Mining
Differentially expressed genes.
First, genes were selected, for which all four replicate hybridizations gave an expression difference of more than 0.2 or all four less than −0.2 (meaning all replicates in the same direction; on a normal scale, this is a 15% increase or decrease of transcript level). This resulted in 5211 genes. Second, genes were selected by ANOVA using the plug-in “Time series analysis” from BRB Array Tools (version 3.4 and R version 2.2.0; http://linus.nci.nih.gov/BRB-ArrayTools.html), with “Time” as descriptor. This resulted in 2808 genes, with a false discovery rate < 0.1. All further analyses on significantly differentially expressed genes start with these 2808 genes.
Functional interpretation of significantly differentially expressed gene sets.
Reactome (www.reactome.org; version 25, November 2008) (Vastrik et al., 2007) is used to examine the biological processes that may be affected at the level of reactions and reaction pathways, a level of detail that cannot be obtained by gene ontology annotations, and KEGG and BioCarta pathways. Reactions are grouped into causal chains to form pathways. The tool Skypainter is used to determine per time point that events (reactions and/or reaction pathways) harbor differentially expressed genes (down- and upregulated genes separately), selecting only events with at least three genes. The proportion of these number of differential genes per event (time point/experiment) with total numbers present in the set of 2808 genes is calculated for each pathway.
Functional interpretation of all expressed genes.
T-profiler is used for functional annotation of pathways and processes in the complete data set of genes without preselecting modulated genes. T-profiler uses the unpaired t-test to score the difference between the mean expression level of predefined sets of genes and that of all other genes within the complete data set (Boorsma et al., 2005). Pathways and processes included gene sets from GO terms, KEGG pathways, GenMAPP pathways, Motifs, Broad/MIT pathways (manually curated), gene sets, and gene sets retrieved from literature, i.e., lists of genes modulated by BaP or TCDD in HepG2 cells (Frueh et al., 2001; Hockley et al., 2006; Kim et al., 2006) and a set of genes with known antioxidant response elements (AREs), also known as Nrf2-binding sites (Wang et al., 2007). Per array, t values were determined for each gene set. Only gene sets with an absolute t value > 1 at at least 2 succeeding time points in all four replicate arrays were selected for further analyses.
Network building around TF.
MetaCore (GeneGo, San Diego, CA) is used to examine connectivity between the genes, using their network building algorithms (Dezso et al., 2009; Ekins et al., 2007). To reduce complexity, only 233 genes were used, each with at least a twofold change in three out of the four replicates. Network were analyzed using the “Transcription regulation” algorithm, where the created networks are centered on transcriptional factors. The network is built dynamically, connecting all genes via direct mechanistic interactions based on manually curated literature evidence. It can then be used for local functional enrichment of the network itself and to associate the interconnected cluster of genes with specific processes or diseases.
Linking TF and modulated genes.
For the 30 TF around which significant transcription regulation networks were observed (see above), we examined if they could be involved in regulation of the affected genes mapped on the Reactome reaction pathways. Both an inventory of TF-binding sites in promoter regions of these genes, as retrieved from GEPAS, as well as literature-based knowledge on regulation, as retrieved from the ElDorado database from Genomatix (www.genomatix.com, v-02-2010; Werner, 2007, 2009), were used.
Time series analyses.
For identification of genes, reactions, pathways, or gene sets co-regulated time dependently and correlating with a biological parameter, the software tool “Short Time-series Expression Miner” (STEM, version v1.3.6; http://www.cs.cmu.edu/∼jernst/stem/) was used (Bar-Joseph et al., 2008; Ernst and Bar-Joseph, 2006). STEM finds statistically significant patterns from time series microarray experiments and can compare data sets across experiments. The clustering algorithm assigns each gene, reaction, pathway, or gene set to these profiles based on correlation coefficients (also between the two independent experiments) and permutation analyses (n = 50). For all parameters, log-transformed data were used, normalized such that treated-control ratios ranged from + 1 to −1, with 0 for no change.
RESULTS
In order to investigate the time-dependent changes in cellular responses upon exposure to BaP by bioinformatics tools dedicated to time series analyses, first, the phenotypic markers and then the gene expression profiles will be described.
Effects on BaP Levels, DNA Adduct Levels, Cell Cycle, and Apoptosis
Clearance by HepG2 cells of BaP from the culture medium was fast, with 1% remaining after 30–36 h, levels being not detectable after 48 h (see Supplementary DataSupplementary Data). BaP-DNA adduct levels increased between 3 and 12 h and were stable until 18/30 h, whereafter they slowly declined (Fig. 1). For the second experiment, adduct levels reached a lower maximum and faster decreased than for the first experiment. Adduct levels stabilized to 15–25% of the highest level, even after 60 h, suggesting a diminished repair of adducts during the last 12–30 h. Similar results were obtained by in situ staining for γH2AX foci, a marker for double-strand DNA breaks (Zhou et al., 2006). At 9 h, this was induced, at 24 h, it was maximal, and at 48 h, it was decreased to ∼25% (Supplementary DataSupplementary Data). Analyses of cell cycle distributions show no changes up to 12 h. At 18 h, BaP causes a drastic increase in S phase cells, with a concomitant decrease of cells in G1 and G2 phases (Fig. 1). This S phase accumulation lasts 6–18 h, whereafter the number of S phase cells decrease with a concomitant increase of G2/M-phase cells. Similar as for BaP-DNA adducts, the effects on cell cycle in the second experiment last shorter than in the first experiment. After about 54 h, the distribution of BaP-treated cells becomes similar to that of the controls. In the control cultures, the percentage of cells in S phase decreases after 18–24 h together with an increase of G1 phase cells. This may be because of depletion of the culture medium for essential compounds. Apoptotic cells levels did not increase because of BaP treatment, varied between 0.5 and 3% for control and BaP cultures, and were slightly higher in the first experiment (data not shown). The differences in adduct and cell cycle profiles between the two independent experiments are representative for the biological variability, which can exist in the HepG2 model. Large expression changes for CYP1A1 and NQO1 have been reported in HepG2 cells during 10 passages (Wilkening and Bader, 2003); this may explain the variability that we observed. Despite that, HepG2 are a very popular in vitro liver model for toxicologists and pharmacologists.
Levels of DNA adducts and cell cycle distribution in HepG2 cells following exposure to BaP. Cells are analyzed at various times points after starting the exposure; this was performed for two independent experiments (circles and triangles are for BaP-treated and vehicle control cells, respectively. closed and open symbols are for experiments 1 and 2, respectively). Adducts for vehicle control cells are not indicated, as they were almost zero.
Effects on Gene Expression
Differentially expressed genes.
Gene expression profiles were determined by DNA microarrays, which resulted in 48 hybridizations, four per time point. Differentially expressed genes were identified by reproducible similar changes over the four replicates, combined with ANOVA, resulting in 2808 differentially expressed genes (Supplementary Data “modulated genes.xls”). The same approach on randomized arrays resulted in not a single significantly affected gene, which indicates that few false positives were obtained. Visualization of the log base 2 expression ratios for these 2808 genes and all arrays in a heat map shows a high reproducibility in the expression profiles between the various replicate experiments and arrays (Supplementary Data). A few genes were verified by quantitative reverse transcriptase PCR (RT-PCR), namely CYP1A1, CYP1B1, CDKN1A, DDB2, ERCC1, and POLD4. With both methods, the patterns are similar, although the effects on gene expression as measured by RT-PCR are larger (Supplementary Data).
Pathway analyses by Reactome.
In order to examine the responses that may occur at the level of reactions and reaction pathways, the sets of differentially expressed genes per time point were analyzed in Reactome. An overview of the reaction pathways is in Figure 2 (Supplementary Data “Reactome and TF data”), thereby focusing on reaction pathways/processes that are considered most relevant for chemical carcinogenesis. Functional annotation based on overrepresentation analyses using GO terms and other gene sets was also performed using FatiGO+ (http://babelomics.bioinfo.cipf.es) (Al-Shahrour et al., 2005). Results added limited knowledge and thus are not presented.
Proportions of genes in reaction pathways from Reactome affected in HepG2 cells by BaP treatment, as observed in two independent experiments. At each time point/experiment, the proportion of differentially expressed genes in a Reactome pathway was calculated relative to the number of genes in that pathway as present in the set of 2808 modulated genes. Positive scores indicate the proportion of upregulated genes, negative scores that of downregulated genes.
Pathway analyses by T-profiler.
T-profiler was applied to identify affected functions using the complete data set, thus without preselecting modulated genes. Besides gene sets based on GO annotation and so on, also gene sets retrieved from literature were incorporated. Following an analysis per time point, in total, 155 gene sets were identified as being modulated (Supplementary Data “T-profiler data”). Hierarchical clustering was used to group these gene sets (Fig. 3). The main functional annotations of these clusters show upregulation of TP53 and NRF2 signaling pathways, apoptosis-related gene sets, and responses induced by genotoxic stressors. Downregulation consists of metabolisms of amines, sterols and lipids, (regulation of) transcription, and nucleosome assembly. This confirms the Reactome data; although T-profiler presents less in-depth information, it provides a better statistical evaluation.
Hierarchical clustering of gene sets that are significantly affected in HepG2 cells by treatment with BaP. The significantly affected gene sets were identified with the functional annotation tool T-profiler (Boorsma et al., 2005). The t values provided by T-profiler for two independent experiments were used for clustering with centroid linkage and Euclidean correlation. The tree was cut at 14 clusters, and the main annotations for these are presented. Positive and negative scores mean up- and downregulation of the gene set, respectively.
TF analyses.
The network-building algorithm on transcription regulation from MetaCore (GeneGo) was used to examine whether the modulated genes are connected to TF. Analyses of the 233 genes that showed at least a twofold change in three out of the four replicates revealed 30 significant networks centered around TF (Supplementary Data“MetaCore Networks.xls” and “Metacore Report 233 genes.doc”). The number of genes in each network varies from 10 to 48. An overview of the gene expression changes in the networks for these TF reaction pathways is shown in Figure 4 (Supplementary Data“MetaCore Networks.xls”). Again clear temporal differences are observed between the various TF networks, most showing an early or late upregulation, and only two showing a downregulation.
Proportions of genes in significantly affected transcription regulation networks centered around transcription factors as retrieved from MetaCore for gene expression changes induced in HepG2 cells by BaP treatment. At each time point/experiment, an inventory was made for the presence of differentially expressed in these transcription regulation networks (see Supplementary Data “MetaCore Networks.xls”), and the proportion was calculated relative to the number of genes in that network as present in the set of 233 modulated genes. Positive scores indicate the proportion of upregulated genes, negative scores of downregulated genes.
In order to examine whether these TF can be linked with the modulated genes in the Reactome pathways, an inventory was made of the TF-binding sites in promoter regions of these genes as well as of literature-based knowledge on regulation of the gene by the TF. The majority of the affected genes can be regulated by NF-κB (73% of the genes), followed in descending order by c-MYC, SRF, AP1, and E2F1 (from 42 to 34%) and TP53, CREB1, NRF2, and ESR1 (from 11 to 7%) (see Table 1 and Supplementary Data“Reactome and TF data”). All other TF hardly regulate any of the genes affected in the Reactome pathways. Many TF can also regulate a few to many of the other TF, sometimes also vice versa, thereby forming a complex regulatory network (see Supplementary Data).
Overview of the Percentage of Modulated Genes (MGs) in the Reactome Pathways That Can Be Regulated by a Transcription Factor (TF)
| Reactome pathways | Total No. of genes in pathway | No. of MGs in pathway | % of genes regulated by the TF | ||||||||
| NF-κB | c-MYC | SRF | AP1 | E2F1 | TP53 | CREB1 | NRF2 | ESR1 | |||
| Apoptosis | |||||||||||
| Regulation_of_Apoptosis | 44 | 9 | 67 | 44 | 33 | 67 | 11 | 0 | 11 | 0 | 0 |
| Extrinsic_Pathway_for_Apoptosis | 12 | 6 | 50 | 33 | 50 | 17 | 33 | 33 | 0 | 0 | 0 |
| Intrinsic_Pathway_for_Apoptosis | 27 | 8 | 75 | 25 | 25 | 13 | 13 | 50 | 25 | 0 | 13 |
| Apoptotic_execution_phase | 46 | 13 | 85 | 38 | 8 | 31 | 31 | 0 | 0 | 8 | 23 |
| Cell cycle checkpoints | |||||||||||
| G1/S_DNA_Damage_Checkpoints | 51 | 12 | 83 | 33 | 42 | 58 | 33 | 8 | 8 | 0 | 0 |
| G2/M_Checkpoints | 43 | 16 | 69 | 19 | 31 | 25 | 56 | 6 | 0 | 0 | 13 |
| Mitotic_Spindle_Checkpoint | 19 | 4 | 75 | 50 | 50 | 25 | 25 | 25 | 0 | 0 | 25 |
| Cell cycle, mitotic | |||||||||||
| G1_Phase | 14 | 3 | 67 | 33 | 67 | 0 | 33 | 33 | 0 | 0 | 0 |
| G1/S_Transition | 96 | 26 | 69 | 27 | 31 | 31 | 42 | 8 | 4 | 0 | 4 |
| S_Phase | 97 | 24 | 75 | 25 | 42 | 38 | 42 | 4 | 4 | 8 | 4 |
| G2/M_Transition | 14 | 5 | 80 | 60 | 60 | 40 | 20 | 40 | 0 | 0 | 0 |
| M_Phase | 82 | 19 | 79 | 37 | 53 | 11 | 37 | 11 | 16 | 5 | 11 |
| M/G1_Transition | 62 | 15 | 73 | 27 | 27 | 40 | 33 | 0 | 7 | 0 | 7 |
| APC/C-mediated_degradation_of_cell_cycle_proteins | 144 | 28 | 71 | 36 | 43 | 50 | 14 | 14 | 7 | 0 | 0 |
| DNA repair | |||||||||||
| Base_Excision_Repair | 18 | 4 | 75 | 25 | 50 | 75 | 25 | 0 | 0 | 25 | 0 |
| Global_Genomic_NER_(GG-NER) | 33 | 7 | 86 | 71 | 43 | 43 | 43 | 14 | 0 | 29 | 0 |
| Transcription-coupled_NER_(TC-NER) | 44 | 8 | 88 | 63 | 38 | 38 | 50 | 13 | 25 | 25 | 13 |
| Metabolism | |||||||||||
| Nucleotide_metabolism | 69 | 18 | 83 | 50 | 33 | 61 | 50 | 17 | 17 | 6 | 11 |
| Lipid_and_lipoprotein_metabolism | 131 | 39 | 69 | 31 | 41 | 44 | 38 | 0 | 21 | 8 | 13 |
| Metabolism_of_amino_acids | 158 | 39 | 69 | 41 | 38 | 46 | 41 | 0 | 15 | 3 | 10 |
| Metabolism_of_carbohydrates | 89 | 21 | 62 | 43 | 57 | 62 | 38 | 0 | 14 | 24 | 5 |
| Integration_of_energy_metabolism | 53 | 8 | 75 | 38 | 25 | 38 | 25 | 0 | 13 | 25 | 0 |
| Biological oxidations | |||||||||||
| Functionalization_of_compounds | 66 | 13 | 46 | 38 | 46 | 15 | 8 | 0 | 31 | 0 | 15 |
| Phase_II_conjugation | 47 | 17 | 88 | 65 | 76 | 29 | 29 | 12 | 0 | 24 | 0 |
| Gene expression | |||||||||||
| Transcription | 129 | 24 | 63 | 42 | 38 | 42 | 42 | 4 | 21 | 0 | 13 |
| Metabolism_of_noncoding_RNA | 23 | 3 | 100 | 67 | 33 | 0 | 67 | 0 | 0 | 0 | 0 |
| mRNA_Processing | 32 | 7 | 71 | 57 | 43 | 14 | 57 | 14 | 29 | 0 | 14 |
| Translation | 119 | 6 | 50 | 17 | 67 | 67 | 33 | 0 | 17 | 33 | 0 |
| Posttranslational_protein_modification | 40 | 6 | 83 | 83 | 0 | 67 | 17 | 17 | 17 | 17 | 17 |
| Average | 73 | 42 | 41 | 37 | 34 | 11 | 10 | 8 | 7 | ||
| Reactome pathways | Total No. of genes in pathway | No. of MGs in pathway | % of genes regulated by the TF | ||||||||
| NF-κB | c-MYC | SRF | AP1 | E2F1 | TP53 | CREB1 | NRF2 | ESR1 | |||
| Apoptosis | |||||||||||
| Regulation_of_Apoptosis | 44 | 9 | 67 | 44 | 33 | 67 | 11 | 0 | 11 | 0 | 0 |
| Extrinsic_Pathway_for_Apoptosis | 12 | 6 | 50 | 33 | 50 | 17 | 33 | 33 | 0 | 0 | 0 |
| Intrinsic_Pathway_for_Apoptosis | 27 | 8 | 75 | 25 | 25 | 13 | 13 | 50 | 25 | 0 | 13 |
| Apoptotic_execution_phase | 46 | 13 | 85 | 38 | 8 | 31 | 31 | 0 | 0 | 8 | 23 |
| Cell cycle checkpoints | |||||||||||
| G1/S_DNA_Damage_Checkpoints | 51 | 12 | 83 | 33 | 42 | 58 | 33 | 8 | 8 | 0 | 0 |
| G2/M_Checkpoints | 43 | 16 | 69 | 19 | 31 | 25 | 56 | 6 | 0 | 0 | 13 |
| Mitotic_Spindle_Checkpoint | 19 | 4 | 75 | 50 | 50 | 25 | 25 | 25 | 0 | 0 | 25 |
| Cell cycle, mitotic | |||||||||||
| G1_Phase | 14 | 3 | 67 | 33 | 67 | 0 | 33 | 33 | 0 | 0 | 0 |
| G1/S_Transition | 96 | 26 | 69 | 27 | 31 | 31 | 42 | 8 | 4 | 0 | 4 |
| S_Phase | 97 | 24 | 75 | 25 | 42 | 38 | 42 | 4 | 4 | 8 | 4 |
| G2/M_Transition | 14 | 5 | 80 | 60 | 60 | 40 | 20 | 40 | 0 | 0 | 0 |
| M_Phase | 82 | 19 | 79 | 37 | 53 | 11 | 37 | 11 | 16 | 5 | 11 |
| M/G1_Transition | 62 | 15 | 73 | 27 | 27 | 40 | 33 | 0 | 7 | 0 | 7 |
| APC/C-mediated_degradation_of_cell_cycle_proteins | 144 | 28 | 71 | 36 | 43 | 50 | 14 | 14 | 7 | 0 | 0 |
| DNA repair | |||||||||||
| Base_Excision_Repair | 18 | 4 | 75 | 25 | 50 | 75 | 25 | 0 | 0 | 25 | 0 |
| Global_Genomic_NER_(GG-NER) | 33 | 7 | 86 | 71 | 43 | 43 | 43 | 14 | 0 | 29 | 0 |
| Transcription-coupled_NER_(TC-NER) | 44 | 8 | 88 | 63 | 38 | 38 | 50 | 13 | 25 | 25 | 13 |
| Metabolism | |||||||||||
| Nucleotide_metabolism | 69 | 18 | 83 | 50 | 33 | 61 | 50 | 17 | 17 | 6 | 11 |
| Lipid_and_lipoprotein_metabolism | 131 | 39 | 69 | 31 | 41 | 44 | 38 | 0 | 21 | 8 | 13 |
| Metabolism_of_amino_acids | 158 | 39 | 69 | 41 | 38 | 46 | 41 | 0 | 15 | 3 | 10 |
| Metabolism_of_carbohydrates | 89 | 21 | 62 | 43 | 57 | 62 | 38 | 0 | 14 | 24 | 5 |
| Integration_of_energy_metabolism | 53 | 8 | 75 | 38 | 25 | 38 | 25 | 0 | 13 | 25 | 0 |
| Biological oxidations | |||||||||||
| Functionalization_of_compounds | 66 | 13 | 46 | 38 | 46 | 15 | 8 | 0 | 31 | 0 | 15 |
| Phase_II_conjugation | 47 | 17 | 88 | 65 | 76 | 29 | 29 | 12 | 0 | 24 | 0 |
| Gene expression | |||||||||||
| Transcription | 129 | 24 | 63 | 42 | 38 | 42 | 42 | 4 | 21 | 0 | 13 |
| Metabolism_of_noncoding_RNA | 23 | 3 | 100 | 67 | 33 | 0 | 67 | 0 | 0 | 0 | 0 |
| mRNA_Processing | 32 | 7 | 71 | 57 | 43 | 14 | 57 | 14 | 29 | 0 | 14 |
| Translation | 119 | 6 | 50 | 17 | 67 | 67 | 33 | 0 | 17 | 33 | 0 |
| Posttranslational_protein_modification | 40 | 6 | 83 | 83 | 0 | 67 | 17 | 17 | 17 | 17 | 17 |
| Average | 73 | 42 | 41 | 37 | 34 | 11 | 10 | 8 | 7 | ||
Overview of the Percentage of Modulated Genes (MGs) in the Reactome Pathways That Can Be Regulated by a Transcription Factor (TF)
| Reactome pathways | Total No. of genes in pathway | No. of MGs in pathway | % of genes regulated by the TF | ||||||||
| NF-κB | c-MYC | SRF | AP1 | E2F1 | TP53 | CREB1 | NRF2 | ESR1 | |||
| Apoptosis | |||||||||||
| Regulation_of_Apoptosis | 44 | 9 | 67 | 44 | 33 | 67 | 11 | 0 | 11 | 0 | 0 |
| Extrinsic_Pathway_for_Apoptosis | 12 | 6 | 50 | 33 | 50 | 17 | 33 | 33 | 0 | 0 | 0 |
| Intrinsic_Pathway_for_Apoptosis | 27 | 8 | 75 | 25 | 25 | 13 | 13 | 50 | 25 | 0 | 13 |
| Apoptotic_execution_phase | 46 | 13 | 85 | 38 | 8 | 31 | 31 | 0 | 0 | 8 | 23 |
| Cell cycle checkpoints | |||||||||||
| G1/S_DNA_Damage_Checkpoints | 51 | 12 | 83 | 33 | 42 | 58 | 33 | 8 | 8 | 0 | 0 |
| G2/M_Checkpoints | 43 | 16 | 69 | 19 | 31 | 25 | 56 | 6 | 0 | 0 | 13 |
| Mitotic_Spindle_Checkpoint | 19 | 4 | 75 | 50 | 50 | 25 | 25 | 25 | 0 | 0 | 25 |
| Cell cycle, mitotic | |||||||||||
| G1_Phase | 14 | 3 | 67 | 33 | 67 | 0 | 33 | 33 | 0 | 0 | 0 |
| G1/S_Transition | 96 | 26 | 69 | 27 | 31 | 31 | 42 | 8 | 4 | 0 | 4 |
| S_Phase | 97 | 24 | 75 | 25 | 42 | 38 | 42 | 4 | 4 | 8 | 4 |
| G2/M_Transition | 14 | 5 | 80 | 60 | 60 | 40 | 20 | 40 | 0 | 0 | 0 |
| M_Phase | 82 | 19 | 79 | 37 | 53 | 11 | 37 | 11 | 16 | 5 | 11 |
| M/G1_Transition | 62 | 15 | 73 | 27 | 27 | 40 | 33 | 0 | 7 | 0 | 7 |
| APC/C-mediated_degradation_of_cell_cycle_proteins | 144 | 28 | 71 | 36 | 43 | 50 | 14 | 14 | 7 | 0 | 0 |
| DNA repair | |||||||||||
| Base_Excision_Repair | 18 | 4 | 75 | 25 | 50 | 75 | 25 | 0 | 0 | 25 | 0 |
| Global_Genomic_NER_(GG-NER) | 33 | 7 | 86 | 71 | 43 | 43 | 43 | 14 | 0 | 29 | 0 |
| Transcription-coupled_NER_(TC-NER) | 44 | 8 | 88 | 63 | 38 | 38 | 50 | 13 | 25 | 25 | 13 |
| Metabolism | |||||||||||
| Nucleotide_metabolism | 69 | 18 | 83 | 50 | 33 | 61 | 50 | 17 | 17 | 6 | 11 |
| Lipid_and_lipoprotein_metabolism | 131 | 39 | 69 | 31 | 41 | 44 | 38 | 0 | 21 | 8 | 13 |
| Metabolism_of_amino_acids | 158 | 39 | 69 | 41 | 38 | 46 | 41 | 0 | 15 | 3 | 10 |
| Metabolism_of_carbohydrates | 89 | 21 | 62 | 43 | 57 | 62 | 38 | 0 | 14 | 24 | 5 |
| Integration_of_energy_metabolism | 53 | 8 | 75 | 38 | 25 | 38 | 25 | 0 | 13 | 25 | 0 |
| Biological oxidations | |||||||||||
| Functionalization_of_compounds | 66 | 13 | 46 | 38 | 46 | 15 | 8 | 0 | 31 | 0 | 15 |
| Phase_II_conjugation | 47 | 17 | 88 | 65 | 76 | 29 | 29 | 12 | 0 | 24 | 0 |
| Gene expression | |||||||||||
| Transcription | 129 | 24 | 63 | 42 | 38 | 42 | 42 | 4 | 21 | 0 | 13 |
| Metabolism_of_noncoding_RNA | 23 | 3 | 100 | 67 | 33 | 0 | 67 | 0 | 0 | 0 | 0 |
| mRNA_Processing | 32 | 7 | 71 | 57 | 43 | 14 | 57 | 14 | 29 | 0 | 14 |
| Translation | 119 | 6 | 50 | 17 | 67 | 67 | 33 | 0 | 17 | 33 | 0 |
| Posttranslational_protein_modification | 40 | 6 | 83 | 83 | 0 | 67 | 17 | 17 | 17 | 17 | 17 |
| Average | 73 | 42 | 41 | 37 | 34 | 11 | 10 | 8 | 7 | ||
| Reactome pathways | Total No. of genes in pathway | No. of MGs in pathway | % of genes regulated by the TF | ||||||||
| NF-κB | c-MYC | SRF | AP1 | E2F1 | TP53 | CREB1 | NRF2 | ESR1 | |||
| Apoptosis | |||||||||||
| Regulation_of_Apoptosis | 44 | 9 | 67 | 44 | 33 | 67 | 11 | 0 | 11 | 0 | 0 |
| Extrinsic_Pathway_for_Apoptosis | 12 | 6 | 50 | 33 | 50 | 17 | 33 | 33 | 0 | 0 | 0 |
| Intrinsic_Pathway_for_Apoptosis | 27 | 8 | 75 | 25 | 25 | 13 | 13 | 50 | 25 | 0 | 13 |
| Apoptotic_execution_phase | 46 | 13 | 85 | 38 | 8 | 31 | 31 | 0 | 0 | 8 | 23 |
| Cell cycle checkpoints | |||||||||||
| G1/S_DNA_Damage_Checkpoints | 51 | 12 | 83 | 33 | 42 | 58 | 33 | 8 | 8 | 0 | 0 |
| G2/M_Checkpoints | 43 | 16 | 69 | 19 | 31 | 25 | 56 | 6 | 0 | 0 | 13 |
| Mitotic_Spindle_Checkpoint | 19 | 4 | 75 | 50 | 50 | 25 | 25 | 25 | 0 | 0 | 25 |
| Cell cycle, mitotic | |||||||||||
| G1_Phase | 14 | 3 | 67 | 33 | 67 | 0 | 33 | 33 | 0 | 0 | 0 |
| G1/S_Transition | 96 | 26 | 69 | 27 | 31 | 31 | 42 | 8 | 4 | 0 | 4 |
| S_Phase | 97 | 24 | 75 | 25 | 42 | 38 | 42 | 4 | 4 | 8 | 4 |
| G2/M_Transition | 14 | 5 | 80 | 60 | 60 | 40 | 20 | 40 | 0 | 0 | 0 |
| M_Phase | 82 | 19 | 79 | 37 | 53 | 11 | 37 | 11 | 16 | 5 | 11 |
| M/G1_Transition | 62 | 15 | 73 | 27 | 27 | 40 | 33 | 0 | 7 | 0 | 7 |
| APC/C-mediated_degradation_of_cell_cycle_proteins | 144 | 28 | 71 | 36 | 43 | 50 | 14 | 14 | 7 | 0 | 0 |
| DNA repair | |||||||||||
| Base_Excision_Repair | 18 | 4 | 75 | 25 | 50 | 75 | 25 | 0 | 0 | 25 | 0 |
| Global_Genomic_NER_(GG-NER) | 33 | 7 | 86 | 71 | 43 | 43 | 43 | 14 | 0 | 29 | 0 |
| Transcription-coupled_NER_(TC-NER) | 44 | 8 | 88 | 63 | 38 | 38 | 50 | 13 | 25 | 25 | 13 |
| Metabolism | |||||||||||
| Nucleotide_metabolism | 69 | 18 | 83 | 50 | 33 | 61 | 50 | 17 | 17 | 6 | 11 |
| Lipid_and_lipoprotein_metabolism | 131 | 39 | 69 | 31 | 41 | 44 | 38 | 0 | 21 | 8 | 13 |
| Metabolism_of_amino_acids | 158 | 39 | 69 | 41 | 38 | 46 | 41 | 0 | 15 | 3 | 10 |
| Metabolism_of_carbohydrates | 89 | 21 | 62 | 43 | 57 | 62 | 38 | 0 | 14 | 24 | 5 |
| Integration_of_energy_metabolism | 53 | 8 | 75 | 38 | 25 | 38 | 25 | 0 | 13 | 25 | 0 |
| Biological oxidations | |||||||||||
| Functionalization_of_compounds | 66 | 13 | 46 | 38 | 46 | 15 | 8 | 0 | 31 | 0 | 15 |
| Phase_II_conjugation | 47 | 17 | 88 | 65 | 76 | 29 | 29 | 12 | 0 | 24 | 0 |
| Gene expression | |||||||||||
| Transcription | 129 | 24 | 63 | 42 | 38 | 42 | 42 | 4 | 21 | 0 | 13 |
| Metabolism_of_noncoding_RNA | 23 | 3 | 100 | 67 | 33 | 0 | 67 | 0 | 0 | 0 | 0 |
| mRNA_Processing | 32 | 7 | 71 | 57 | 43 | 14 | 57 | 14 | 29 | 0 | 14 |
| Translation | 119 | 6 | 50 | 17 | 67 | 67 | 33 | 0 | 17 | 33 | 0 |
| Posttranslational_protein_modification | 40 | 6 | 83 | 83 | 0 | 67 | 17 | 17 | 17 | 17 | 17 |
| Average | 73 | 42 | 41 | 37 | 34 | 11 | 10 | 8 | 7 | ||
Time-related effects with the functional data.
Numerous different approaches can be used to identify co-regulated genes with time-related patterns in time series experiments on microarray studies, including hierarchical clustering, or dedicated tools, such as CAGED (Ramoni et al., 2002), GQL (Costa et al., 2005), and STEM (Bar-Joseph et al., 2008; Ernst and Bar-Joseph, 2006). All these tools have been applied on our gene expression data set. Although STEM resulted in the most comprehensible information, this mainly confirmed the data from the functional annotation analyses mentioned above, which is thus not presented here.
STEM was also used to uncover time-related patterns in the reactions, pathways, biological processes, and TF networks as identified by Reactome, T-profiler, and MetaCore (Figs. 2–4), in combination with the phenotypical parameters (BaP levels in medium, BaP-DNA adducts, and various phases of the cell cycle). In these analyses, experiments 1 and 2 were used independently. Both positive and inverted correlations are observed, showing temporal changes in pathways and reactions in DNA repair, DNA damage response, cell cycle, and general metabolisms. These are summarized and visualized in a heat map in Figure 5 and further discussed below.
Pathways, gene sets, and TF networks showing a time-dependent correlation or inverse correlation with BaP in medium, BaP-DNA adduct levels, and cell cycle phases. The correlations were identified with the time series analysis tool STEM, using the data from Figures 1–4 and the same scoring scheme.
DISCUSSION
Chemical carcinogens may cause a multitude of effects in cells. These range from altering the cellular metabolic capacity to inducing DNA damage and repair, mutagenesis, and influencing cell growth and viability. As many of these processes are exerted through changing the expression of genes at the level of messenger RNA (mRNA), genome-wide transcriptome profiling over an extensive period and multiple time points will provide new insights in TF-networks related to cellular responses to chemical carcinogens.
In our study, transcriptome responses in HepG2 cells were time dependently compared with BaP exposure levels, DNA adduct levels, and distribution of cells over various phases of the cell cycle. The bioinformatic analyses focused on identification of the TF that may be involved in the sequential changes in gene expression profiles and can be related with the biological effects.
In order to identify the TF that are most relevant for the changes in mRNA profiles, the network-building algorithm on transcription regulation from MetaCore was used to examine whether the modulated genes are connected to TF. This revealed significant networks centered around 30 different TF. Sixteen of these TF can be linked with carcinogenesis based on queering Entrez Gene (www.ncbi.nlm.nih.gov/gene/), GeneCards (www.genecards.org) and UniProt (www.uniprot.org), many of them are proto-oncogenes and tumor suppressor genes.
To identify the TF that may be involved in the sequential changes in gene expression profiles and can be related with the biological effects, we have correlated temporal changes in phenotypical/biological parameters (BaP levels in medium, BaP-DNA adducts, and various phases of the cell cycle) with those in gene sets from Reactome reactions, biological pathways, and TF networks, as identified by the time series analyses tool STEM.
The responses in phenotypical/biological parameters have different temporal profiles (see Fig. 1 and Supplementary Data) and can be divided in:
(1) Immediate response to changes in BaP levels; immediately after adding BaP to the cells the levels drop quickly, which is most likely because of metabolism by CYP450 enzymes.
(2) Early response to changes in BaP-DNA adducts: between 3 and 6 h, the levels drastically increase, remain stable for several whereafter they decrease again.
(3) Medium response to effects on S and G1 phase of the cell cycle; between 15 and 18 h, these changes take place, remain stable for several whereafter they return to normal values again.
(4) Medium/late responses to effects on G2 phase of the cell cycle; after 15 h, the G2 levels decrease followed by a drastic increase at 18 to 30 h, and eventually a slow decrease until 60 h posttreatment.
The temporal correlations of these phenotypical/biological parameters with gene expression changes for reactions, pathways, and TF networks as revealed by the time series analyses tool STEM are presented in Figure 5. These will be further discussed.
Immediate Responses: Correlations with BaP Levels
Three pathways/processes show a positive or an inverted correlation with BaP levels. The correlation with TCDD responsive genes is as expected because also PAHs like BaP activate AhR (Ma, 2001). The correlation with TF networks around the proto-oncogenes AFX1 (involved in the regulation of the insulin signaling pathway) and E2A (involved in regulation of immunoglobulin gene expression) is unexpected. For both, no involvement in toxicological responses is known. An inverted correlation with BaP levels is observed with various diverse gene sets, namely “React_474-metabolism of carbohydrates,” “T-profiler-11-oxidoreductase activity,” and “T-profiler-5-various.”
Early Responses: Correlations with BaP-DNA Adducts
BaP-DNA adduct levels show a positive correlation with seven T-profiler and Reactome gene sets (DNA repair–related and DNA damage and stress response–related gene sets) and 10 TF networks and with one gene set an inverted correlation (T-profiler-8-chromatin).
The 10 TF networks whose genes correlate with BaP-DNA adducts include many TF, which are known to be involved in stress responses, cell proliferation, and apoptosis, such as AP1, NF-κB, c-MYC, p53, and E2F1. In agreement with this, many of these TF can regulate the Reactome gene sets that correlate with BaP-DNA adducts (Table 1). This suggests that these TF, or more likely, combinations thereof, are involved in regulating the early transcriptional responses for following exposure to BaP and the induction of DNA damage.
Very interesting are the effects on DNA repair pathways: both BER and transcription-coupled nucleotide excision repair (TC-NER) correlate with DNA adducts and are induced from 9/12 h until 48 h (Fig. 2). ROS are expected to be formed during the metabolic activation of BaP. ROS leads to DNA damage, such as 8-OHdG, which is repaired by BER (Izumi et al., 2003). Although at 30/36 h most BaP has been consumed, its metabolites and intermediates are still present and thus may cause an extended oxidative stress. Indeed, the induction of BER agrees very well with the period that the gene set “T-profiler_3-NRF2-regulated gene sets + motifs” and the NRF2 network are induced. The transcription factor Nrf2 is activated by ROS and induces the expression of many genes with antioxidant responsive elements in their promoter regions (Wang et al., 2007). The BER genes modulated by BaP are FEN1, MBD4, POLD4, and SMUG1, some of which are also affected by ultraviolet light (Garinis et al., 2005).
Repair of bulky DNA damages—like BaP-DNA adducts—is by NER (Braithwaite et al., 1998), and most of that is strand specific: the transcribed strand is faster repaired than the nontranscribed strand (Dreij et al., 2005; Mellon, 2005). This is well reflected by the gene expression alterations we observed: the TC-NER reaction pathway is upregulated for a shorter period than the GC-NER reaction pathway, namely from 9–42 to 9–60 h, respectively (Fig. 2). Six out of the nine NER genes modulated by BaP are in both GG-NER and TC-NER. DDB2 is specific for GG-NER and upregulated from 9 to 60 h; POLR2A and POLR2D are specific for TC-NER and upregulated at early time points or downregulated at late time points, respectively. Effects by BaP on DDB2 have been described before in some cell types but not in HepG2 cells (Hockley et al., 2006; Mahadevan et al., 2005; Schreck et al., 2009). Taken together, our data suggest that the DNA damage repair system is regulated to some extent at the level of transcription. This is novel as regulation of DNA damage repair pathways has been considered to occur mostly at posttranslational level through protein modifications rather than gene expression (Branzei and Foiani, 2008; Huang and D’Andrea, 2006).
The correlation of BaP-DNA adduct levels with many pathways representing the DNA damage response is as expected because of PAHs like BaP cause genotoxic damage and induce the DNA damage response pathway (Hockley et al., 2008; Jeffy et al., 2000). Although many pathways/networks related to apoptosis are induced, BaP causes only a marginal induction of apoptosis at the doses used here. Apparently, the apoptotic machinery is activated, but BaP-induced damages are repaired sufficiently to halt that process before execution. At higher BaP doses (30μM at 24 h) and by other PAHs, apoptosis is induced in HepG2 cells (Staal et al., 2007).
The set of genes identified by Hockley et al. (2006) to be affected by BaP in HepG2 cells is part of “T-profiler-1-ARE-Nrf2, many genotoxic stressors (e.g., UV, BaP), oxidative stress, apoptosis.” The largest overlap in affected genes for that study with ours is at 9–30 h.
Medium Responses: Correlations with Cells in S and G1 Phase of the Cell Cycle
Numbers of cells in G1 phase of the cell cycle are decreased from 18 to 24/36 h. None of the TF networks correlates over time with cells in S and G1 phase of the cell cycle. The TF that regulate the genes in the Reactome pathways “React_1590-G1 phase” and “React_2253-global genomic NER (GG-NER),” which both inversely correlate with cells in G1 phase, are mainly NF-κB, SRF, c-MYC, AP1, E2F1, and TP53 (Table 1). These are the same as for the reaction pathways correlating with DNA adduct levels and suggest that both the early and medium responses are mediated by the same complex network of TF.
G1-cell count shows a positive correlation with “T-profiler-8-chromatin” (which has an inverted relation with BaP-DNA adducts) and an inverse correlation with three others, namely “React_1590-G1 phase,” “React_2253-global genomic NER (GG-NER),” and “T-profiler-12-p53, genotoxic stressors.” All these also correlate with the number of cells in S phase, but positively. When the S phase cell counts are high and those for G1 phase cell are low, the DNA adduct levels are increased. This may explain their correlations with gene sets for DNA damage repair and DNA damage response.
The concomitant increase of cells in S phase with a decrease of cells in G1 phase (Fig. 1) suggests an arrest of cells in S phase and/or a stimulation of G1 cells into S phase. BaP is known to block DNA synthesis at replication forks, thereby causing an S phase arrest (Bi et al., 2005; Hockley et al., 2008; Jeffy et al., 2000; Staal et al., 2007).
Medium/Late Responses: Correlations with Cells in G2 Phase of the Cell Cycle
Three pathways/processes/reactions/networks have been shown to correlate with G2/M cells, which are decreased from 18 to 24/36 h and then increased from ∼24/36 to 54 h (Fig. 1), namely with “T-profiler-10-(regulation of) immune response, amino acid metabolism, steroid metabolism, lipid metabolism,” “T-profiler-7-developmental process,” and “TF-RARα.” Within the T-profiler-10 cluster, many gene sets relate to metabolism and these also show this downregulation followed by an upregulation (see Fig. 3). The activity of RARα (retinoic acid receptor α), a proto-oncogene in Promyelocytic Leukemia (Melnick and Licht, 1999), is inhibited by AhR (Widerak et al., 2006), which is activated by BaP. The link with metabolism related gene sets is not clear, however. This suggests that in the period when DNA adduct and numbers of cells S phase are high and when G1 and G2/M cell levels are low, the general metabolism of the cells is downregulated. When thereafter S phase levels decrease and G2/M increase, the demand of the cells for general metabolism is high. Once cell cycle distribution is back to baseline and DNA adducts have mostly been repaired, metabolic processes return to their normal levels.
G2/M cell levels show an inverted correlation with four pathways. Noteworthy, the inverted correlation is mainly with reactions involved in gene expression, transcription, and translation. Most of these show an initial upregulation followed by a downregulation thereafter. Also the order is remarkable: the induction period increases from 6 h for “React_11052-metabolism of non-coding RNA,” to ∼27 h for “React_1014-translation” and 39 h for “React_1788-transcription” (see Fig. 2). This suggests that the regulation of these reaction pathways is strictly organized. G2/M cell levels also correlate inversely with “React_1725-M/G1 transition,” which is not understood.
Reactome Pathways and TF
The inventory of the TF linked with the modulated genes in the Reactome pathways showed that 73% of the affected genes can be regulated by NF-κB, followed in descending order by c-MYC, SRF, AP1, and E2F1 (each of these can regulate at least 34% of the affected genes) (Table 1). This is observed for almost every Reactome pathway, irrespective of whether it shows an early, intermediate, or late transcriptional response. Based on the current data, it is not possible to indicate which TF or complex network of TF is responsible for the successive transcriptomic responses following BaP exposure. All these TF and many of the TF from the list of 30 can also regulate a few to many of the other TF, sometimes also vice versa, thereby forming a complex regulatory network (see Supplementary Data). Key regulators in the TF network appear to be p53, NF-κB, c-FOS, c-JUN, and c-MYC. Indeed, cross talk, interactions, or collaboration are not uncommon among these TF (Kensler et al., 2010; Schneider et al., 2010; Shyu et al., 2008; You and Mak, 2005).
These all indicates that a complex network of a few TF is responsible for the majority of the transcriptional changes induced by BaP. This network hardly changes over time, despite that the transcriptional profiles clearly alter, suggesting that also other regulatory mechanisms are involved. Fine tuning of the transcriptional responses in a time-dependent manner may be by other TF, such as p53, NRF2, and CREB1. Other processes may also be involved in that, such as posttranscriptional regulation at the level of miRNA or posttranslational regulation at the level of (de)phosphorylation, ubiquitinylation, and many others.
NF-κB is a transcription regulator that is activated by various intra- and extracellular stimuli, such as cytokines, oxidant-free radicals, ultraviolet irradiation, and bacterial or viral products, and also by BaP (Bolotina et al., 2007). Activated NF-κB translocates into the nucleus and stimulates the expression of genes involved in a wide variety of biological functions. NF-κB is a pleiotropic TF, which is present in almost all cell types and is involved in many biological processed, such as inflammation, immunity, differentiation, cell growth, tumorigenesis, and apoptosis (www.genecards.org).
c-Myc is a TF that is believed to regulate expression of 15% of all genes through binding on Enhancer Box sequences and recruiting histone acetyltransferases (Cotterman et al., 2008). The protein encoded by this gene is a multifunctional nuclear phosphoprotein that plays a role in cell cycle progression, apoptosis, and cellular transformation and can be affected by BaP (Fields et al., 2004).
Serum response factor (SRF) is a TF that binds to the serum response element in the promoter region of target genes. This protein regulates the activity of many immediate-early genes, e.g., c-FOS, and thereby participates in cell cycle regulation, apoptosis, cell growth, and cell differentiation. In literature, no data could be found indicating that BaP or other PAHs activate SRF (www.genecards.org).
The activator protein 1, AP-1, is a TF that is a heterodimeric protein composed of proteins belonging to the c-FOS, c-JUN, ATF, and JDP families. It regulates gene expression in response to a variety of stimuli, including cytokines, growth factors, stress, and bacterial and viral infections, and also by BaP (Bolotina et al., 2007). AP-1 in turn controls a number of cellular processes, including differentiation, proliferation, and apoptosis.
E2F1 is a member of the E2F family of TF. The E2F family plays a crucial role in the control of cell cycle. It can mediate both cell proliferation and p53-dependent/independent apoptosis. E2F-1 binds preferentially RB1 protein, in a cell cycle–dependent manner. In literature, no data could be found indicating that BaP or other PAHs activate E2F1.
Final Remarks
We believe that the current detailed time series study provides an unprecedented insight in time-dependent interactions of the affected processes and pathways caused by the full carcinogen BaP. All these new observations indicate that the cellular transcriptome is strictly regulated, with main responses once DNA damages are generated, and adapts to the changing needs of cells following acute exposure to a carcinogen. Effects of BaP on gene expression changes as measured by DNA microarrays in HepG2 cells have been published before (Hockley et al., 2006, 2007; Staal et al., 2006, 2007). These papers, however, did not identify the TF that might be involved in regulation of the affected genes. Thus, many of our findings were not previously reported. We were able to conduct more in-depth analyses because we investigated many time points over a long period, compared transcription profiles with several phenotypic data, focused on TF and TF networks, and applied software dedicated for time series microarray studies.
CONCLUSIONS
Based on all the information presented in this paper, we conclude/hypothesize that a complex network of a few TF is responsible for the majority of the transcriptional changes induced by BaP. This network consisting of NF-κB, c-MYC, SRF, AP1, and E2F1 hardly changes over time, despite that the transcriptional profiles clearly alter. This suggests that also other regulatory mechanisms are involved.
FUNDING
Netherlands Genomics Initiative.
Part of the analyses was performed using BRB ArrayTools developed by Dr Richard Simon and Amy Peng Lam.
Comments