Abstract
Conventional biochemical and molecular techniques identified previously several genes whose expression is regulated by the aryl hydrocarbon receptor (AHR). We sought to map the complete spectrum of AHR-dependent genes in male adult liver using expression arrays to contrast mRNA profiles in Ahr-null mice (Ahr–/–) with those in mice with wild-type AHR (Ahr+/+). Transcript profiles were determined both in untreated mice and in mice treated 19 h earlier with 1000 μg/kg 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD). Expression of 456 ProbeSets was significantly altered by TCDD in an AHR-dependent manner, including members of the classic AHRE-I gene battery, such as Cyp1a1, Cyp1a2, Cyp1b1, and Nqo1. In the absence of exogenous ligand, AHR status alone affected expression of 392 ProbeSets, suggesting that the AHR has multiple functions in normal physiology. In Ahr–/– mice, only 32 ProbeSets exhibited responses to TCDD, indicating that the AHR is required for virtually all transcriptional responses to dioxin exposure in liver. The flavin-containing monooxygenases, Fmo2 and Fmo3, considered previously to be uninducible, were highly induced by TCDD in an AHR-dependent manner. The estrogen receptor α as well as two estrogen-receptor-related genes (α and γ) exhibit AHR-dependent expression, thereby extending cross-talk opportunities between the intensively studied AHR and estrogen receptor pathways. p53 binding sites are over-represented in genes down-regulated by TCDD, suggesting that TCDD inhibits p53 transcriptional activity. Overall, our study identifies a wide range of genes that depend on the AHR, either for constitutive expression or for response to TCDD.
Initial studies of the aryl hydrocarbon receptor (AHR) focused on its roles in regulating the induction of CYP1 enzymes (Nebert et al., 2004) and mediating toxicity of dioxin-like chemicals (Okey et al., 2005). More recently, the creation of mice with altered AHR signaling revealed phenotypic changes that implicate the AHR in multiple aspects of growth, development, differentiation, and physiology, irrespective of exposure to toxic environmental chemicals (Fernandez-Salguero et al., 1995; Lahvis et al., 2000; Bunger et al., 2003; Walisser et al., 2004a,b). The AHR is a member of the basic helix-loop-helix PAS superfamily and is located in the cytoplasm in association with chaperone proteins such as heat shock protein 90 and XAP2. The AHR is activated by binding to TCDD, translocates to the nucleus, and dimerizes with another basic helix-loop-helix protein, ARNT. The activated AHR/ARNT heterodimer complex interacts with AH-responsive elements and activates the expression of AHR target genes (Nebert et al., 2004).
Mice in which the Ahr gene has been knocked out (Ahr–/–) are extraordinarily resistant to TCDD toxicity (Mimura et al., 1997; Peters et al., 1999; Bunger et al., 2003). Major toxic effects in mice such as hepatic toxicity, thymic atrophy, and cleft palate formation (Bunger et al., 2003) require that the AHR have an intact nuclear translocation/transactivation domain. Rats with a large deletion in the AHR transactivation domain also are highly resistant to lethality from TCDD (Okey et al., 2005). Moreover, mice that are hypomorphic for the AHR's essential dimerization partner, ARNT, are highly resistant to TCDD toxicity (Walisser et al., 2004b). This composite evidence strongly points to AHR-mediated transcriptional effects as the mechanism of dioxin toxicity.
The AHR acts as a ligand-dependent transcription factor, but the spectrum of genes regulated by the AHR is incompletely defined, and the specific genes whose AHR-mediated dysregulation by dioxins leads to major forms of dioxin toxicity are largely unknown. Expression arrays can be highly effective tools for identifying the suite of genes whose expression is altered by xenobiotic chemicals. Several previous array studies examined the responses of various cells or tissues to TCDD or to other xenobiotic chemicals that act as AHR ligands (Puga et al., 2000; Martinez et al., 2002; Guo et al., 2004; Karyala et al., 2004; Boverhof et al., 2005; Fletcher et al., 2005). The majority of these studies were performed in cell culture rather than in whole animals.
Very few of the expression changes in previous studies have been shown to be direct AHR-mediated responses, either by promoter analysis or by determination of AHR-dependence. The possibility remains that some changes in gene expression might be mediated by pathways other than the AHR, or that they might represent a cascade of secondary effects attendant to the onset of dioxin toxicity after long exposure to TCDD. One strategy to assess AHR-dependence is to compare gene expression profiles of cells or tissues that have the wild-type AHR with those of Ahr-null cells in a matched genetic background. This has been done in vitro in two array experiments on aortic smooth muscle cells, derived from Ahr–/– versus Ahr+/+ mice, treated in culture with TCDD (Guo et al., 2004) or benzo[a]pyrene (Karyala et al., 2004) but not previously in tissues harvested from Ahr-null animals treated in vivo. We used expression arrays to identify the batteries of genes whose expression in vivo is affected by the AHR status only, by TCDD only, or by the combination of TCDD and the AHR. This is the first transcriptomic analysis of tissue from Ahr–/– animals.
Materials and Methods
TCDD. TCDD was purchased from the UFA-Oil Institute (Ufa, Russia) and was >99% pure as determined by gas chromatography-mass spectrometry.
Animals and In Vivo Treatment. Male Ahr–/– mice (10 weeks old) in a C57BL/6J background were obtained from The Jackson Laboratory (Bar Harbor, ME). Wild-type (Ahr+/+) male C57BL/6J mice (15 weeks old) were bred at the National Public Health Institute (Kuopio, Finland) from stock originally obtained from The Jackson Laboratory. Mice were given a single dose of 1000 μg/kg TCDD or corn oil vehicle by gavage. TCDD initially was dissolved in ether and added to corn oil; the ether subsequently was evaporated off. Liver was harvested 19 h after treatment, sliced, snap-frozen, and stored in liquid nitrogen until homogenization. There were three TCDD-treated and three control mice in the Ahr–/– groups and six TCDD-treated and five control mice in the Ahr+/+ groups (Fig. 1).
RNA Isolation and Expression Array Studies. Total RNA was extracted using RNeasy kits (QIAGEN, Valencia, CA) according to the manufacturer's instructions. DNase (QIAGEN) was added to the RNeasy elution column as recommended by the manufacturer. RNA yield was quantified by UV spectrophotometry, and RNA integrity was verified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). The isolated RNA was then assayed on Affymetrix MOE430-2 arrays (Affymetrix, Santa Clara, CA) at The Centre for Applied Genomics at The Hospital for Sick Children (Toronto, Canada) following standard manufacturer's protocols.
Preprocessing of Array Data. Array data were loaded into the R statistical package (version 2.0.0) using the affy package (version 1.5.8; http://www.bioconductor.org) (Gautier et al., 2004) of the BioConductor open-source project (Gentleman et al., 2004). Data were investigated for spatial and distributional homogeneity, then preprocessed with a sequence-specific version of the RMA algorithm (Irizarry et al., 2003) termed GCRMA, as implemented in BioConductor (version 1.1.3). The data then were written to disk and parsed with Perl scripts into a custom-built Oracle database partially derived from the MAGE-OM-compliant RAD schema (Stoeckert et al., 2001). The normalized data were then associated with updated sequence annotation as detailed elsewhere (P. C. Boutros, N. Tijet, R. Pohjanvirta, and A. B. Okey, manuuscript in preparation). In brief, the annotation process involved comprehensive alignment of each target sequence to species-specific expressed sequence tag databases, followed by a poll-based cluster-assignment algorithm to determine the best matching to UniGene clusters (Wheeler et al., 2004).
Significance Testing of Array Data. Normalized, annotated data were written out from the Oracle database and significance-tested with a general-linear model using the limma package (version 1.8.14) in BioConductor. The following linear model was fit to each individual ProbeSet on the MOE430-2 array: Y = Basal + AHR + TCDD + AHR-TCDD + Batch. Here, Y refers to the expression level of a single ProbeSet; Basal refers to the underlying basal expression level across all animals; AHR captures “AHR-dependent, TCDD-independent” expression changes; TCDD captures “TCDD-dependent, AHR-independent” expression changes; AHR-TCDD captures “AHR-dependent, TCDD-dependent” expression changes; and Batch captures the effects of the hybridization batch. After fitting the linear model, we used empirical Bayes moderation of the standard error (Smyth et al., 2003), followed by a false-discovery rate correction for multiple testing (Efron and Tibshirani, 2002). To identify differentially expressed genes, we used a nested F test as implemented in limma on the AHR, TCDD, and AHR-TCDD effects. This test considers first whether a ProbeSet displays any differential expression and then assesses which contrasts contribute to that differential expression. ProbeSets were deemed differentially expressed at the p = 10–3 significance level. For each effect, we extracted estimates of differential expression (“M”) for all ProbeSets and parsed these values back into the Oracle database for storage and downstream processing. Complete output files are available in Supplemental Table S2.
Experimental design for high-throughput transcriptional profiling of the effects of the AHR, of TCDD, and their interaction. Wild-type C57BL/6 and Ahr–/– mice in a C57BL/6 genetic background were treated either with 1000 μg/kg TCDD or corn-oil vehicle (control). The number of independent animals studied, each on a separate array, is given in brackets for each group. The three effects fitted to a general linear model are illustrated. The AHR-effect refers to “AHR-dependent, TCDD-independent” changes in expression profiles. The TCDD-effect refers to “AHR-independent, TCDD-dependent” changes in expression profiles. The AHR-TCDD interaction captures “AHR-dependent, TCDD-dependent” changes in expression profiles.
Transcription-Factor Binding Site Analysis. To understand the regulatory networks underlying AHR-dependent expression, we performed two separate transcription-factor binding-site (TFBS) analyses. Both analyses used the University of California at Santa Cruz's Genome Browser Database build mm5 of the mouse genome, from which regulatory sequences were extracted using BioPerl-based scripts (Stajich et al., 2002). Transcriptional start sites (TSSs) were taken from the University of California at Santa Cruz's Genome Browser Database (Karolchik et al., 2003). For some genes, the annotated TSSs are uncertain or in error, including Cyp1a1 (Sun et al., 2004); for Cyp1a1, we used the manually annotated TSS of Chr9:57,902,425, courtesy of Dr. Patricia Harper. To identify the genomic region encoding a specific ProbeSet, we mapped the UniGene clusters annotated for each ProbeSet to the RefSeq annotation using the LocusLink-UniGene mapping made available by the National Center for Biotechnology Information (Wheeler et al., 2004).
First, to identify the known AHR binding sites, we scanned the region from –5000 to +1000 bp relative to each TSS for AHRE-I and AHRE-II binding sites, as described previously (Boutros et al., 2004).
Second, to characterize the patterns of transcription factor binding sites across the set of putative AHR-regulated genes, we used a library-based TFBS search. These library-based searches compare a set of known TFBSs with the promoter regions of a set of coregulated genes and then apply statistical tests to identify TFBSs that are found more often (enriched) or less often (depleted) than expected by chance in the set of coregulated genes. Using the Clover software (Frith et al., 2004), we tested separately the sets of genes up- and down-regulated in both the AHR and AHR-TCDD contrasts. Statistical testing used both mononucleotide and dinucleotide randomizations of the parental sequences (1000 iterations), randomization of scoring matrices, and two separate background sequence sets (murine CpG islands and human chromosome 21). Two separate libraries of TFBSs were scanned: the JASPAR database (Sandelin et al., 2004), and a custom-designed database that included three variants of the AHRE-I consensus sequence, the AHRE-II sequence, the AnoC element, and a Cyp1a1 negative regulatory element. This search scanned the region from –4000 to +1000 bp relative to the annotated TSS for each gene. The sequences are the following (Boutros et al., 2004): AHRE-1(Core) = GCTCG; AHRE-I(Extended) = TNGCGTG; AHRE-I(Full) = [T|G]NGCGTG[A|C][G|C]A; AHRE-II = CATG{N6}C[T|A]TG; AnoC = [C|T]GCG[C|T]GCGC[C|A]GC; and Cyp1a1 negative regulatory element = G[G|T]GCTCTG[G|C][G|C][G|A][G|A][T|A]CA[G|A][G|A][G|T]C[C|A][C|A].
Real-Time qRT-PCR. mRNA levels for selected genes were quantified by real-time qRT-PCR. See Supplemental Table S1 for primer and probe sequences and details of the procedure.
Functional Analysis of Differentially Expressed Genes. Ontological analysis used build 140 of the GO-Miner software package (Zeeberg et al., 2003), which uses the Fisher's exact test to identify significantly enriched functional categories (GO terms) in the lists of differentially expressed genes.
Chromosomal Location of Differentially Expressed Genes. To identify spatial patterns of gene expression across the genome, we developed a novel spatial-search algorithm. For each chromosome, a 400,000-bp region of genomic sequence was checked for two criteria: 1) at least four genes in this region were present on the array and 2) there was at least one differentially expressed gene in the effect being considered. If these criteria were met, then the probability of statistically significant enrichment in this region was calculated from the hypergeometric distribution. The window was then stepped in 200,000-bp increments across the chromosome. At the end of the process, p values were adjusted with a Bonferroni correction for multiple-testing, and significant spatial clustering was reported at the p < 0.05 level. To help visualize these patterns, we plotted the genomic locations of all genes that were responsive to the AHR or AHR-TCDD effects using the Geneplotter package of the BioConductor open-source project (Gentleman et al., 2004).
Results
Genotype Confirmation. The Ahr–/– mice we used originally were generated by replacing exon 2 with the neomycin-resistance gene (Schmidt et al., 1996). Before array studies, we confirmed the Ahr genotypes of all mice using primers localized in exon 1 and exon 4. As expected, the amplified product from mRNA of each Ahr–/– mouse differed in size from that in the Ahr+/+ wild-type mice (Supplemental Fig. S1).
Overall Patterns of Gene Expression Detected by Arrays. To obtain the maximum information and usefulness from our array experiments, we used a two-factor, two-level design: the factors were genotype (Ahr+/+ or Ahr–/–) and treatment (TCDD-treated or control). As shown in Fig. 1, we characterized three distinct sets of genes: 1) those affected by Ahr genotype independent of TCDD; 2) those responsive to TCDD in an AHR-dependent manner; and 3) those responsive to TCDD in an AHR-independent manner. The linear modeling approach we used increases the statistical power of the study by using all 17 animals in the estimation of each of the three contrasts rather than using multiple pairwise comparisons, each involving fewer animals. From the Venn diagram (Fig. 2), it is clear that the AHR is the critical factor governing changes in gene expression. More than 90% of ProbeSets that exhibited differential expression do so solely in response to either the AHR genotype (244/706 = 35%) or the AHR genotype in combination with TCDD administration (430/706 = 61%). Furthermore, there is a strong directionality to the effect of TCDD. Fully 75% of TCDD-responsive ProbeSets are up-regulated (Fig. 3).
AHR Effects Independent of TCDD. Most studies of the AHR's role in gene expression have focused on AHR-mediated transcriptional enhancement in response to foreign agonists such as TCDD. However, there is a large suite of genes whose hepatic expression is clearly affected by the AHR status independent of exposure to any exogenous AHR ligand. Table 1 and Fig. 2 reveal a total of 392 ProbeSets whose expression levels differ significantly between Ahr–/– and Ahr+/+ mice independent of TCDD treatment. The AHR seems to influence hepatic expression of nearly as many ProbeSets in the absence of exogenous ligands (392) as it does after TCDD administration (456) (Table 1). Although not all of these ProbeSets are necessarily regulated directly by the AHR, the impact of AHR status on constitutive expression of numerous genes suggests multiple and diverse roles for the AHR in normal physiology in addition to the AHR's ability to mediate gene expression in response to xenobiotic ligands. We harvested tissue 19 h after TCDD administration; thus, some changes in gene expression could be AHR-dependent but be secondary, downstream responses to other AHR-dependent changes rather than being primary responses directly regulated by the AHR.
Number of ProbeSets whose expression is significantly affected by the AHR, by TCDD, or by the AHR-TCDD interaction
Wild-type Ahr+/+ mice and Ahr-/- mice were treated, and array experiments were performed as described under Materials and Methods. Cells in the table indicate the number ProbeSets whose expression was affected at a significance level of p < 0.001.
Overlap of effects. ProbeSets were identified as differentially expressed by first normalizing the raw expression data, then fitting a linear model to the three effects, followed by Bayesian moderation of standard error and a false-discovery rate correction. ProbeSets were deemed significantly differentially expressed by a nested F test at the p < 0.001 level, which identified 706 ProbeSets from among the 45,101 on the array. The Venn diagram shows the overlap between effects in terms of the number of ProbeSets.
The magnitude of the AHR effect on some genes was very large. For example, in untreated mice, the mRNA for Serpina12, a proteinase inhibitor, was 220-fold more abundant in liver from Ahr+/+ mice than from Ahr–/– mice (Table 2 and Fig. 4). Constitutive expression of numerous other genes also was substantially higher in Ahr+/+ mice than in Ahr–/– mice (Table 2).
AHR effects independent of TCDD
The 25 ProbeSets that are most up-regulated and the 25 most down-regulated are shown. Both transient (UniGene) and stable (LocusLink) gene identifiers are given. ProbeSets that carry no LocusLink could be mapped to a UniGene cluster, but that cluster did not correspond to a RefSeq gene. ProbeSets that carry no name and no UniGene cluster identifier could not be mapped to any cluster of transcribed sequences and may represent repetitive or nonspecific sequences. In some cases, the array contains multiple ProbeSets that map to the same gene; thus, the same gene name can appear at multiple places in the table. In cases in which genomic sequence is available, the number of AHRE-I and AHRE-II elements found in the region -5000 to +1000 bp relative to the TSS are given; see Materials and Methods for the sequence of each element. “M” values indicate the magnitude in log2-space of each contrast. For example, an effect of -3 refers to an 8-fold reduction in expression, whereas an effect of +2 refers to a 4-fold increase in mRNA levels. The abbreviated names of genes that are discussed in the text are shown in boldface type in parentheses. A full list of all responsive genes is available in Supplemental Table S2.
Cyp1a2 represents a gene whose expression is well known to be inducible by TCDD via the AHR (Nebert et al., 2004). In addition to its AHR-regulated inducibility by TCDD, the Cyp1a2 gene also displayed significantly higher expression (4-fold) in Ahr+/+ mice than in Ahr–/– mice, independent of TCDD administration (Table 2); higher constitutive hepatic CYP1A2 mRNA expression in Ahr+/+ mice than in Ahr–/– mice was previously observed (Shimada et al., 2002). It is possible that some endogenous ligand interacts with the AHR to support “constitutive” expression of Cyp1a2 in mice bearing the wild-type AHR. However, if such an endogenous ligand exists, it provokes a selective response, because in the absence of an exogenous ligand, the expression of other prototypical AHR-regulated genes such as Cyp1a1 and Cyp1b1 was not higher in Ahr+/+ mice than in Ahr–/– mice. Hectd2, a gene recently reported to be induced by AHR ligands (Hayes et al., 2005), behaved very much like Cyp1a2; that is, in addition to being inducible by TCDD in an AHR-dependent fashion, constitutive Hectd2 expression also was substantially higher in mice with wild-type AHR than in Ahr–/– mice (Table 2).
Directionality of effect overlap. Differentially expressed ProbeSets were classified by the direction of change and show the profound directionality of the three effects. The AHR effect is evenly distributed between up- and down-regulation (44% of ProbeSets up-regulated). In contrast, the TCDD effect is skewed toward up-regulation (75% of ProbeSets up-regulated) as is the AHR-TCDD interaction (74% up-regulated). Note: ProbeSets regulated in opposite directions will not be represented in an intersection of the Venn diagrams; hence, the values in Figs. 2 and 3 are not additive on a subset-by-subset basis.
Effects of the AHR on hepatic expression in the absence of exogenous ligand are not confined to up-regulation (i.e., higher expression in Ahr+/+ mice than in Ahr–/– mice). In fact, presence of the AHR is associated with lowered expression of significantly more genes than it “up-regulates” (Table 1 and Fig. 3). Two cytochrome P450 mRNAs—CYP17A1 and CYP2B20—seem to be dramatically suppressed by the presence of the AHR in the absence of any exogenous ligand (Table 2 and Fig. 4). In the absence of treatment, mRNA levels for metallothionein 1 and metallothionein 2 are 40- to 64-fold lower in Ahr+/+ mice than in Ahr–/– mice (Table 2 and Fig. 4); mechanistically, this may not be caused by direct suppression of metallothionein gene expression by the AHR (see Discussion).
AHR-Dependent Effects of TCDD. The expression levels of 456 ProbeSets were altered by TCDD in an AHR-dependent manner; that is, these ProbeSets showed differential responses to TCDD between Ahr+/+ versus Ahr–/– mice (Table 1 and Fig. 3). There is substantial overlap between the TCDD-dependent and the -independent effects of the AHR (145/703 ProbeSets differentially expressed in both effects). These effects usually are in accord (125/145 common ProbeSets); it is particularly rare for the AHR to control the up-regulation of a gene's basal levels but the down-regulation of its response to TCDD (3/145 common probe sets). Basal down-regulation combined with TCDD-induced up-regulation is somewhat more common (17/145 common ProbeSets) and includes genes such as Cyp2a4 and the known TCDD-responsive Gstm3.
To assess the validity of the array data and the overall treatment and sample preparation, we examined the expression of genes that are well known to be up-regulated by TCDD in animals that have wild-type AHR. Classic AHR-regulated genes such as Cyp1a1, Cyp1b1 (Nebert et al., 2004), Nqo1, and Tiparp displayed robust AHR-dependent induction by TCDD (Table 3), confirming that the experimental protocol was sound. However, outside of the prototypical AH gene battery, there is only a moderate overlap in responsive genes detected in our high-throughput study compared with several other array studies of response to AHR ligands in varied tissues and cell types. There also is little consistency among the sets of genes identified across the previous array studies themselves (data not shown).
Responses of ProbeSets that represent genes in the classic AH gene battery to AHR status or to AHR-TCDD interaction
Data are derived from array experiments as described under Materials and Methods. “AHR” represents the effect of having wild-type AHR (versus Ahr-/-), whereas “AHR-TCDD” represents the effect of TCDD treatment in mice with wild-type AHR (Ahr+/+ mice). As in Table 2, “M” represents the magnitude in log2-space for each contrast. There are multiple entries in the table for some genes because the array contains multiple ProbeSets for these genes.
Several genes that were not known previously to respond to AHR ligands demonstrated AHR-dependent induction by TCDD of 20-fold or more. The genes whose expression was most highly up-regulated (Table 4) represent a variety of functional categories, including the flavin-containing monooxygenases Fmo2 (30-fold) and Fmo3 (80-fold), two neuronal proteins; and Purkinje cell protein 4-like (Pcp4l1, 80-fold), neuronal pentraxin 1 (Nptx1; 80-fold), and tubulin α8 (Tuba8; 30-fold).
AHR-dependent effects of TCDD
The 25 ProbeSets that are most up-regulated and the 24 most down-regulated are shown. Both transient (UniGene) and stable(LocusLink) gene identifiers are given. ProbeSets that carry no LocusLink could be mapped to a UniGene cluster, but that cluster did not correspond to a RefSeq gene. ProbeSets that carry no name and no UniGene cluster identifier could not be mapped to any cluster of transcribed sequences and may represent repetitive or nonspecific sequences. In some cases, the array contains multiple ProbeSets that map to the same gene; thus, the same gene name can appear at multiple places in the table. In cases in which genomic sequence is available, the number of AHRE-I and AHRE-II elements found in the region -5000 to +1000 bp relative to the TSS are given; see Materials and Methods for the sequence of each element. “M” values indicate the magnitude in log2-space of each contrast. For example, an effect of -3 refers to an 8-fold reduction in expression, whereas an effect of +2 refers to a 4-fold increase in mRNA levels. The abbreviated names of genes that are discussed in the text are shown in boldface type in parentheses. A full list of all responsive genes is available in Supplemental Table S2.
Confirmation of mRNA expression levels by real-time qRT-PCR. RNA used in the array studies was assayed for levels of mRNA for a selected set of the most responsive genes by real-time qRT-PCR as described in Supplemental Data. Bars represent the mean of three animals for Ahr–/– mice and of five or six animals for mice with wild-type AHR. Data are normalized so that for each gene, the treatment group with the highest expression is set at 100%, and expression in other groups is relative to this. Error bars represent S.D. from the mean. RIKEN represents RIKEN gene 4432416J03Rik.
As shown in Table 1, the AHR-TCDD interaction encompasses approximately 3-fold more ProbeSets that are induced than are repressed. Nevertheless, the expression of many genes was strongly suppressed by TCDD in an AHR-dependent manner (Table 4). The mRNA whose expression was down-regulated to the greatest extent (20-fold) in an AHR-dependent fashion is the proteinase inhibitor Serpina7. Down-regulation of Serpina7 by TCDD is in sharp contrast to the response of Serpina1, which is highly induced by TCDD (Table 4), and Serpina12, whose constitutive expression is dramatically higher in Ahr+/+ mice than in Ahr–/– mice (Table 2) and which is moderately down-regulated by TCDD in an AHR-dependent fashion (Table 4 and Fig. 4).
TCDD Has Few AHR-Independent Effects. One of the most striking findings of our study is that very few genes respond to TCDD in the absence of the AHR. Only 32 genes showed statistically significant AHR-independent responses to TCDD, and the maximum difference in expression between TCDD-treated and untreated mice in the absence of the AHR was approximately 2-fold (Table 5). It was anticipated that TCDD would mediate gene expression primarily through the AHR. Our broad interrogation of the mouse transcriptome confirms that almost all transcriptomic effects of TCDD do indeed require the AHR.
Complete list of all ProbeSets that displayed an AHR-independent response to TCDD
The right-hand column indicates the AHR-independent effect of TCDD. As shown in the other two columns, some of the genes that exhibited an AHR-independent effect of TCDD also displayed a TCDD-independent AHR-effect (left-hand column) or an AHR-TCDD interaction effect (middle column). ProbeSets that carry no LocusLink could be mapped to a UniGene cluster, but that cluster did not correspond to a RefSeq gene. ProbeSets that carry no name and no UniGene cluster identifier could not be mapped to any cluster of transcribed sequences and may represent repetitive or nonspecific sequences. In some cases, the array contains multiple ProbeSets that map to the same gene; thus, the same gene name can appear at multiple places in the table. In cases in which genomic sequence is available, the number of AHRE-I and AHRE-II elements found in the region -5000 to +1000 bp relative to the TSS are given; see Materials and Methods for the sequence of each element. “M” values indicate the magnitude in log2-space of each contrast. For example, an effect of -3 refers to an 8-fold reduction in expression, whereas an effect of +2 refers to a 4-fold increase in mRNA levels. The abbreviated names of genes that are discussed in the text are shown in boldface type in parentheses. A full list of all responsive genes is available in Supplemental Table S2.
Confirmation of mRNA Levels by Real-Time Quantitative RT-PCR. Array methods for measuring mRNA levels are becoming highly reliable. Nevertheless, we believed it prudent to use real-time qRT-PCR as an independent method to evaluate mRNA levels for selected genes whose expression, by array analyses, seemed to be related to AHR status and/or TCDD treatment. We used Cyp1a1 as a benchmark because CYP1A1 mRNA is known to be inducible via the AHR pathway. As expected, the array results and the qRT-PCR results showed that CYP1A1 mRNA is strongly up-regulated by TCDD in Ahr+/+ mice, whereas CYP1A1 mRNA is undetectable in Ahr–/– mice or in untreated Ahr+/+ mice (Fig. 4, top row). Nptx1 and Fmo3 behave very much like CYP1A1: each of these mRNAs is highly up-regulated by TCDD in Ahr+/+ mice, whereas in Ahr–/– mice mRNA levels are very low in either the constitutive state or in TCDD-treated animals (Fig. 4; 2nd row). Our results show a strong AHR dependence for induction of Nptx1 and Fmo3 (Table 4 and Fig. 4).
Real-time qRT-PCR also confirms array results which show that mRNA levels for metallothionein 1 (Mt1), CYP17A1, and CYP2B20 are dramatically higher in Ahr–/– mice than in Ahr+/+ mice (Fig. 4, 3rd row); none of these genes shows a significant response to TCDD in either Ahr+/+ mice or Ahr–/– mice.
For Tuba8, constitutive expression is much lower in Ahr+/+ mice than in Ahr–/– mice. TCDD treatment in Ahr+/+ mice essentially “restores” Tuba8 mRNA to levels similar to those in Ahr–/– mice (Fig. 4, row 4). It is possible that the AHR suppresses Tuba8 when no exogenous ligand is present and that TCDD interferes with this suppression.
Functional Categorization of Responsive Genes by Ontological Analysis. High-throughput expression experiments often identify large numbers of genes whose expression is altered by biological status or by the treatment of interest. It can be difficult to discern functionally relevant patterns by simple inspection of large data sets. To characterize the gene lists, we used ontological analysis to identify functional categories that were enriched in each of the effects studied: AHR (Table 6), AHR-TCDD (Table 7), and TCDD (Supplemental Table S3).
Ontological analysis of the AHR effect on transcript expression
Ontological analysis used the GOMiner software package (Zeeberg et al., 2003), which performs a two-tailed Fisher's exact test to determine over-representation of GO categories in the gene list relative to the array as a whole. The GO term and name are given, along with the total number of genes in the category and the numbers of over- and underexpressed genes that have this annotation. p values, shown in parentheses, are calculated for enrichment of the term as a whole and for enrichment of the term among over- and underexpressed genes alone. GO categories shown were selected from those showing statistical significance (P < 0.001) in the AHR effect.
Ontological analysis of the AHR-TCDD interaction on transcript expression
Ontological analysis used the GOMiner software package (Zeeberg et al., 2003), which performs a two-tailed Fisher's exact test to determine over-representation of GO categories in the gene list relative to the array as a whole. The GO term and name are given, along with the total number of genes in the category and the numbers of over- and underexpressed genes that have this annotation. p values, shown in parentheses, are calculated for enrichment of the term as a whole and for enrichment of the term among over- and underexpressed genes alone. GO categories shown were selected from those showing statistical significance (P < 0.001) in the AHR effect. GO categories shown were selected from those showing statistical significance (p < 0.001) in the AHR-TCDD effect.
The AHR effect is quite functionally coherent (Table 6). It includes genes for numerous proteins involved in electron transport, drug metabolism, and defense against xenobiotic toxicants, including 13 cytochrome P450s, 4 glutathione-S-transferases, 2 thioredoxins, and 2 metallothioneins. It also includes three separate aquaporins (water channels), which may indicate a perturbation in osmotic homeostasis within the cell, consistent with the large number of ion channels that also show altered expression. Three of the eight genes on the array that have roles in aromatic amino acid catabolism are down-regulated, and three of nine negative regulators of cell adhesion are up-regulated.
The AHR-TCDD interaction (Table 7) shows the clearest patterns. There is profound up-regulation of genes involved in electron transport, detoxification, and ribosome assembly and structure. There also is a strong down-regulation of genes involved in amino acid metabolism in general, and particularly, catabolism.
The AHR-independent TCDD effect (Supplemental Table S4) displayed little functional coherency—genes that are affected by TCDD in the absence of the AHR do not fit any particular category.
Transcription Factor Binding Sites. Classic TCDD-inducible/AHR-dependent genes such as Cyp1a1 contain an AHRE-I motif that binds the ligand-AHR-ARNT complex and leads to enhanced transcription. A second response element, AHRE-II, has been discovered that seems to mediate the induction of rat CYP1A2 (Sogawa et al., 2004) and several other genes (Boutros et al., 2004) by AHR agonists. We searched for AHRE-I and AHRE-II motifs in the 5′-flanking sequence of all genes whose expression was significantly influenced by the AHR or TCDD in the array experiment.
The core AHRE-I sequence (GCGTG) is found in nearly all responding genes listed in Tables 2, 3, 4, 5. The majority of these core sequences probably represent random occurrence of this pentanucleotide which will occur approximately once every 256 bases (45 for two strands and two orientations) throughout the genome. Neither the core AHRE-I or AHRE-II binding sites nor any of the longer AHRE-I sequences that were examined in an earlier phylogenetic footprinting study (Boutros et al., 2004) showed significant enrichment in any condition. This is consistent with the lack of substantial overlap between the genes identified in previous microarray analyses and in our current direct test of AHR activation. It is possible that the 19-h exposure to TCDD in our array study resulted in some secondary responses that could obscure enrichment in both AHRE-I and AHRE-II binding sites. This hypothesis is supported by an enrichment of AHRE-I motifs in the 25 ProbeSets most highly induced by TCDD in an AHR-dependent fashion (AHR-TCDD up-regulated list). This suggests that the AHRE-I motif is solely involved in TCDD-dependent induction through the AHR, rather than repression or ligand-independent induction. The presence of confounding secondary effects would be particularly relevant for the AHRE-II because AHRE-II was found to be generally depleted in liver-expressed genes.
To characterize patterns of other transcription-factor binding sites in addition to AHRE-I and AHRE-II motifs, we used a library-based TFBS search. Genes were separated according to their direction of change for both the AHR effect and the AHR-TCDD effect, yielding four lists of coexpressed genes. We identified two basic groups of results. First, a series of TFBSs were found to be consistently enriched or depleted in all four sets of coexpressed genes tested (Table 8). These motifs seem to represent liver-specific transcription factors, such as HNF3. When such motifs are under-represented, they may represent nonhepatic transcription factors (e.g., AHRE-II) or transcription factors antagonized by the AHR. The latter class might include E2F motif, in which the AHR can corepress E2F-dependent expression through interactions with rubidium (Marlowe et al., 2004). Second, several TFBSs were found to be enriched or depleted in only one or two of the sets of coexpressed genes (Table 9). These motifs may represent particular mechanisms of AHR regulation. For instance, genes repressed by TCDD through the AHR are enriched for p53 binding motifs, indicating a potential cross-talk between these pathways.
Transcription factor binding sites common to all conditions
p values for transcription factor binding sites that were found to be over-represented in all four of the effects studied. Here, AHR(↓) refers to genes with negative AHR effects, whereas AHR-TCDD(↑) refers to genes with positive interactions. The ID number refers to the JASPAR identifier of the matrix used for testing the specific TFBSs. The p values given here were obtained by dinucleotide randomization studies (1000 iterations). Of the 15 TFBSs here, only five are under-represented (MA0054, MA0103, MA0036, and MA0024), whereas 10 are over-represented. Many of these TFBSs are likely to be liver-specific rather than AHR-specific.
AHR-specific transcription factor binding sites identified
A selection of statistically significant TFBSs identified as depleted (i.e., HLF) or enriched (all others) in at least one contrast. ↓ refers to genes with negative effects, whereas ↑ refers to genes with positive effects. p values were derived from 1000 iterations of dinucleotide randomization. The enrichment of SP1 sites in genes induced by TCDD in an AHR-dependent fashion (such as Cyp1a1) is supportive of the known direct interactions between the AHR and SP1. On the other hand, p53-regulated genes generally are proapoptotic/cell-defense genes; thus, the TCDD-dependent down-regulation of such genes is surprising and may contribute to some of the acute toxic effects of TCDD exposure.
Chromosomal Locations of Responsive Genes. It is possible that genes which respond to a particular regulatory molecule such as the AHR might occur in clusters or “hotspots” in the genome (Cohen et al., 2000). We mapped the chromosomal location for named genes as depicted in Fig. 5. A sliding-window analysis identified some hotspots of differential expression. One notable hotspot for the AHR effect is on chromosome 8, where carboxylesterase-2 and -5 are in close proximity to a RIKEN gene (2210023G05Rik); this enrichment is statistically significant at p < 10–3. The AHR-TCDD effect is much more spatially clustered than the AHR effect. One notable cluster for AHR-TCDD involves genes encoding four glutathione S-transferase proteins on chromosome 3 (Gstm1, Gstm2, Gstm3, and Gstm4); this group of glutathione S-transferase genes may have arisen by gene duplication and remain in proximity on chromosome 3. Another cluster, on chromosome 6, involves three unrelated genes: a solute carrier (Slc6a13), a protease (Usp18), and α-tubulin (Tuba8) with p < 10–3.
Chromosomal distribution of genes that respond to the AHR or to the AHR-TCDD interaction. Genes that exhibited a statistically significant response to AHR or to the AHR-TCDD interaction were mapped to the genome using the geneplotter package of the BioConductor open-source project. The chromosomes that contain genes encoding the AHR, ARNT, and prototypical AHR-responsive genes (Cyp1a1, Cyp1a2, Cyp1b1, Cyp2s1, Nqo1, Ugt1a6, and Ahrr) are indicated. Red bars represent positions of genes that are up-regulated, and blue bars indicate down-regulated genes. White bars represent the locations of all other genes contained on that chromosome that are not affected by AHR status or by AHR-TCDD interaction.
Although there are some significant clusters, it is clear from data in Fig. 5 that genes whose expression is affected by the AHR in a ligand-dependent or -independent manner are widely dispersed across the mouse genome, with the exception of sex chromosomes. The sex chromosomes contain few genes that seem to be regulated by the AHR, at least in liver; there are no hepatically expressed AHR-regulated genes on the Y chromosome and only a handful on the X chromosome. Although the X chromosome generally is gene-rich, it shows minimal differential expression related to AHR, with nonrandom under-representation for both for the AHR effect (p = 0.004) and the AHR-TCDD interaction (p = 0.002).
Discussion
Our array experiments in Ahr–/– versus Ahr+/+ mice are intended to comprehensively define the separate sets of genes whose expression is significantly affected by AHR status per se or by the interaction of TCDD with AHR in mouse liver. This is a prelude to the longer-term goal of determining which AHR-regulated genes are central to normal physiology/development and which genes are dysregulated by dioxins in a manner that evokes toxicity. Our functional analyses (Tables 6 and 7) reveal that AHR status and TCDD exposure influence the expression of genes that carry out a wide variety of functions.
Genes Related to Reproduction, Growth, and Development.Ahr–/– mice exhibit reduced fecundity (Lahvis et al., 2000). The reasons are not clear, but female reproduction is susceptible to disruption by alterations in estrogen function. One ongoing challenge in dioxin research is to mechanistically rationalize the extensive cross-talk between estrogen-signaling pathways and dioxin-signaling pathways (Beischlag and Perdew, 2005). Our study shows that expression of both Esr1 (ERα) itself and two of its closely related proteins, Esrra and Esrrg, can be regulated by the AHR. With the extensive intertwining and antagonism of estrogen and AHR pathways, it seems that they are reciprocally regulated. That is, the AHR modulates Esr1 levels, thereby potentially dampening estrogen effects. We postulate that basal AHR levels are, in turn, partly regulated by estrogens through Esr1 or one of the Esr-related receptors. In the absence of TCDD, hepatic Esr1 levels are 3.5-fold lower in Ahr+/+ than in Ahr–/– mice. The estrogen-receptor-related α gene (Esrra) also is moderately down-regulated by the AHR in a TCDD-independent manner (Supplemental Table S2). Moreover, our data show that estrogen-receptor-related γ (Esrrg) is dramatically induced (11-fold) by TCDD in an AHR-dependent fashion (Table 4). The function of Esrrg is largely unknown. Esrrg does not bind estrogenic ligands but can homodimerize and occupy a wide range of estrogen response elements, potentially antagonizing the action of the major ligand-dependent estrogen receptors, Esr1 (ERα) and Esr2 (ERβ) (Razzaque et al., 2004). Thus, the high degree of Esrrg induction by TCDD constitutes a novel mechanism by which TCDD interferes with estrogen signaling in development and reproduction.
Our study also reveals that the constitutive level of CYP17A1 mRNA is 30-fold higher in Ahr–/– mice than in Ahr+/+ mice (Table 2 and Fig. 4). CYP17 enzymes are essential for biosynthesis of both glucocorticoids and sex steroid precursors. It is possible that excessive hepatic CYP17A1 in livers of Ahr–/– mice disrupts homeostasis of reproductive steroids and impairs fecundity, because this enzyme has 17α-hydroxylase activity and 17,20-lyase activity, as well as activity as a squalene monooxygenase, which can alter cholesterol biosynthesis. TCDD or the dioxin-like PCB126 has been reported to decrease CYP17A expression in rat testis (Mebus et al., 1987), but our results show that hepatic CYP17A1 expression is unaltered by TCDD exposure for 19 h, regardless of Ahr genotype.
Genes Related to Cell Cycle, Cell Growth, Differentiation, and Apoptosis. Evidence continues to expand that the AHR plays fundamental roles in the cell cycle and apoptosis (Hanlon et al., 2005). We found neuronal pentraxin I (Nptx1) to be induced more than 60-fold in livers of Ahr+/+ mice treated with TCDD. In the nervous system, Nptx1 is induced by hypoxia and is involved in apoptotic cell death in neurons undergoing stress. Nptx1 is expressed in liver but its function there is unclear. Induction of Nptx1 by TCDD may indicate that Nptx1 is broadly responsive to a wide variety of environmental stressors.
Proteinase Inhibitors. Proteinase inhibitors play major roles in tissue maintenance by reducing degradation of the extracellular matrix by enzymes such as matrix metalloproteinases (MMPs). Up-regulation of proteinase inhibitors in TCDD-exposed tissues may be beneficial because TCDD can induce MMPs. We found that MMP-24 was up-regulated 2-fold by TCDD in livers of Ahr+/+ mice (Supplemental Table S2), but this was accompanied by a 30-fold induction of the proteinase inhibitor Serpine1 (Table 4). Our experiments revealed that Serpina12 mRNA is dramatically reduced in livers of Ahr–/– mice compared with Ahr+/+ mice, but the effect of TCDD is modest (Table 2 and Fig. 4). It is possible that some of the liver pathology that afflicts Ahr–/– mice (Walisser et al., 2004b) is related to deficiency of Serpina12.
Genes Related to Toxicity or to Defense from Toxicity of Xenobiotic Chemicals. Our array experiment shows that TCDD elicits the expected very high induction of mRNAs for CYP1A1, CYP1A2, CYP1B1, and Nqo1 in Ahr+/+ mice along with induction of the conjugating enzymes Gstm1 and Ugt1a6. We also found (Table 4) a remarkably high induction of the flavin monooxygenase enzymes Fmo2 and Fmo3. Flavin monooxygenase enzymes play important roles in detoxification of many foreign chemicals, including psychoactive drugs, pesticides, and dietary-derived compounds, but previously were believed not to be inducible in mammals (Cashman, 2002). Flavin monooxygenase enzyme induction may be a further adaptive mechanism by which the AHR fosters protection from xenobiotic chemicals.
Constitutive PEPCK (Pck1) expression was much higher in livers of Ahr–/– mice than in livers of Ahr+/+ mice. PEPCK is a key enzyme in gluconeogenesis, and studies in rodents have shown that acute toxicity of TCDD may be associated with reduced PEPCK activity, together with appetite suppression and inhibition of gluconeogenesis (Viluksela et al., 1999). We found that TCDD treatment decreased Pck1 mRNA levels in Ahr+/+ mice but not in Ahr–/– mice (Fig. 4), lending support to the potential AHR-dependent role of PEPCK in dioxin toxicity.
Another striking finding in our array experiments was the profoundly lower expression of mRNAs for the metallothioneins Mt1 and Mt2 in Ahr+/+ mice than in Ahr–/– mice (Table 2). Metallothioneins are considered a major defense mechanism against DNA damage and apoptosis (Cherian et al., 2003). We speculate that the high constitutive levels of Mt1 and Mt2 in Ahr–/– mice reflect a state of oxidative stress in livers of these animals compared with their wild-type counterparts caused by an insufficiency of other protective mechanisms against oxidative stress.
Our functional analysis revealed broad-scale dysregulation of genes involved in protein synthesis. A total of 19 ribosomal transcripts had their expression moderately up-regulated by TCDD in an AHR-dependent manner, and eight are involved in ribosome biogenesis, suggesting that TCDD is stimulating the assembly of new ribosomes. Expression of 12 genes involved in amino acid catabolism is altered by TCDD in an AHR-dependent manner, and fully eight of these are down-regulated. Furthermore, we see that Eif3s6 (formerly Int-6), a critical member of the EIF3 translation initiation complex and its interacting protein, Eif3s6ip, both are up-regulated by TCDD in an AHR-dependent manner. These two interacting proteins mediate cross-talk between the translation initiation complex, the COP9 signalosome, and the 26S proteasome. These observations suggest that TCDD initiates a broad cellular effort to increase protein translation by increasing ribosome assembly and decreasing amino acid catabolism, possibly through Eif3-mediated changes in protein degradation rates. In fact, increased plasma amino acid levels have been observed in TCDD-treated rats (Viluksela et al., 1999).
AHR-Mediated Alteration of Gene Expression: Possible Mechanisms for Up-Regulation and Down-Regulation. As shown in Tables 2, 4, and 5, the core sequence for AHRE-I is “available” in nearly all genes that exhibited up-regulation in response to AHR status or to the AHR-TCDD interaction. The “full” AHRE-I sequence is found in 12 of the 25 genes that were most highly up-regulated by the AHR-TCDD interaction but in only one of the genes that were most greatly down-regulated (Table 4).
The mere presence of an AHRE motif is not sufficient to establish that this motif is responsible for regulating the expression of a particular gene. In addition to AHRE motifs, we mapped other TFBSs and found a set of regulatory elements common to all four gene lists (Table 8) that includes 10 enriched and five depleted TFBSs. The identity of several of the enriched TFBSs strongly supports the hypothesis that these are general liver-specific promoter elements (e.g., HFH-3). An alternative mechanism has emerged by which the AHR or TCDD might down-regulate gene expression. Ray and Swanson (2004) reported that TCDD suppresses the expression of p53 and p16INK4a by provoking methylation of promoters in these genes in an AHR-dependent fashion. In our study, p53 binding sites are over-represented in genes down-regulated by the TCDD in an AHR-dependent fashion (Table 9).
To sum up, our key findings are that the AHR affects the expression of large, separable gene batteries in the presence or absence of TCDD. Furthermore, virtually all effects of TCDD on gene expression require the AHR. Key dysregulated pathways include estrogen signaling, energy metabolism, and protein synthesis. Almost as many genes are responsive to AHR status per se as are responsive to the AHR-TCDD interaction, which supports the concept that the AHR plays important roles in normal development and physiology, not just in response to toxic environmental chemicals. Array studies provide menus of genes that are plausible candidates to be involved in particular biological or toxicological outcomes. The challenge now is to link AHR-mediated expression of specific genes and gene batteries to those phenotypic outcomes.
Acknowledgments
We thank the BioConductor user-support group for outstanding support, Dr. Gordon K. Smyth for assistance in fitting linear models, Dr. Rafal Kustra for assistance in formulating the spatial-assessment algorithm, Virpi Tiihonen and Arja Tamminen for excellent technical assistance, and Dr. Philip Seeman for use of the real-time polymerase chain reaction instrument.
Footnotes
- Received September 7, 2005.
- Accepted October 7, 2005.
Supported by grant MOP-57903 from the Canadian Institutes of Health Research (to A.B.O.) and by grants 200980 and 211120 from the Academy of Finland (to R.P.). I.D.M. was supported by a fellowship from the Natural Sciences and Engineering Research Council of Canada.
N.T. and P.C.B. contributed equally to this work.
Article, publication date, and citation information can be found at http://molpharm.aspetjournals.org.
doi:10.1124/mol.105.018705.
ABBREVIATIONS: AHR, aryl hydrocarbon receptor; ARNT, aryl hydrocarbon receptor nuclear translocator; AHRE-I, aryl hydrocarbon response element (also known as DRE or XRE-I); AHRE-II, aryl hydrocarbon response element-II (also known as XRE-II); GO, gene ontology; TFBS, transcription factor binding site; TCDD, 2,3,7,8-tetrachlorodibenzo-p-dioxin; TSS, transcription start site; bp, base pair; Mt1, metallothionein 1; ER, estrogen receptor; MMP, metalloproteinase; qRT-PCR, quantitative reverse transcription-polymerase chain reaction.
↵
The online version of this article (available at http://molpharm.aspetjournals.org) contains supplemental material.
- The American Society for Pharmacology and Experimental Therapeutics
