Transcriptomic responses to 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) in liver: Comparison of rat and mouse

Background Mouse and rat models are mainstays in pharmacology, toxicology and drug development – but differences between strains and between species complicate data interpretation and application to human health. Dioxin-like polyhalogenated aromatic hydrocarbons represent a major class of environmentally and economically relevant toxicants. In mammals dioxin exposure leads to a broad spectrum of adverse affects, including hepatotoxicity of varying severity. Several studies have shown that dioxins extensively alter hepatic mRNA levels. Surprisingly, though, analysis of a limited portion of the transcriptome revealed that rat and mouse responses diverge greatly (Boverhof et al. Toxicol Sci 94:398–416, 2006). Results We employed oligonucleotide arrays to compare the response of 8,125 rat and mouse orthologs. We confirmed that there is limited inter-species overlap in dioxin-responsive genes. Rat-specific and mouse-specific genes are enriched for specific functional groups which differ between species, conceivably accounting for species-specificities in liver histopathology. While no evidence for the involvement of copy-number variation was found, extensive inter-species variation in the transcriptional-regulatory network was identified; Nr2f1 and Fos emerged as candidates to explain species-specific and species-independent responses, respectively. Conclusion Our results suggest that a small core of genes is responsible for mediating the similar features of dioxin hepatotoxicity in rats and mice but non-overlapping pathways are simultaneously at play to result in distinctive histopathological outcomes. The extreme divergence between mouse and rat transcriptomic responses appears to reflect divergent transcriptional-regulatory networks. Taken together, these data suggest that both rat and mouse models should be used to screen the acute hepatotoxic effects of drugs and toxic compounds.


Supplementary Figure 2: Quality Assessment of Rat Arrays
The rat microarray data were subjected to an extensive quality-control procedure prior to conducting downstream analyses. We assessed the spread in intensity distributions before (A) and after (B) GCRMA-pre-processing. We verified that each sample displayed similar RNA degradation characteristics, as indicated by having similar slopes in a plot of Probe intensity vs. position along the transcription (C). Finally we showed that biological replicates clustered together when using the unnormalized data.

Analysis
To assess the variability of the response to TCDD between mice and rats we took the pre-processed and linearly-modeled data and selected all ProbeSets with evidence for differential mRNA abundances in at least one species. We mapped homologs between the two species using the Homologene database. To determine if genes showed similar trends in their profiles we plotted the fold-change in log 2 space (M-values) for all homologs. A) At p adjusted < 0.05 the two profiles are correlated (Spearman's rho = 0.28, p < 2.2 x 10 -16 ), showing similar trends in direction. B) At p adjusted < 0.001 the two profiles are also correlated (Spearman's rho = 0.31, p = 6.12 x 10 -6 ), showing that this finding is independent of the p-value threshold selected.

Supplementary Figure 4: Parameter Sensitivity of Cross-Species Clustering
Some microarray analyses have been described to be sensitive to the selection of specific statistical parameters. For example gene-ontology enrichments may only occur if a gene-list is thresholded at p < 0.01 but not p < 0.05 or p < 0.001. Because threshold selection is ultimately arbitrary, it is important that major conclusions are independent of parameter choices. We demonstrated that mouse and rat orthologs cluster according to TCDD-treatment status ( Figure 2B) by selecting genes that were statistically significant at p adjusted < 0.01. To verify parameter-independence we selected genes at p adjusted < 0.05 (A) and p adjusted < 0.001 (B) and subjected them to the same median-centering, rootmean-square-scaling, and divisive hierarchical clustering as was done previously. In each case the mouse (red) and rat (pink) TCDD-treated animals cluster separately from the mouse (dark blue) and rat (light blue) control animals. Clustering of the genes dimension (B) shows that common-responsive genes (yellow) group together and are few in number relative to rat-specific (red) and mouse-specific (blue) responders.

Supplementary Figure 5: Cross-Species Clustering With Unscaled Data
The use of appropriately normalized data is critical for higher-order statistical analyses.
For example, when we co-clustered rat and mouse data ( Figure 2B) the data were scaled within species prior to executing the pattern-recognition analysis. This step was necessary because the individual ProbeSet-intensities for a given gene can vary between species as a result of the region of the gene targeted, the presence of crosshybridizing transcripts, and the use of different quantile-normalization targets. We explored the effect of omitting this step from the analysis using sets of genes selected at p adjusted < 0.05 (A), p adjusted < 0.01 (B), and p adjusted < 0.001 (C). In each case the TCDD-treated (pink) and control (light blue) rats clustered together, separate from the TCDDtreated (red) and control (blue) mice. Clustering of the genes dimension (C) shows that common-responsive genes (yellow) group together and are few in number relative to ratspecific (red) and mouse-specific (blue) responders.

Supplementary Figure 6: Cross-Species Clustering With Global Scaling
The use of appropriately normalized data is critical for higher-order statistical analyses.
For example, when we co-clustered rat and mouse data ( Figure 2B) the data were scaled within species prior to executing the pattern-recognition analysis. This step was necessary because the individual ProbeSet-intensities for a given gene can vary between species as a result of the region of the gene targeted, the presence of crosshybridizing transcripts, and the use of different quantile-normalization targets. We explored the effect of omitting this step from the analysis (Supplementary Figure 5), and found that removing it eliminated our ability to determine toxicological effects. Next we explored the effect of using a global, per-gene scaling procedure across all animals, rather than by species. Tets of genes were selected at p adjusted < 0.05 (A), p adjusted < 0.01 (B), and p adjusted < 0.001 (C) and subjected to divisive hierarchical clustering, as before.
In each case the TCDD-treated (pink) and control (light blue) rats clustered together, separate from the TCDD-treated (red) and control (blue) mice. Clustering of the genes dimension (C) shows that common-responsive genes (yellow) group together and are few in number relative to rat-specific (red) and mouse-specific (blue) responders.

Supplementary Figure 7: Parameter Sensitivity of Overlap Analysis
Some microarray analyses have been described to be sensitive to the selection of specific statistical parameters. For example gene-ontology enrichments may only occur if a gene-list is thresholded at p < 0.01 but not p < 0.05 or p < 0.001. Because threshold selection is ultimately arbitrary, it is important that major conclusions are independent of parameter choices. We demonstrated that only a small fraction of mouse and rat orthologs respond similarly to TCDD by selecting genes that were statistically significant at p adjusted < 0.01 ( Figure 3A). To verify this conclusion we selected significance thresholds of p adjusted < 0.05 (A) and p adjusted < 0.001 (B). In each case the number of common responses ranges from 4-22% of all TCDD-responsive genes.

Analysis
To test if genes with strong evidence (low p-values) for differential expression in one species would have similarly strong evidence (low p-values) in the other we plotted the adjusted p-values for the two species in log 10 space. A) When we selected genes that were differentially expressed in either species at a moderate-stringency threshold (p adjusted < 0.05) these values are strongly anti-correlated (Spearman's rho = -0.54, p < 2.2 x 10 -16 ). B) When we selected genes that were differentially expressed in either species at a high-stringency threshold (p adjusted < 0.001) we again observed a strong anticorrelation (Spearman's rho = -0.43, p = 5.25 x 10 -11 ). Taken together these data suggest that genes that have strong evidence for differential expression in one species are generally likely to have weak or no evidence in the other, and that this conclusion is independent of the p-value threshold used to select genes.