In vivo – in vitro toxicogenomic comparison of TCDD-elicited gene expression in Hepa1c1c7 mouse hepatoma cells and C57BL/6 hepatic tissue

Abstract

Background

In vitro systems have inherent limitations in their ability to model whole organism gene responses, which must be identified and appropriately considered when developing predictive biomarkers of in vivo toxicity. Systematic comparison of in vitro and in vivo temporal gene expression profiles were conducted to assess the ability of Hepa1c1c7 mouse hepatoma cells to model hepatic responses in C57BL/6 mice following treatment with 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD).

Results

Gene expression analysis and functional gene annotation indicate that Hepa1c1c7 cells appropriately modeled the induction of xenobiotic metabolism genes in vivo. However, responses associated with cell cycle progression and proliferation were unique to Hepa1c1c7 cells, consistent with the cell cycle arrest effects of TCDD on rapidly dividing cells. In contrast, lipid metabolism and immune responses, representative of whole organism effects in vivo, were not replicated in Hepa1c1c7 cells.

Conclusion

These results identified inherent differences in TCDD-mediated gene expression responses between these models and highlighted the limitations of in vitro systems in modeling whole organism responses, and additionally identified potential predictive biomarkers of toxicity.

Background

Advances in microarray and related technologies continue to revolutionize biomedical research and are being incorporated into toxicology and risk assessment. These technologies not only facilitate a more comprehensive elucidation of the mechanisms of toxicity, but also support mechanistically-based quantitative risk assessment [1–5]. In addition, these technologies are being used to develop predictive toxicity screening assays to screen drug candidates with adverse characteristics earlier in the development pipeline in order to prioritize resources and maximize successes in clinical trials [6–8]. Comparable screening strategies are also being proposed to rank and prioritize commercial chemicals, natural products, and environmental contaminants that warrant further toxicological investigation. Traditionally, rodent models or surrogates for ecologically-relevant species are typically used in regulatory testing. However, public and regulatory pressure, especially in Europe, seek to minimize the use of animals in testing [9]. Similar policies in the US, such as the ICCVAM Authorization Act of 2000, provide guidelines to facilitate the regulatory acceptance of alternative testing methods. These initiatives combined with the need to assess an expanding list of drug candidates and commercial chemicals for toxicity, have increased demand for the development and implementation of high-throughput in vitro screening assays that are predictive of toxicity in humans and ecologically-relevant species.

Various in vitro hepatic models including the isolated perfused liver, precision cut liver slices, isolated primary liver cells and a number of immortalized liver cell lines, have been used as animal alternatives [10]. In addition to providing a renewable model, in vitro systems are a cost-effective alternative and are amenable to high-throughput screening. These models, particularly immortalized cell lines, also allow for more in-depth biochemical and molecular investigations, such as over-expression, knock-down, activation or inhibition strategies, thus further elucidating mechanisms of action. However, inherent limitations in the ability of cell cultures to model whole organism responses must also be considered when identifying putative biomarkers for high-throughput toxicity screening assays, and elucidating relevant mechanisms of toxicity that support quantitative risk assessment. Despite several in vitro toxicogenomic reports [11–13], few have systematically examined the ability of in vitro systems to predict in vivo gene expression profiles in response to chemical treatment [10, 14].

2,3,7,8-Tetrachlorodibenzo-p-dioxin (TCDD) is a widespread environmental contaminant that elicits a number of adverse effects including tumor promotion, teratogenesis, hepatotoxicity, and immunotoxicity as well as the induction of several metabolizing enzymes [15]. Many, if not all of these effects, are due to alterations in gene expression mediated by the aryl hydrocarbon receptor (AhR), a basic-helix-loop-helix-PAS (bHLH-PAS) transcription factor [15, 16]. Ligand binding to the cytoplasmic AhR complex triggers the dissociation of interacting proteins and results in the translocation of the ligand-bound AhR to the nucleus where it heterodimerizes with the aryl hydrocarbon receptor nuclear translocator (ARNT), another member of the bHLH-PAS family. The heterodimer then binds specific DNA elements, termed dioxin response elements (DREs), within the regulatory regions of target genes leading to changes in expression that ultimately result in the observed responses [17]. Although the role of AhR is well established, the gene regulatory pathways responsible for toxicity are poorly understood and warrant further investigation to assess the potential risks to humans and ecologically relevant species.

Hepa1c1c7 cells and C57BL/6 mice are well-established models routinely used to examine the mechanisms of action of TCDD and related compounds. In this study, TCDD-elicited temporal gene expression effects were systematically compared in order to assess the ability of Hepa1c1c7 cells to replicate C57BL/6 hepatic tissue responses. Our results indicate that several phase I and II metabolizing enzyme responses are aptly reproduced. However, many responses were model-specific and reflect inherent in vitro and in vivo differences that must be considered in mechanistic studies and during the selection of biomarkers for developing toxicity screening assays.

Results

In vitro microarray data analysis

Temporal gene expression profiles were assessed in Hepa1c1c7 wild type cells following treatment with 10 nM TCDD using cDNA microarrays with 13,362 spotted features. Empirical Bayes analysis of the in vitro time course data identified 331 features representing 285 unique genes with a P1(t) value greater than 0.9999 at one or more time points, and differential expression greater than ± 1.5 fold relative to time-matched vehicle controls. The number of differentially regulated genes gradually increased from 1 to 24 hrs, followed by a slight decrease at 48 hrs (Figure 1A). In vitro dose-response data performed at 12 hrs with TCDD covering 6 different concentrations (0.001, 0.01, 0.1, 1.0, 10 and 100 nM), identified 181 features representing 155 unique genes (P1(t) > 0.9999 and an absolute fold change > 1.5 at one or more doses; Figure 1B). Complete in vitro time course and dose-response data are available in Additional file 1 and 2, respectively.

Figure 1

Number of genes differentially regulated (P 1(t) > 0.9999 and Ifold changel > 1.5-fold) as measured by microarray analysis for the (A) time course and (B) and dose-response studies in mouse hepatoma Hepa1c1c7 cells. For the time course study, cells were treated with 10 nM TCDD and harvested at 1, 2, 4, 8, 12, 24 or 48 hrs after treatment. Cells for the 12 hr dose-response study were treated with 0.001, 0.01, 0.1, 1.0, 10 and 100 nM of TCDD

As a control, the gene expression effects elicited by 10 nM TCDD in ARNT-deficient c4 Hepa1c1c7 mutants [18] were examined at 1 and 24 hrs (data not shown). Only ATPase, H⁺ transporting, V1 subunit E-like 2 isoform 2 (Atp6v1e2) and SUMO/sentrin specific peptidase 6 (Senp6) exhibited a significant change in expression using the same criteria (P1(t) > 0.9999 and an absolute fold change > 1.5). Neither Atp6v1e2 nor Senp6 were among the active genes in wild-type Hepa1c1c7 cells or in C57BL/6 liver samples [19]. These results provide further evidence that the AhR/ARNT signaling pathway mediates TCDD-elicited gene expression responses, which are consistent with in vivo microarray results with AhR knockout mice [20].

Hierarchical clustering of the genes expressed in Hepa1c1c7 time course indicate that 2 and 4 hrs were most similar, as were 8 and 12 hrs, and 24 and 48 hrs, while the 1 hr time point was segregated (Figure 2A). A strong dose-response relationship was also evident with clusters sequentially branching out with increasing concentration (Figure 2B). At 12 hrs, 117 genes were differentially expressed with 112 exhibiting a dose-dependent response. Moreover, the fold changes measured in both the time course and dose-response studies using 10 nM TCDD were comparable. For example, xanthine dehydrogenase (Xdh) and NAD(P)H dehydrogenase, quinone 1 (Nqo1) were induced 2.39- and 4.89-fold respectively in the time course and 2.93- and 4.71-fold in the dose-response study. There is a strong correlation (R = 0.97) between the differentially expressed genes at 12 hrs in the time course with the differentially regulated genes in the dose-response study at 10 nM, demonstrating the reproducibility between independent studies and providing further evidence that these genes are regulated by TCDD.

Figure 2

Hierarchical clustering of the differentially regulated gene lists for A) temporal and B) dose-response microarray studies in mouse hepatoma Hepa1c1c7 cells. The results illustrate time and dose-dependent clustering patterns. From the A) temporal results, the early (2 hr and 4 hr), intermediate (8 hr and 12 hr) and late (24 hr and 48 hr) time points cluster separately while the 1 hr time point clusters alone. Results from the B) dose-response show that highest doses clustered together, while the remaining doses branched out in a dose-dependent manner

The list of temporally regulated genes was subjected to k-means clustering using the standard correlation distance metrics. Five k-means clusters best characterized the dataset and identified clusters representing A) up-regulated early and sustained, B) up-regulated intermediate and sustained, C) up-regulated intermediate, D) up-regulated immediate and E) down-regulated late (Figure 3). These were comparable to the k-means clusters identified in hepatic tissue of C57BL/6 mice following treatment with 30 μg/kg TCDD [19]. Although, no discernable functional category is over-represented in any one cluster, the sustained up-regulation of early (Cluster A) and intermediate (Cluster B) responding genes include classic TCDD-responsive genes such as cytochrome P450, family 1, subfamily a, polypeptide 1 (Cyp1a1), Xdh and Nqo1. Many down-regulated late genes were associated with cell cycle regulation such as myelocytomatosis oncogene (Myc). Additionally, targets of Myc, including cyclin D1 and ornithine decarboxylase (Odc1), were also down-regulated suggesting a mechanism for cell cycle arrest [21–23], a common in vitro response to TCDD.

Figure 3

K-means clustering of temporally differentially regulated genes in vitro. Five k-mean clusters corresponding to (A) up-early and sustained, (B) up-intermediate and sustained, (C) up-regulated intermediate, (D) up-regulated immediate and (E) down-regulated late. Time and expression ratio are indicated on the x- and y-axis respectively. The color of individual gene expression profiles reflects the expression ratio observed at 24 hrs

Classification of gene expression responses for common regulated genes

Using the same filtering criteria (P1(t) > 0.9999 and an absolute fold change > 1.5), 678 features representing 619 unique genes were differentially expressed as previously reported in a time course study conducted in hepatic tissue from C57BL/6 mice orally gavaged with 30 μg/kg TCDD [19]. The number of responsive in vivo genes and their temporal expression patterns closely paralleled the results from this in vitro study. The fewest number of active genes was observed at 2 hrs, followed by a large increase at 4 hrs, which was sustained to 72 hrs. However, the substantial increase in expressed in vivo genes at 168 hrs was attributed to triglyceride accumulation and immune cell infiltration, which was not observed in Hepa1c1c7 cells. This list of 619 of in vivo genes served as the basis for subsequent comparisons against TCDD-elicited in vitro responses.

Comparison of in vitro and in vivo differentially expressed gene lists identified common and model specific responses (Figure 4A). TCDD treatment resulted in a total of 838 regulated genes in either model and with 67 common to both. TCDD elicited 218 gene expression changes unique to Hepa1c1c7 cells while 552 genes were specific to C57BL/6 hepatic samples. Although 67 genes were regulated in both models, not all possessed similar temporal patterns of expression. Contingency analysis using a 2 × 2 table and the χ² test resulted in a p-value < 0.001 (α = 0.05) illustrate a statistically significant association between the lists of differentially regulated genes in vitro and in vivo. Further stratification revealed genes that were either induced in both models (class I), repressed in both models (class II), induced in vivo while repressed in vitro (class III), or repressed in vivo while induced in vitro (class IV; Figure 4B). Genes regulated in a similar fashion in both models (classes I and II) accounted for 49 of the 67 common active genes, while the remaining genes exhibited divergent expression profiles (classes III and IV). Hierarchical clustering of the temporal expression values for the 67 overlapping genes identified the same four classes (Figure 4C). The pattern across model and time illustrates that the earliest time points (i.e. 1 hr in vitro and 2 hr in vivo time points) cluster together while the remaining clusters branch into in vitro or in vivo clusters according to time. These results suggest that potential biomarkers of acute TCDD-mediated responses may best be predicted by the immediate-early in vitro gene responses.

Figure 4

Comparison of common significant in vitro and in vivo TCDD elicited time-dependent gene expression changes. A) 285 differentially regulated in vitro genes and 619 differentially regulated in vivo genes were identified, with 67 genes common to both studies. B) The temporal gene expression profiles from both studies were categorized into (I) induced in both, (II) repressed in both, (III) induced in vivo and repressed in vitro, and (IV) repressed in vitro and induced in vivo. C) Hierarchical clustering identified similar classification groups. Clustering across both time and model, separated samples from in vitro and in vivo, with the exception of the early time points from both studies (1 hr in vitro and 2 hr in vivo), which clustered together. * identifies in vitro time points

In vitro and in vivo induced genes (class I) include xenobiotic and oxidoreductase enzymes such as abhydrolase domain containing 6 (Abhd6), Cyp1a1, dehydrogenase/reductase (SDR family) member 3 (Dhrs3), Nqo1, prostaglandin-endoperoxide synthase 1 (Ptgs1), UDP-glucose dehydrogenase (Ugdh) and Xdh (Table 1). These genes have previously been reported to be TCDD-responsive [19, 24], with Cyp1a1 and Nqo1 being members of the "AhR gene battery" [25]. Glutathione S-transferase, alpha 4 (Gsta4) was also induced in vitro and in vivo, 1.7- and 2.0-fold respectively, consistent with TCDD-mediated induction of phase I and II metabolizing enzymes. Of the 35 genes responding similarly in both models, approximately 71% of were similarly up-regulated (class I) while the remaining genes were repressed across both models (class II). Repressed class II genes include minichromosome maintenance deficient 6 (Mcm6), glycerol kinase (Gyk) and ficolin A (Fcna) (repressed 1.6-, 1.6- and 1.7-fold in vitro, respectively). Overall, repressed genes did not share any common discernable biological function.

Table 1 Classification of common differentially regulated temporal gene expression responses to TCDD in both in vitro and in vivo models

Forty-two of the 67 common differentially expressed genes were dose responsive at 12 and 24 hrs in vitro and in vivo, respectively, further suggesting the role of the AhR in mediating these responses. Microarray-based EC₅₀ values spanned at least 3 orders of magnitude ranging from 0.05 μg/kg to >150 μg/kg in vivo, and 0.00118 nM to 2.4 nM in vitro (Table 1). Cyp1a1, the prototypical marker of TCDD exposure, had EC₅₀ values of 0.05 μg/kg and 0.014 nM, in vivo and in vitro respectively, and was induced 38-fold in both time course studies. Complete data sets for the in vivo time course and dose-responses experiments are available in Additional file 3 and 4.

Of the 67 overlapping genes, 18 exhibited divergent temporal profiles (classes III and IV). Class III contains 12 genes induced in vivo but repressed in vitro, while 6 were repressed in vivo and induced in vitro (class IV). Example genes include Myc (class III) and B-cell translocation gene 2 (Btg2, class IV) which are both involved in regulating cell cycle progression [23, 26–30]. Myc was induced 3.7-fold in vivo and repressed 2.2-fold in vitro, while Btg2 was repressed 1.8-fold in vivo and induced 1.5-fold in vitro.

In addition to the regulated genes common to both models, 218 in vitro- and 559 in vivo-specific genes were identified. Many of the unique in vitro responses are involved in cell cycle regulation, including cyclins D1 and B2 (Table 2). Cyclin D1, which complexes with cyclin-dependent kinase 4 (Cdk4) to regulate the progression from G₁ to S phase [31, 32], was down-regulated early and repressed 1.7-fold to 48 hrs. Furthermore, cyclin B2 and cell division cycle 2 homolog A (Cdc2a) which interact to form an active kinase required for G₂ promotion, were down-regulated, 1.8-fold and 1.5-fold, respectively. In addition to cell cycle related genes, UDP glucuronosyltransferase 1 family, polypeptide A2 (Ugt1a2), a phase II metabolizing enzyme, was induced 2.8-fold in vitro, but not significantly regulated in vivo.

Table 2 Examples of TCDD-elicited gene expression responses unique to Hepa1c1c7 cells

Analysis of the C57BL/6 hepatic time course identified 552 unique genes that were solely regulated in vivo. This included TCDD induced transcripts for microsomal epoxide hydrolase 1 (Ephx1) and carbonyl reductase 3 (Cbr3) which both function as xenobiotic metabolizing enzymes. Notch gene homolog 1 (Notch1) and growth arrest specific 1 (Gas1) which are both associated with development and differentiation but serve undetermined roles in the liver, were also induced by TCDD (Table 3). Genes related to immune cell accumulation were also specific to the in vivo study, coincident with immune cell accumulation at 168 hr as determined by histopathological examination [19].

Table 3 Examples of TCDD-elicited gene expression responses unique to C57BL/6 hepatic tissue

Comparison of basal gene expression levels in Hepa1c1c7 cells and hepatic tissue

In order to further investigate differences in gene expression levels, Hepa1c1c7 cells and C57BL/6 liver samples were directly compared by competitive hybridization on the same array, to identify basal gene expression level differences. Subsequent linear regression analysis of the mean normalized signal intensities from the untreated samples resulted in a correlation value of R = 0.75 (Figure 5), which is consistent with basal gene expression comparisons of various in vitro rat hepatic systems against whole livers, where correlation values decreased between liver slices (R = 0.97), primary cells (R = 0.85), BRL3A (R = 0.3) and NRL clone 9 (R = 0.32) rat liver cell lines [10]. Overall, the correlation illustrates reasonable concordance in basal gene expression levels between the two models. However, data points which deviate from the fitted line indicate differences in the basal expression of individual genes between the Hepa1c1c7 cells and hepatic tissue from C57BL/6 mice. Although there are differences, they may be negligible if the TCDD-elicited responses are conserved in vitro and in vivo. Complete microarray data for the untreated comparisons are available in Additional file 5.

Figure 5

Comparison of Hepalclc7 cell and C57BL/6 hepatic tissue basal gene expression. Untreated samples from Hepalclc7 cells and hepatic tissue from immature ovariectomized C57BL/6 mice taken at 0 hrs were competitively hybridized to the 13,362 feature cDNA microarray. Log2 normalized signal intensities were plotted for in vitro versus in vivo data to generate the correlation coefficient. The linear correlation coefficient R was 0.75 between in vitro and in vivo models

The relative basal expression of the 67 common active features was further investigated (Figure 5). In general, class I (i.e. induced in both models) genes fell close to the regression line, indicating that the basal expression of induced genes were comparable as were their in vitro and in vivo responses to TCDD. In contrast, basal expression levels of class III genes (i.e. induced in vivo while repressed in vitro) were generally higher in the Hepa1c1c7 cells, while levels in class II and IV (i.e. repressed in both models and repressed in vivo while induced in vitro, respectively) genes were scattered around the fitted linear line in Figure 5.

Quantitative real-time PCR verification of microarray responses

In total, 14 in vitro and 24 in vivo responsive genes representing common and model-specific genes were verified by quantitative real-time PCR (QRTPCR) (see Additional file 6). Of the selected genes regulated in both models, all displayed temporal patterns comparable to the microarray data (Figure 6). For example, Xdh, Myc and fatty acid binding protein (Fabp5) exhibited good agreement in fold change and temporal expression pattern when comparing microarray and QRTPCR data. However, significant data compression was evident when comparing in vitro and in vivo Cyp1a1 induction by QRTPCR, although in vitro and in vivo microarray induction levels were comparable. Previous studies suggest this is likely due to the limited fluorescence intensity range (0 – 65,535) of microarrays resulting in signal saturation and compression of the true magnitude of induction of transcript levels [33, 34]. Cross hybridization of homologous probes to a given target sequence on the microarray may also be a contributing factor, especially in comparison to other, more gene-specific measurement techniques [35].

Figure 6

Quantitative real-time PCR verification of in vitro and in vivo microarray results. The same RNA used for cDNA microarray analysis was examined by QRTPCR. All fold changes were calculated relative to time-matched vehicle controls. Bars (left axis) and line (right axis) represent data obtained by QRTPCR and cDNA microarrays, respectively. Genes are indicated by official gene symbols, and results are the average of four biological replicates. Classes refer to the respective classification categories as illustrated in Figure 4B. Error bars represent the standard error of measurement for the average fold change. *p < 0.05 for QRTPCR

Discussion

Microarrays have become an invaluable tool in toxicogenomics for comprehensively characterizing gene expression responses following treatment with an environmental contaminant, commercial chemical, natural product or drug as well as for investigating complex mixtures relevant to human and wildlife exposures. An emerging consensus suggests that toxicogenomics will accelerate drug development and significantly improve quantitative risk assessments [36, 37]. In addition, toxicogenomics supports the development and refinement of predictive in vitro high-throughput toxicity screening assays that can be used as alternatives to traditional in vivo testing. Ideally, in vitro high-throughput toxicity screens can be used to rank and prioritize drug candidates, environmental contaminants, and commercial chemicals, which warrant further development or testing. Although in vitro responses are assumed to reflect a subset of comparable in vivo responses, few studies have completed a comprehensive and systematic comparison. This study closely examined two well-established models, and comprehensively compared the TCDD-elicited gene expression to assess the predictive value of in vitro systems.

Comparative analysis of Hepa1c1c7 cell and hepatic C57BL/6 microarray data identified 67 differentially expressed genes co-regulated by TCDD. Four classes based on their temporal expression patterns were identified (Figure 4B and 4C), with 42 of the 67 common regulated genes exhibiting dose-response characteristics in both models. In vitro EC₅₀ values ranged from 0.001182 nM to 2.4 nM, while in vivo the values ranged from 0.05 μg/kg to >150 μg/kg. The wide range of EC₅₀ values illustrate the varying sensitivity of regulated genes to TCDD in both models.

Hepa1c1c7 cells and hepatic tissue from C57BL/6 mice are the prototypical models used to investigate the mechanisms of action of TCDD and other related compounds and both exhibited the classic induction of phase I and II metabolizing enzymes including Cyp1a1 and Nqo1 [38, 39]. Gsta4 and Xdh were also up-regulated in both models further demonstrating Hepa1c1c7 cells as a suitable model for investigating TCDD-regulated induction of xenobiotic metabolizing genes. In addition to these genes, the responses of Nqo1, Ugdh and Tnfaip2 were also conserved across models and were categorized as class I genes (similarly induced in both models; Figure 4B and 4C). However, Gsta2 was induced in vivo while no significant effect was detected in Hepa1c1c7 cells, and Ugt1a2 was induced in vitro but not differentially expressed in C57BL/6 hepatic tissue. Although many phase I and II metabolizing enzyme responses were conserved, differences exist that may limit Hepa1c1c7 cells from accurately modeling the full spectrum of in vivo hepatic responses elicited by TCDD.

A direct comparison of untreated Hepa1c1c7 cells and C57BL/6 hepatic tissue was performed to further investigate innate differences between the two models. Comparison of the normalized signal intensities revealed a good correlation (R = 0.75) between in vitro and in vivo basal expression levels (Figure 5). This illustrates that many genes are basally expressed to similar levels in both models as illustrated by the cluster of class I (similarly induced genes) closely surrounding the fitted line. Although a correlation exists, there are still differences in basal expression which may be associated with the origins of the models (i.e. normal hepatic tissue versus hepatoma derived Hepa1c1c7 cells), as well as the inability of in vitro systems to effectively model complex interactions between different cell types (e.g. Kupffer and stellate cells). For example, Myc, a G₁ to S phase cell cycle regulator [23, 26–29], was repressed in vitro while being induced in vivo and the model-specific responses may be related to difference in basal expression levels between the two models (Table 1). The levels of Myc transcripts in untreated Hepa1c1c7 cells were higher relative to untreated C57BL/6 hepatic tissue, consistent with the proliferative state of the in vitro system (data not shown). Examination of other class III genes suggests that they are more highly expressed in vitro when compared to in vivo (Figure 5). Consequently, differences in basal expression may be a factor contributing to divergent in vitro – in vivo responses. Another possible source for the model-specific responses may be related to DNA methylation status of the promoter region of TCDD-responsive genes in either model. DNA methylation results in gene silencing [40, 41] and a previous study with Hepa1c1c7 has shown that TCDD-elicited gene expression responses are influenced by DNA methylation status [42]. The differing methylation states between the in vitro and in vivo systems may further contribute to the model-specific gene expression responses.

Many in vitro specific gene expression responses elicited by TCDD were associated with cell cycle progression and cell cycle arrest. Myc and its downstream target, cyclin D1, which forms a kinase complex with Cdk4 [43, 44] were both repressed by TCDD. In contrast, Cdkn1a, an inhibitor of cyclin-dependent kinase 2 (Cdk2)-cyclin E complex kinase activity [43], was induced. Inactivation of the Cdk2-cyclin E complex prevents the phosphorylation of pRb resulting in cell cycle arrest during G₁. Additionally, the in vitro induction of Btg2 suggests an alternative mechanism for cell cycle arrest during the G₂ phase. Constitutively active BTG2 in human leukemia U937 cells, induces G₂/M cell cycle arrest by inhibiting the formation of the cyclin B1 and Cdc2 complex, thereby inhibiting the active kinase function of the complex [30]. Collectively, these results corroborate and extend previous in vitro TCDD-mediated cell cycle arrest studies [45–48].

TCDD treatment resulted in a number of divergent gene responses across both models as represented by classes III and IV (Figures 4B and 4C). Genes related to immune cell accumulation, including major histocompatibility complex (MHC) molecules were only observed in vivo, and are likely a response to hepatic damage mediated by ROS or fatty accumulation and therefore independent of direct AhR action [19]. This is characteristic of the complex interaction between different cell types responding to liver injury that cannot be modeled in homogenous cultures of cells.

Pharmacokinetics may also contribute to response differences between the two models. Hepa1c1c7 cells were directly treated, whereas in vivo, TCDD must first be delivered to the liver and targeted cells prior to eliciting its effects. Additionally, C57BL/6 studies were able to be carried out to 168 hrs following TCDD treatment, while in vitro studies were limited to 48 hrs to minimize potentially confounding effects due to cell confluency. However, early responses associated with classes I and II (induced or repressed in both models; Figure 4B and 4C) are well conserved and exhibit comparable levels of induction or repression in both models. Hierarchical clustering of the common active genes (Figure 4C) illustrates gene induction occurs early while gene repression occurs later in both models. Clustering across both time and model revealed that gene expression profiles at 1 hr in vitro and 2 hr in vivo were most similar. This clustering pattern implies that early in vitro responses may accurately model early in vivo gene expression effects.

Conclusion

Comparative analysis of global gene expression from Hepa1c1c7 cells and hepatic tissue from C57BL/6 mice identified several model-specific responses to TCDD that should be considered when extrapolating in vitro results to potential in vivo effects. Despite these differences, immortalized cells as well as other emerging in vitro systems (e.g., primary cells, stem cells and 3-D culture systems) provide valuable mechanistic information that supports the further development of high-throughput toxicity screening assays. However, the relevance of in vitro responses requires complementary in vivo verification. Furthermore, comparative studies exploiting other in vitro and in vivo systems, different structurally diverse ligands and other relevant model species will not only corroborate the relevance of the mechanisms, but will also support more appropriate extrapolations between rodent studies and potential effects in humans and ecologically-relevant species.

Methods

Culture and treatment of cell lines

Hepa1c1c7 wild-type and c4 ARNT-deficient cell lines (gifts from O. Hankinson, University of California, Los Angeles, CA) were maintained in phenol-red free DMEM/F12 media (Invitrogen, Carlsbad, CA) supplemented with 5% fetal bovine serum (FBS) (Hyclone, Logan, UT), 2.5 μg/mL amphotericin B (Invitrogen), 2.5 μg/mL amphotericin B (Invitrogen), 50 μg/mL gentamycin (Invitrogen), 100 U/mL penicillin and 100 μg/mL streptomycin (Invitrogen). 1 × 10⁶ cells were seeded into T175 culture flasks (Sarstedt, Newton, NC) and incubated under standard conditions (5% CO₂, 37°C). Time course studies were performed with wild-type and c4 mutant cells where both were dosed with either 10 nM TCDD (provided by S. Safe, Texas A&M University, College Station, TX) or DMSO (Sigma, St. Louis, MO) vehicle and harvested at 1, 2, 4, 8, 12, 24 or 48 hrs. Additional untreated control cells were harvested at the time of dosing (i.e. 0 hrs). For the dose-response study, wild-type cells were treated with DMSO vehicle or 0.001, 0.01, 0.1, 1.0, 10 or 100 nM TCDD and harvested at 12 hrs. The treatment and harvesting regimen for cell culture studies are illustrated in Additional file 7.

Animal treatment

The handling and treatment of female C57BL/6 mice has been previously described [19]. Briefly, immature ovariectomized mice were orally gavaged with 30 μg/kg TCDD for the time course study and sacrificed at 2, 4, 8, 12, 18, 24 72 or 168 hrs after treatment. For the dose-response study, mice were treated with 0.001, 0.01, 0.1, 1, 10, 100 or 300 μg/kg TCDD and sacrificed 24 hrs after dosing. Animals were sacrificed by cervical dislocation and tissue samples were removed, weighed, flash frozen in liquid nitrogen and stored at -80°C until further use.

RNA isolation

Cells were harvested by scraping in 2.0 mL of Trizol Reagent (Invitrogen). Frozen liver samples (approximately 70 mg) were transferred to 1.0 mL of Trizol Reagent and homogenized in a Mixer Mill 300 tissue homogenizer (Retsch, Germany). Total RNA from each study was isolated according to the manufacturer's protocol with an additional acid phenol:chloroform extraction. Isolated RNA was resuspended in The RNA Storage Solution (Ambion Inc., Austin, TX), quantified (A₂₆₀), and assessed for purity by determining the A₂₆₀/A₂₈₀ ratio and by visual inspection of 1.0 μg on a denaturing gel.

Microarray experimental design

Changes in gene expression were assessed using customized cDNA microarrays containing 13,362 features representing 8,284 unique genes. For the time course study, TCDD-treated samples were compared to time-matched vehicle controls using an independent reference design [49]. In this design, treated Hepa1c1c7 cell or hepatic tissue samples were compared to the corresponding time-matched vehicle control with two independent labelings (dye swaps; Additional file 8). Four replicates of this design were performed, each using independent cell culture samples or different animals. Dose-response changes in gene expression were analyzed using a common reference design in which samples from TCDD-treated cells or mice were co-hybridized with a common vehicle reference (i.e. independent DMSO treated Hepa1c1c7 cell samples, hepatic samples from independent sesame oil treated C57BL/6 mice) using two independent labelings (Additional file 8). Four replicates with two independent labelings were performed for both in vitro and in vivo samples. Co-hybridizations of untreated Hepa1c1c7 cells and hepatic tissue from C57BL/6 mice were performed to investigate differences in basal gene expression levels between models (Additional file 8). Four replicates were performed with two independent labelings per sample (dye swap).

More detailed protocols regarding the microarray assay, including microarray preparation, labeling of the cDNA probe, sample hybridization and washing can be obtained from the dbZach website [50]. Briefly, polymerase chain reaction (PCR) amplified cDNAs were robotically arrayed onto epoxy-coated glass slides (Schott-Nexterion, Duryea, PA) using an Omnigrid arrayer (GeneMachines, San Carlos, CA) equipped with 48 (4 × 12) Chipmaker 2 pins (Telechem) at Michigan State University's Research Technology Support Facility [51]. Total RNA (30 μg) was reverse transcribed in the presence of Cy3- or Cy5-deoxyuridine triphosphate (dUTP) to create fluorescence-labeled cDNA, which was purified using a Qiagen PCR kit (Qiagen, Valencia, CA). Cy3 and Cy5 samples were mixed, vacuum dried and resuspended in 48 μL of hybridization buffer (40% formamide, 4× SSC, 1% sodium dodecyl sulfate [SDS]) with 20 μg polydA and 20 μg of mouse COT-1 DNA (Invitrogen) as competitor. This probe mixture was heated at 95°C for 3 min and hybridized on the array under a 22 × 60 mm LifterSlip (Erie Scientific Company, Portsmouth, NH) in a light-protected and humidified hybridization chamber (Corning Inc., Corning, NY) for 18–24 hrs in a 42°C water bath. Slides were then washed, dried by centrifugation and scanned at 635 nm (Cy5) and 532 nm (Cy3) on an Affymetrix 428 Array Scanner (Santa Clara, CA). Images were analyzed for feature and background intensities using GenePix Pro 5.0 (Molecular Devices, Union City, CA).

Microarray data quality assurance, normalization and analysis

Microarray data were first passed through a quality assurance protocol prior to further analysis to ensure consistently high quality data throughout the dose-response and time course studies prior to normalization and further analysis [52]. All the collected data were then normalized using a semi-parametric approach [53]. Empirical Bayes analysis was used to calculate posterior probabilities (P1(t) value) of activity on a per gene and time point or dose group basis using the model-based t-value [54]. The data were filtered using a P1(t) cutoff of 0.9999 and ± 1.5 fold change to identify the most robust changes in gene expression and to obtain an initial subset of differentially regulated genes for further investigation and data interpretation. Subsequent analysis included agglomerative hierarchical and k-means clustering using the standard correlation distance metric implemented in GeneSpring 6.0 (Silicon Genetics, Redwood City, CA). Functional categorization of differentially regulated genes were mined and statistically analyzed from Gene Ontology [55] using GOMiner [56].

Quantitative real-time PCR analysis

For each sample, 1.0 μg of total RNA was reverse transcribed by Superscript II using an anchored oligo-dT primer as described by the manufacturer (Invitrogen). The cDNA (1.0 μL) was used as a template in a 30 μL PCR reaction containing 0.1 μM of forward and reverse gene-specific primers designed using Primer 3 [57], 3 mM MgCl₂, 1.0 mM dNTPs, 0.025 IU AmpliTaq Gold, and 1× SYBR Green PCR buffer (Applied Biosystems, Foster City, CA). PCR amplification was conducted in MicroAmp Optical 96-well reaction plates (Applied Biosystems) on an Applied Biosystems PRISM 7000 Sequence Detection System under the following conditions: initial denaturation and enzyme activation for 10 min at 95°C, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. A dissociation protocol was performed to assess the specificity of the primers and the uniformity of the PCR-generated products. Each plate contained duplicate standards of purified PCR products of known template concentration covering 7 orders of magnitude to interpolate relative template concentrations of the samples from the standard curves of log copy number versus threshold cycle (Ct). No template controls (NTC) were also included on each plate. Samples with a Ct value within 2 standard deviations of the mean Ct values for the NTCs were considered below the limits of detection. The copy number of each unknown sample for each gene was standardized to the geometric mean of three house-keeping genes (β-actin, Gapd and Hprt) to control for differences in RNA loading, quality, and cDNA synthesis. For graphing purposes, the relative expression levels were scaled such that the expression level of the time-matched control group was equal to 1. Statistical analysis was performed with SAS 8.02 (SAS Institute, Cary, NC). Data were analyzed by analysis of variance (ANOVA) followed by Tukey's post hoc test. Differences between treatment groups were considered significant when p < 0.05. Official gene names and symbols, RefSeq and Entrez Gene IDs, forward and reverse primer sequences, and amplicon sizes are listed in Table 4.

Table 4 Gene names and primer sequences for QRTPCR