Abstract

The number of legacy chemicals without toxicity reference values combined with the rate of new chemical development is overwhelming the capacity of the traditional risk assessment paradigm. More efficient approaches are needed to quantitatively estimate chemical risks. In this study, rats were dosed orally with multiple doses of six chemicals for 5 days and 2, 4, and 13 weeks. Target organs were analyzed for traditional histological and organ weight changes and transcriptional changes using microarrays. Histological and organ weight changes in this study and the tumor incidences in the original cancer bioassays were analyzed using benchmark dose (BMD) methods to identify noncancer and cancer points of departure. The dose-response changes in gene expression were also analyzed using BMD methods and the responses grouped based on signaling pathways. A comparison of transcriptional BMD values for the most sensitive pathway with BMD values for the noncancer and cancer apical endpoints showed a high degree of correlation at all time points. When the analysis included data from an earlier study with eight additional chemicals, transcriptional BMD values for the most sensitive pathway were significantly correlated with noncancer (r = 0.827, p = 0.0031) and cancer-related (r = 0.940, p = 0.0002) BMD values at 13 weeks. The average ratio of apical-to-transcriptional BMD values was less than two, suggesting that for the current chemicals, transcriptional perturbation did not occur at significantly lower doses than apical responses. Based on our results, we propose a practical framework for application of transcriptomic data to chemical risk assessment.

The process for developing a chemical risk assessment is time and resource intensive. Depending on the nature of the endpoints evaluated, data generation, literature review, data analysis, and review and publication of the assessment can take many years to complete. As a result, only 553 substances have published toxicological reviews in the U.S. Environmental Protection Agency’s (EPA) Integrated Risk Information System (IRIS), and 613 pesticides have been evaluated by the Office of Pesticides Programs (OPP) from 1991 through 2008. By comparison, 2000–2500 premanufacturing notifications for new chemicals are reviewed by the EPA each year (EPA, 1997). These numbers provide a clear disparity between chemicals with published toxicity values and new submissions. The public health consequences of using untested chemicals can be significant. For polluted sites across the United States, it is not uncommon to find that contaminants identified in potential exposure media have no published toxicity values. These chemicals cannot be considered quantitatively in the overall hazard index calculation (EPA, 1989). Toxicologically uncharacterized chemicals are essentially treated as not posing any potential health hazard, which could lead to an underestimation of human health risk.

In the absence of available human data, a typical EPA noncancer risk assessment incorporates repeat-dose toxicity information from animal studies of various exposure durations (e.g., short term, subchronic, and chronic). The amenable dose-response data are fit to a series of statistical models, and a benchmark dose (BMD) and the associated lower confidence bound (BMDL) are estimated using the best-fitting model. Uncertainty factors are then applied to the BMDL value (i.e., the point of departure) to take into account potential differences in cross-species extrapolation, pharmacodynamic and pharmacokinetic variability, sensitive subpopulations, and exposure duration. The end result is a reference value, which is defined as an estimate of an exposure to the human population (including susceptible subgroups) that is likely to be without an appreciable risk of adverse health effects over a lifetime (EPA, 2002).

A similar approach is undertaken for cancer-related endpoints. In the absence of human data, tumor incidence and multiplicity are evaluated from chronic rodent cancer bioassays. The bioassays are typically 2 years in duration and performed at two or three dose levels. A BMD and associated BMDL are estimated for the amount of chemical causing a specified increase in tumor incidence, e.g., 10% (EPA 2005). Both the cancer slope factor and the amount of the chemical corresponding to risk levels in humans (e.g., from 1 in 10,000 to 1 in 1,000,000) are estimated based on linear extrapolation from the BMDL. The inherent limitation for performing such analyses is the availability of in-life toxicity data that inform potential toxicological hazards and associated dose response.

The application of transcriptomic data in risk assessment has primarily focused on informing the mode of action of a chemical (EPA, 2009) and using the information as part of a weight of evidence. In a previous study, histological, organ weight, and transcriptomic changes were evaluated in dose response following a subchronic exposure to five chemicals (Thomas et al., 2011). The apical (i.e., histological and organ weight) and transcriptional changes were analyzed using a BMD approach, and the results showed good concordance between the apical and transcriptional BMD values when the transcriptional changes were summarized based on either biological processes (Thomas et al., 2011) or signaling pathways (Thomas et al., 2012). Many of the signaling pathways that correlated with the apical responses could be used to inform key events because they have been implicated in noncancer and cancer disease pathogenesis.

Although transcriptomic data have potential utility for defining mode of action in a weight of evidence, the use of transcriptomic data in this context has not yet addressed the need for more rapid risk assessments on larger numbers of chemicals. In this study, rats were orally exposed to multiple doses of six chemicals for 5 days and 2, 4, or 13 weeks. The six chemicals were selected because they had previously published toxicological reviews in the IRIS database, including results of cancer bioassays or subchronic repeat-dose toxicity evaluations. The target tissues were analyzed for traditional histological and organ weight changes and transcriptional changes using microarrays. The dose-response changes in gene expression were analyzed using a BMD approach and the responses grouped based on signaling pathways. The results were analyzed in two ways. The first examined the relative sensitivity and the overall correlation between the BMD values for cancer or noncancer apical effects of a chemical and the transcriptional BMD values for the most sensitive signaling pathway. The second evaluated the temporal stability of the correlation between apical and transcriptional BMD values to assess the duration of exposure required to estimate the point of departure. The results showed that transcriptional changes in the most sensitive pathway can be used to estimate both a cancer and noncancer BMD value for application to chemical risk assessment. The transcriptional changes can be assessed at any time between 5 days and 13 weeks.

MATERIALS AND METHODS

Chemicals.

1,2,4-Tribromobenzene (TRBZ; CAS No. 615-54-3; Purity 95%; Catalog No. 1617321) and 2,3,4,6-tetrachlorophenol (TTCP; CAS No. 58-90-2; Purity 98%; Catalog No. 1618611) were purchased from International Laboratories USA (San Bruno, CA). Bromobenzene (BRBZ; CAS No. 108-86-1; Purity 99%; Catalog No. AC10668-0010) and 4,4′-Methylenebis(N,N-dimethyl) benzenamine (MDMB; CAS No. 101-61-1; Purity 99+%; Catalog No. AC12674-1000) were purchased from Acros Organics through Fisher Scientific (Pittsburg, PA). Hydrazobenzene (HZBZ; CAS No. 122-66-7; Catalog No. 126721) was purchased from Sigma-Aldrich (St Louis, MO). N-Nitrosodiphenylamine (NDPA; CAS No. 86-30-6; Purity 97+%; Catalog No. D0899) was purchased from TCI America through VWR Scientific (Radnor, PA). The corn oil (Catalog No. 700000–134) and olive oil (Catalog No. 95034–796) vehicles were purchased from VWR Scientific.

Animals and treatment.

The animal exposures for all the chemicals have been described previously (Dodd et al., 2012a,, b,, c,, d,, e,, f). For TRBZ and TTCP, 4- to 5-week-old male Sprague Dawley (Rat Crl:CD(SD)) rats were purchased from Charles River Laboratories (Raleigh, NC). For BRBZ and HZBZ, 4- to 5-week-old male F344/DuCrl rats were purchased from Charles River Laboratories (Kingston, NY). For MDMB and NDPA, 4- to 5-week-old female F344/DuCrl rats were purchased from Charles River Laboratories. The chosen strain, sex, and route of exposure were those that displayed the critical effect in the previous IRIS toxicological reviews for these chemicals. Upon arrival, the rats were acclimated for 2 weeks prior to initiation of the study. The rats were weighed and randomized into treatment groups to ensure the mean body weight in each group was approximately the same. Group sizes were 10 rats per concentration per time point. An additional 1 or 2 rats were assigned to the control group per time point when available to be used to assure a group size of at least 10 for evaluation of biological end points. Animals were ear tagged and housed two or fewer rats per cage in shoebox-style cages separated by treatment group. Alpha-dri cellulose bedding (Shepard Specialty Papers, Kalamazoo, MI) was used. Animals had access to reverse osmosis water (Hydro Systems, Durham, NC) and pellet NIH-07 certified feed (Zeigler Brothers, Gardners, PA) ad libitum. The animal room was kept within the standard temperature and humidity parameters (64°F–79°F with 30–70% relative humidity) and standard light cycle (0700–1900h).

Animal exposures were initiated at 5–6 weeks of age and were performed by the route and dose listed in Table 1. For each chemical, animals were exposed at five dose levels for 5 days and 2, 4, or 13 weeks. A matched vehicle control was run concurrently with each exposure. The doses for the exposures overlapped those used in the IRIS toxicological reviews. Gavage exposures were administered 7 days per week in a volume of 5ml/kg of corn oil or olive oil. Exposures in the feed were administered continuously 7 days per week.

TABLE 1

Chemical Treatment Details and Rodent Model Used in the Study

Chemical Abbreviation Dosesa Rodent model Target tissue 
1,2,4-Tribromobenzene TRBZ 2.5, 5, 10, 25, and 75mg/kgb Male Sprague Dawley rats Liver 
Bromobenzene BRBZ 25, 100, 200, 300, and 400mg/kgb Male F344 rats Liver 
2,3,4,6-Tetrachlorophenol TTCP 10, 25, 50, 100, and 200mg/kgb Male Sprague Dawley rats Liver 
4,4′-Methylenebis (N,N-dimethyl) benzenamine MDMB 50, 200, 375, 500, and 750 ppmc Male F344 rats Bladderd 
N-Nitrosodiphenylamine NDPA 250, 1000, 2000, 3000, and 4000 ppmc Female F344 rats Thyroidd 
Hydrazobenzene HZBZ 5, 20, 80, 200, and 300 ppmc Male F344 rats Liverd 
Chemical Abbreviation Dosesa Rodent model Target tissue 
1,2,4-Tribromobenzene TRBZ 2.5, 5, 10, 25, and 75mg/kgb Male Sprague Dawley rats Liver 
Bromobenzene BRBZ 25, 100, 200, 300, and 400mg/kgb Male F344 rats Liver 
2,3,4,6-Tetrachlorophenol TTCP 10, 25, 50, 100, and 200mg/kgb Male Sprague Dawley rats Liver 
4,4′-Methylenebis (N,N-dimethyl) benzenamine MDMB 50, 200, 375, 500, and 750 ppmc Male F344 rats Bladderd 
N-Nitrosodiphenylamine NDPA 250, 1000, 2000, 3000, and 4000 ppmc Female F344 rats Thyroidd 
Hydrazobenzene HZBZ 5, 20, 80, 200, and 300 ppmc Male F344 rats Liverd 

Note. aUnderlined doses used in previous rodent subchronic or chronic studies.

bGavage dosing, corn oil vehicle, 5ml/kg (BRBZ and TRBZ); gavage dosing, olive oil vehicle, 5 ml/kg (TTCP).

cFeed dosing.

dHave rodent cancer bioassay data.

Animal use in this study was approved by the Institutional Animal Use and Care Committee of The Hamner Institutes for Health Sciences and was conducted in accordance with the National Institutes of Health guidelines for the care and use of laboratory animals. Animals were housed in fully accredited American Association for Accreditation of Laboratory Animal Care (AAALAC) facilities.

Necropsy and histology.

Animal necropsies occurred at the specified time points. For the gavage exposures, necropsies occurred within a few hours following the last gavage administration. Animals were weighed and anesthetized with a lethal ip injection of sodium pentobarbital. A cardiac puncture was performed to collect blood samples, and the animal was then exsanguinated via transection of the abdominal aorta. Following gross examination for abnormalities, the target organs were removed. For the liver, slices from each of the lobes were placed in RNAlater (Ambion, Austin, TX), separate slices were flash frozen in liquid nitrogen, and slices were preserved in 10% neutral buffered formalin. For the thyroid, the gland was transversely sectioned. The top portion was preserved in 10% neutral buffered formalin, and the bottom portion was minced and placed in RNAlater. For the bladder, the organ was longitudinally bisected, a slice was placed in 10% neutral buffered formalin, and the remainder of the organ was minced and placed in RNAlater. The formalin-fixed tissues were embedded in paraffin, sectioned at 5 µm, and stained with hematoxylin and eosin. Histological changes were assessed by an accredited pathologist and are described in detail in previous publications (Dodd et al., 2012a, b, c, d, e, f). Typically, 10 rats were evaluated per concentration per time point.

Gene expression microarray measurements.

Total RNA was isolated from six rats per dose per time point using either flash frozen (TRBZ, TTCP, and HZBZ) or RNAlater preserved tissues (BRBZ, MDMB, and NDPA). Tissues were homogenized in either TRIzol (LifeTechnologies, Carlsbad, CA)(BRBZ) or RLT Lysis Buffer (Qiagen, Valencia, CA)(TRBZ, TTCP, MDMB, NDPA, and HZBZ). After homogenization, the samples were stored at −80°C until the remainder of the RNA isolation procedure was performed. Samples homogenized in TRIzol were thawed, centrifuged to remove debris, and loaded onto Phase Lock Gel tubes (5 PRIME, Gaithersburg, MD). The Phase Lock Gel tubes were centrifuged, and the aqueous phase was used for RNA purification. The samples that were homogenized in RLT buffer were centrifuged to remove debris, and the supernatant was used for RNA purification. RNA purification was performed using a QIAcube (Qiagen) with the standard RNA Protocol and RNeasy mini kits.

Following purification, the quantity of RNA was determined spectometrically and the integrity of the RNA was evaluated with the Agilent 2100 Bioanalyzer (Palo Alto, CA). The RNA from the five animals with the highest RNA integrity numbers were used for microarray analysis. Double-stranded cDNA was synthesized from 150ng of total RNA, and biotin-labeled cRNA was transcribed using the IVT Express Kit (Affymetrix, Santa Clara, CA). A total of 15 µg of labeled cRNA was fragmented, and a hybridization cocktail was prepared, plated, and loaded into an Affymetrix GeneTitan system. Hybridization was performed for 16h on HT Rat230+PM microarrays. The GeneTitan system was used to perform the microarray hybridization, washing, and scanning. The gene expression results have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (Accession No.: GSE45892).

Transcriptional BMD analysis.

The analysis of the gene expression microarray data was performed as described previously (Thomas et al., 2007) with some modification. The analysis was performed using BMDExpress (Version 1.41)(Yang et al., 2007). The microarray data were processed using Robust Multi-Array Average normalization with a log2 transformation (Irizarry et al., 2003). The normalized intensity values for each probe set were fit as continuous data to a series of four different dose-response models—linear, 2° polynomial, 3° polynomial, and power models. Each model was run assuming constant variance, and the benchmark response (BMR) factor was set to 1.349 multiplied by the SD in the control animals to estimate a BMD with a 10% increase in tail area (Thomas et al., 2007). For model selection, a nested likelihood ratio test was performed on the linear, 2° polynomial, and 3° polynomial models. If the more complex model provided a significantly improved fit (p < 0.05), the more complex model was selected. If the more complex model did not provide a significantly improved fit (p ≥ 0.05), the simpler model was selected (Posada and Buckley, 2004). The Akaike information criterion (AIC) for the selected polynomial model was then compared with the AIC for the power model. The model with the lowest AIC (Akaike, 1973) was selected as the final, best fitting model and was used to calculate a BMD and BMDL. To avoid model extrapolation and any potential bias from probe sets with poorly fitting models, probe sets with a BMD value greater than the highest dose or a goodness-of-fit p value < 0.1 were removed from further analysis.

Affymetrix probe sets were converted into unique genes based on their NCBI Entrez Gene ID. Promiscuous probe sets (i.e., those that interrogate more than one gene) were removed from the analysis. When two or more probe sets were associated with a single gene, the BMD and BMDL values for the individual probe sets were averaged to obtain single BMD and BMDL values. The Entrez Gene identifiers were then matched to their corresponding pathways using the GeneGo Metacore database (Minibase version 6.11.41105_7884, GeneGo, St Joseph, MI). The database consisted of a total of 744 pathways. Pathways with fewer than five genes with BMD value less than the highest dose and a goodness-of-fit p value > 0.1 were removed from the analysis. Unless otherwise stated, median BMD and BMDL values were used to summarize each pathway. BMD and BMDL values for MDMB, NDPA, and HZBZ in ppm were converted to milligrams per kilogram-day using strain-specific subchronic food intake factors calculated based on recommended biological values from the EPA (1988).

Gene expression data from our previous study (National Center for Biotechnology Information Gene Expression Omnibus Accession Nos: GSE18858 and GSE17933) (Thomas et al., 2011) were also reanalyzed using the updated transcriptional BMD methods and GeneGo Metacore database.

BMD analysis of noncancer endpoints.

The organ weight data and treatment-related histological data associated with these studies have been published previously (Dodd et al., 2012a,, b,, c,, d,, e,, f). The organ weight data (absolute and relative) were fit as continuous data to a series of seven models (exponential model 2, exponential model 3, exponential model 4, exponential model 5, linear, polynomial 5°, and power; BMDS software, Version 2.1.1, U.S. EPA). A BMR of 1 SD was used to estimate the BMD and BMDL. An adequate model fit was judged based on the goodness-of-fit p value (p ≥ 0.1), magnitude of the scale residuals in the vicinity of the BMR, and visual inspection of the model fit. In addition to these three criteria for judging adequacy of model fit, a determination was made as to whether the variance across dose groups was homogeneous. Among all the models providing adequate fit to the data set, the model with the lowest AIC was selected.

The treatment-related histological changes were fit as dichotomous data to a series of eight models (gamma, logistic, log logistic, probit, log probit, multistage with 5° polynomial, Weibull, and Hill; BMDS software, Version 2.1.1, U.S. EPA). All goodness-of-fit assessments were evaluated by a likelihood ratio-based p value (compared with the saturated model) and AIC comparison. A BMR of 10% extra risk was used to estimate the BMD and BMDL. BMD and BMDL values for MDMB, NDPA, and HZBZ in ppm were converted to milligrams per kilogram-day using strain-specific subchronic food intake factors calculated based on recommended biological values from the EPA (1988).

BMD analysis of cancer apical endpoints.

The thyroid follicular cell carcinoma and adenoma incidences (combined) for MDMB (NTP, 1979a), the incidences of transitional cell carcinoma of the bladder for NDPA (NTP, 1979b), and the incidences of hepatocellular carcinomas and neoplastic liver nodules for HZBZ (NTP, 1978) from the 2-year rodent cancer bioassays were fit using the multistage cancer model with a BMR of 10% extra risk (BMDS software, Version 2.1.1, U.S. EPA). Goodness of fit was assessed by a likelihood ratio-based p value (compared with the saturated model) and AIC comparison. BMD and BMDL values in ppm were converted to milligrams per kilogram-day using strain-specific chronic EPA food intake factors (EPA, 1988).

Correlation analysis.

The correlation between the apical and transcriptional BMD values was assessed using a standard Pearson’s correlation coefficient. When adjusting for the maximum tolerated dose (MTD), a partial correlation was calculated. Calculations of the correlation coefficients and the partial correlation coefficients were performed in JMP (Version 9.0.0, SAS Institute, Inc., Cary, NC). The statistical significance of the correlation coefficients was assessed using a standard t distribution (Steel and Torrie, 1980).

RESULTS

BMD Analysis of Noncancer Endpoints

The changes in the absolute and relative organ weights and histology were used to obtain BMD values for noncancer effects in the target tissues at each time point. The detailed organ weight and histological changes for each chemical are provided in Supplementary files 1 and 2. From all the histological and organ weight changes for each chemical, the most sensitive adverse response was selected. For TRBZ and BRBZ, changes in absolute liver weight, inflammation, and hepatocytic hypertrophy were the most sensitive adverse apical responses (Table 2). For TRBZ, the BMD values declined over time beginning with 12.4mg/kg-day at 5 days and dropping to 3.8mg/kg-day at 13 weeks. For BRBZ, the BMD values declined from 191mg/kg-day at 5 days to 72mg/kg-day at 4 weeks and then increased to 96.8mg/kg-day at 13 weeks. For TTCP, absolute liver weight and hepatocytic vacuolation were the most sensitive adverse apical responses. Similar to TRBZ, the BMD values declined throughout the study from 91.6mg/kg-day at 5 days to 1.0mg/kg-day at 13 weeks. For MDMB, follicular cell hypertrophy was the most sensitive adverse apical response at all time points. The BMD values for MDMB declined from 22.1mg/kg-day at 5 days to 1.7mg/kg-day at 4 weeks and then increased back to 11.0mg/kg-day at 13 weeks. For NDPA, a change in absolute bladder weight was the most sensitive adverse response at 5 days, whereas diffuse hyperplasia of the transitional epithelium was the most sensitive endpoint at the remaining time points. The BMD values were similar at 5 days and 2 weeks, dropped at 4 weeks, and remained steady out to 13 weeks. Finally, pathological changes related to HZBZ treatment were observed (e.g., hypertrophy, bile duct duplication, and macrovesiculation), but the constellation of effects was not sufficient at any time or dose to be considered adverse (Dodd et al., 2012f)(Supplementary file 2).

TABLE 2

Noncancer Point-of-Departure Values for the Most Sensitive Apical Response

Chemical Target tissue Endpoint BMD (BMDL) (mg/kg-day or ppm)a Route-adjusted BMD (BMDL) (mg/kg-day)b 
5 Days 
TRBZ Liver Absolute liver weight 12.4 (5.5) — 
BRBZ Liver Inflammation 191 (159) — 
TTCP Liver Absolute liver weight 91.6 (53.7) — 
MDMB Thyroid Follicular cell hypertrophy 221 (137) 22.1 (13.7) 
NDPA Bladder Absolute bladder weight 2051 (1614) 232 (182) 
HZBZ Liver None — — 
2 Weeks 
TRBZ Liver Hepatocytic hypertrophy 9.3 (6.0) — 
BRBZ Liver Absolute liver weight 130 (27.2) — 
TTCP Liver Absolute liver weight 13.6 (8.0) — 
MDMB Thyroid Follicular cell hypertrophy 157 (48.4) 15.7 (4.8) 
NDPA Bladder Diffuse transitional epithelial hyperplasia 2220 (1806) 251 (204) 
HZBZ Liver None — — 
4 Weeks 
TRBZ Liver Hepatocytic hypertrophy 5.1 (3.4) — 
BRBZ Liver Absolute liver weight 72 (40.2) — 
TTCP Liver Hepatocytic vacuolation 6.7 (2.3) — 
MDMB Thyroid Follicular cell hypertrophy 17.0 (3.1) 1.7 (0.3) 
NDPA Bladder Diffuse transitional epithelial hyperplasia 1689 (1053) 191 (119) 
HZBZ Liver None — — 
13 Weeks 
TRBZ Liver Absolute liver weight 3.8 (2.0) — 
BRBZ Liver Absolute liver weight 96.8 (61.4) — 
TTCP Liver Hepatocytic vacuolation 1.0 (0.6) — 
MDMB Thyroid Follicular cell hypertrophy 110 (37.9) 11.0 (3.8) 
NDPA Bladder Diffuse transitional epithelial hyperplasia 1567 (971) 177 (110) 
HZBZ Liver None — — 
Chemical Target tissue Endpoint BMD (BMDL) (mg/kg-day or ppm)a Route-adjusted BMD (BMDL) (mg/kg-day)b 
5 Days 
TRBZ Liver Absolute liver weight 12.4 (5.5) — 
BRBZ Liver Inflammation 191 (159) — 
TTCP Liver Absolute liver weight 91.6 (53.7) — 
MDMB Thyroid Follicular cell hypertrophy 221 (137) 22.1 (13.7) 
NDPA Bladder Absolute bladder weight 2051 (1614) 232 (182) 
HZBZ Liver None — — 
2 Weeks 
TRBZ Liver Hepatocytic hypertrophy 9.3 (6.0) — 
BRBZ Liver Absolute liver weight 130 (27.2) — 
TTCP Liver Absolute liver weight 13.6 (8.0) — 
MDMB Thyroid Follicular cell hypertrophy 157 (48.4) 15.7 (4.8) 
NDPA Bladder Diffuse transitional epithelial hyperplasia 2220 (1806) 251 (204) 
HZBZ Liver None — — 
4 Weeks 
TRBZ Liver Hepatocytic hypertrophy 5.1 (3.4) — 
BRBZ Liver Absolute liver weight 72 (40.2) — 
TTCP Liver Hepatocytic vacuolation 6.7 (2.3) — 
MDMB Thyroid Follicular cell hypertrophy 17.0 (3.1) 1.7 (0.3) 
NDPA Bladder Diffuse transitional epithelial hyperplasia 1689 (1053) 191 (119) 
HZBZ Liver None — — 
13 Weeks 
TRBZ Liver Absolute liver weight 3.8 (2.0) — 
BRBZ Liver Absolute liver weight 96.8 (61.4) — 
TTCP Liver Hepatocytic vacuolation 1.0 (0.6) — 
MDMB Thyroid Follicular cell hypertrophy 110 (37.9) 11.0 (3.8) 
NDPA Bladder Diffuse transitional epithelial hyperplasia 1567 (971) 177 (110) 
HZBZ Liver None — — 

Note. aBMD (BMDL).

bRoute-adjusted BMD and BMDL values for MDMB, NDPA, and HZBZ in ppm were converted to millligrams per kilogram-day using strain-specific subchronic EPA food intake factors (EPA, 1988). Male F344 rats, 0.100; Female F344 rats, 0.113.

BMD Analysis of Cancer Endpoints

The incidences of thyroid follicular cell carcinoma and adenoma (combined) for MDMB (NTP, 1979a), the incidences of transitional cell carcinoma of the bladder for NDPA (NTP, 1979b), and the incidences of hepatocellular carcinomas and neoplastic liver nodules for HZBZ (NTP, 1978) from the 2-year rodent cancer bioassays were analyzed to obtain the BMD values corresponding to 10% extra risk (Table 3). In the original bioassays, all three chemicals were run with two doses and a control. For MDMB and NDPA, the fourth order multistage cancer model produced the lowest AIC value and was selected as the best-fitting model. The goodness-of-fit p value was adequate for both models. For HZBZ, a second order multistage cancer model was selected. However, a goodness-of-fit p value could not be calculated because the model used all three parameters to fit the three dose levels.

TABLE 3

Cancer Point-of-Departure Estimates for the Tumor Responses in the 2-Year Rodent Cancer Bioassay

Chemical Model degree of polynomiala Fit p value BMD (BMDL) (ppm)b Route-adjusted BMD (BMDL) (mg/kg-d)c 
MDMB 0.98 381 (283) 30.1 (22.4) 
NDPA 0.85 2002 (1499) 184 (138) 
HZBZ NAd 38.1 (21.1) 3.0 (1.7) 
Chemical Model degree of polynomiala Fit p value BMD (BMDL) (ppm)b Route-adjusted BMD (BMDL) (mg/kg-d)c 
MDMB 0.98 381 (283) 30.1 (22.4) 
NDPA 0.85 2002 (1499) 184 (138) 
HZBZ NAd 38.1 (21.1) 3.0 (1.7) 

Note. aTumor incidence was fit to a multistage cancer model with the degree of polynomial listed in the table.

bBMD, dose at 10% extra risk; BMDL, 95% lower bound on BMD.

cRoute-adjusted BMD and BMDL values in ppm were converted to milligrams per kilogram-day using strain-specific chronic EPA food intake factors (EPA, 1988). Male F344 rats, 0.079; Female F344 rats, 0.092.

dThe “NA” for the p value indicates that the model used all three of their parameters to fit.

BMD Analysis for Transcriptional Changes

The BMD values for each gene were grouped based on canonical signaling pathways and summarized as the median of each pathway. The distribution of pathway transcriptional BMD values provides an indication of the dose at which the majority of signaling pathways become transcriptionally active. The distributions were plotted as a function of exposure duration (Fig. 1). For four of the chemicals (BRBZ, TTCP, NAPD, and HZBZ), the pathway transcriptional BMD values decreased at 2 weeks and then increased at 4 and 13 weeks. For TRBZ, the distributions of pathway transcriptional BMD values generally decreased over time. Conversely, the distributions of pathway transcriptional BMD values increased for MDMB. The full set of pathway transcriptional BMD analysis results is provided in Supplementary files 3–5.

FIG. 1.

Box-and-whisker plots of the median transcriptional BMD values across the signaling pathways for each chemical and exposure duration. (A) TRBZ, (B) BRBZ, (C) TTCP, (D) MDMB, (E) NDPA, and (F) HZBZ. The horizontal lines within the box depict the medians, and the lower and upper edges of the boxes represent the 25th and 75th percentiles (i.e., interquartile range). The diamond in the box contains the mean and the upper and lower 95% confidence interval for the mean. The whiskers represent the range of values 1.5 times the interquartile range. The bracket outside of the box identifies the shortest half, which is the most dense 50% of the observations. A table containing the descriptive statistics for each time point is included below each panel. UQ, upper quartile; LQ, lower quartile; Min, minimum.

FIG. 1.

Box-and-whisker plots of the median transcriptional BMD values across the signaling pathways for each chemical and exposure duration. (A) TRBZ, (B) BRBZ, (C) TTCP, (D) MDMB, (E) NDPA, and (F) HZBZ. The horizontal lines within the box depict the medians, and the lower and upper edges of the boxes represent the 25th and 75th percentiles (i.e., interquartile range). The diamond in the box contains the mean and the upper and lower 95% confidence interval for the mean. The whiskers represent the range of values 1.5 times the interquartile range. The bracket outside of the box identifies the shortest half, which is the most dense 50% of the observations. A table containing the descriptive statistics for each time point is included below each panel. UQ, upper quartile; LQ, lower quartile; Min, minimum.

The most sensitive signaling pathways for each chemical were identified at each time point (Table 4). For each individual chemical, there was no consistency across time points with respect to which pathway was the most sensitive. To further examine this lack of consistency, the pathways at each time point were sorted in descending order based on their median BMD value, and the top 20 most sensitive pathways (i.e., those with the lowest median BMD) were ranked. The number of times a pathway was in the top 20 was totaled, and its rankings were averaged (Table 5). No pathway was consistently in the top 20 at all four time points for any chemical. All chemicals except NDPA and HZBZ had at least one pathway that was in the top 20 at three of the four time points. These results suggest that for most chemicals the most sensitive pathway was not consistent across time points.

TABLE 4

Transcriptional Point-of-Departure Values for the Most Sensitive Signaling Pathway

Chemical No. of genesa Pathway BMD (BMDL) (mg/kg-day or ppm)b Route-adjusted BMD (BMDL) (mg/kg-day)c 
5 Days 
TRBZ 26 (12) Cell cycle_Sister chromatid cohesion 9.2 (6.4) — 
BRBZ 17 (7) HETE and HPETE biosynthesis and metabolism 72.1 (49.0) — 
TTCP 26 (9) Autophagy_Autophagy 59.8 (32.8) — 
MDMB 8 (5) Acetaminophen metabolism 166 (131) 16.6 (13.1) 
NDPA 78 (59) CFTR translational fidelity (class I mutations) 1530 (1153) 173 (130) 
HZBZ 23 (5) CAR-mediated direct regulation of xenobiotic-metabolizing enzymes 69.9 (42.0) 7.0 (4.2) 
2 Weeks 
TRBZ 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 9.1 (6.4) — 
BRBZ 32 (5) Immune response_IL-7 signaling in B lymphocytes 72.6 (49.2) — 
TTCP 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 17.6 (10.8) — 
MDMB 13 (5) Leukotriene 4 biosynthesis and metabolism 226 (141) 22.6 (14.1) 
NDPA 19 (18) Cholesterol biosynthesis 924 (582) 104 (65.7) 
HZBZ 31 (5) Regulation of immune cell differentiation by Notch signaling 77.9 (61.4) 7.8 (6.1) 
4 Weeks 
TRBZ 19 (7) IL-1 beta-dependent CFTR expression 13.5 (8.5) — 
BRBZ 10 (5) Mitochondrial ketone bodies biosynthesis and metabolism 47.6 (27.8) — 
TTCP 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 23.3 (17.1) — 
MDMB 18 (7) Naphthalene metabolism 196 (119) 19.6 (11.9) 
NDPA 52 (9) Nicotine signaling in glutamatergic neurons 1432 (1093) 162 (124) 
HZBZ 17 (10) Immune response_Antigen presentation by MHC class II 62.9 (40.2) 6.3 (4.0) 
13 Weeks 
TRBZ 19 (6) Cholesterol biosynthesis 3.7 (2.4) — 
BRBZ 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 72.0 (57.8) — 
TTCP 15 (6) Methionine metabolism 24.3 (15.0) — 
MDMB 22 (10) Prostaglandin 2 biosynthesis and metabolism 209 (124) 20.9 (12.4) 
NDPA 24 (12) Cell cycle_Role of 14-3-3 proteins in cell cycle regulation 1395 (1077) 158 (122) 
HZBZ 78 (6) CFTR translational fidelity (class I mutations) 36.2 (20.4) 3.6 (2.0) 
Chemical No. of genesa Pathway BMD (BMDL) (mg/kg-day or ppm)b Route-adjusted BMD (BMDL) (mg/kg-day)c 
5 Days 
TRBZ 26 (12) Cell cycle_Sister chromatid cohesion 9.2 (6.4) — 
BRBZ 17 (7) HETE and HPETE biosynthesis and metabolism 72.1 (49.0) — 
TTCP 26 (9) Autophagy_Autophagy 59.8 (32.8) — 
MDMB 8 (5) Acetaminophen metabolism 166 (131) 16.6 (13.1) 
NDPA 78 (59) CFTR translational fidelity (class I mutations) 1530 (1153) 173 (130) 
HZBZ 23 (5) CAR-mediated direct regulation of xenobiotic-metabolizing enzymes 69.9 (42.0) 7.0 (4.2) 
2 Weeks 
TRBZ 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 9.1 (6.4) — 
BRBZ 32 (5) Immune response_IL-7 signaling in B lymphocytes 72.6 (49.2) — 
TTCP 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 17.6 (10.8) — 
MDMB 13 (5) Leukotriene 4 biosynthesis and metabolism 226 (141) 22.6 (14.1) 
NDPA 19 (18) Cholesterol biosynthesis 924 (582) 104 (65.7) 
HZBZ 31 (5) Regulation of immune cell differentiation by Notch signaling 77.9 (61.4) 7.8 (6.1) 
4 Weeks 
TRBZ 19 (7) IL-1 beta-dependent CFTR expression 13.5 (8.5) — 
BRBZ 10 (5) Mitochondrial ketone bodies biosynthesis and metabolism 47.6 (27.8) — 
TTCP 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 23.3 (17.1) — 
MDMB 18 (7) Naphthalene metabolism 196 (119) 19.6 (11.9) 
NDPA 52 (9) Nicotine signaling in glutamatergic neurons 1432 (1093) 162 (124) 
HZBZ 17 (10) Immune response_Antigen presentation by MHC class II 62.9 (40.2) 6.3 (4.0) 
13 Weeks 
TRBZ 19 (6) Cholesterol biosynthesis 3.7 (2.4) — 
BRBZ 9 (5) 2-Naphthylamine and 2-Nitronaphthalene metabolism 72.0 (57.8) — 
TTCP 15 (6) Methionine metabolism 24.3 (15.0) — 
MDMB 22 (10) Prostaglandin 2 biosynthesis and metabolism 209 (124) 20.9 (12.4) 
NDPA 24 (12) Cell cycle_Role of 14-3-3 proteins in cell cycle regulation 1395 (1077) 158 (122) 
HZBZ 78 (6) CFTR translational fidelity (class I mutations) 36.2 (20.4) 3.6 (2.0) 

Note. aTotal number of genes in pathway (number of genes in pathway with BMD and fit p value > 0.1).

bBMD (BMDL).

cRoute-adjusted BMD and BMDL values for MDMB, NDPA, and HZBZ in ppm were converted to milligrams per kilogram-day using strain-specific subchronic EPA food intake factors (EPA, 1988). Male F344 rats, 0.100; Female F344 rats, 0.113.

TABLE 5

Frequency of Signaling Pathways Being in the Top 20 Most Sensitive and Average Rank

Chemical Pathway No. of time points in top 20 Average rank 
TRBZ 2-Naphthylamine and 2-Nitronaphtalene metabolism 
PXR-mediated direct regulation of xenobiotic-metabolizing enzymes 13.3 
CAR-mediated direct regulation of xenobiotic-metabolizing enzymes 3.5 
Butanoate metabolism 7.5 
Propionate metabolism p.2 8.5 
(L)-Alanine, (L)-cysteine, and (L)-methionine metabolism 11.5 
Neurophysiological process_Role of CDK5 in presynaptic signaling 14 
dCTP/dUTP metabolism 15.5 
BRBZ Bile acid biosynthesis 10.3 
CAR-mediated direct regulation of xenobiotic-metabolizing enzymes 11 
Aminoacyl-tRNA biosynthesis in cytoplasm 14.7 
2-Naphthylamine and 2-Nitronaphtalene metabolism 
Influence of low doses of arsenite on glucose stimulated insulin secretion in pancreatic cells 3.5 
Acetaminophen metabolism 4.5 
Mitochondrial ketone bodies biosynthesis and metabolism 
Naphthalene metabolism 
BRBZ metabolism 
Vitamin E (alfa-tocopherol) metabolism 
Mitochondrial unsaturated fatty acid β-oxidation 7.5 
Mitochondrial long-chain fatty acid β-oxidation 
CoA biosynthesis 11 
Estradiol metabolism 12 
Regulation of lipid metabolism_Regulation of fatty acid synthesis: NLTP and EHHADH 12 
Leukotriene 4 biosynthesis and metabolism 15 
Aryl hydrocarbon receptor signaling 17 
TTCP N-Acylethanolamines, HSRL5-transacylation pathway 2.7 
N-Acylethanolamines, Phospholipase A2 pathway 3.7 
2-Naphthylamine and 2-Nitronaphtalene metabolism 
Thiamine metabolism 
1-Naphthylamine and 1-Nitronaphtalene metabolism 4.5 
Vitamin K metabolism 4.5 
BRBZ metabolism 
Selenoamino acid metabolism 
Keratan sulfate metabolism p.1 
Taurine and hypotaurine metabolism 9.5 
Mechanism of Pioglitazone/Glimepiride and Rosiglitazone/Glimepiride cooperative action in diabetes mellitus, Type 2 11 
Acetaminophen metabolism 13 
Glycolysis and gluconeogenesis p. 2 13 
Immune response_IL-12 signaling pathway 13 
MDMB Vitamin K metabolism 3.3 
Benzo[a]pyrene metabolism 4.7 
2-Naphthylamine and 2-Nitronaphtalene metabolism 7.7 
Acetaminophen metabolism 8.3 
1-Naphthylamine and 1-Nitronaphtalene metabolism 
Immune response_IL-15 signaling via JAK-STAT cascade 
Apoptosis and survival_Anti-apoptotic TNFs/NF-kB/IAP pathway 10 
Heme metabolism 12.5 
Prostaglandin 1 biosynthesis and metabolism 14 
Apoptosis and survival_Anti-apoptotic TNFs/NF-kB/Bcl-2 pathway 17.5 
NDPA N-Acylethanolamines. N-Acyltransferase pathway 
Peroxysomal straight-chain fatty acid β-oxidation 
Neurophysiological process_Sweet taste signaling 6.5 
Immune response_Classical complement pathway 10.5 
Glycolysis and gluconeogenesis p. 2 11 
Androstenedione and testosterone biosynthesis and metabolism p.2 17.5 
LRRK2 in neurons in Parkinson’s disease 17.5 
HZBZ Phospholipid metabolism p.3 
Immune response_Antigen presentation by MHC class II 6.5 
Immune response_Antiviral actions of Interferons 13 
LRRK2 and immune function in Parkinson’s disease 13 
Chemical Pathway No. of time points in top 20 Average rank 
TRBZ 2-Naphthylamine and 2-Nitronaphtalene metabolism 
PXR-mediated direct regulation of xenobiotic-metabolizing enzymes 13.3 
CAR-mediated direct regulation of xenobiotic-metabolizing enzymes 3.5 
Butanoate metabolism 7.5 
Propionate metabolism p.2 8.5 
(L)-Alanine, (L)-cysteine, and (L)-methionine metabolism 11.5 
Neurophysiological process_Role of CDK5 in presynaptic signaling 14 
dCTP/dUTP metabolism 15.5 
BRBZ Bile acid biosynthesis 10.3 
CAR-mediated direct regulation of xenobiotic-metabolizing enzymes 11 
Aminoacyl-tRNA biosynthesis in cytoplasm 14.7 
2-Naphthylamine and 2-Nitronaphtalene metabolism 
Influence of low doses of arsenite on glucose stimulated insulin secretion in pancreatic cells 3.5 
Acetaminophen metabolism 4.5 
Mitochondrial ketone bodies biosynthesis and metabolism 
Naphthalene metabolism 
BRBZ metabolism 
Vitamin E (alfa-tocopherol) metabolism 
Mitochondrial unsaturated fatty acid β-oxidation 7.5 
Mitochondrial long-chain fatty acid β-oxidation 
CoA biosynthesis 11 
Estradiol metabolism 12 
Regulation of lipid metabolism_Regulation of fatty acid synthesis: NLTP and EHHADH 12 
Leukotriene 4 biosynthesis and metabolism 15 
Aryl hydrocarbon receptor signaling 17 
TTCP N-Acylethanolamines, HSRL5-transacylation pathway 2.7 
N-Acylethanolamines, Phospholipase A2 pathway 3.7 
2-Naphthylamine and 2-Nitronaphtalene metabolism 
Thiamine metabolism 
1-Naphthylamine and 1-Nitronaphtalene metabolism 4.5 
Vitamin K metabolism 4.5 
BRBZ metabolism 
Selenoamino acid metabolism 
Keratan sulfate metabolism p.1 
Taurine and hypotaurine metabolism 9.5 
Mechanism of Pioglitazone/Glimepiride and Rosiglitazone/Glimepiride cooperative action in diabetes mellitus, Type 2 11 
Acetaminophen metabolism 13 
Glycolysis and gluconeogenesis p. 2 13 
Immune response_IL-12 signaling pathway 13 
MDMB Vitamin K metabolism 3.3 
Benzo[a]pyrene metabolism 4.7 
2-Naphthylamine and 2-Nitronaphtalene metabolism 7.7 
Acetaminophen metabolism 8.3 
1-Naphthylamine and 1-Nitronaphtalene metabolism 
Immune response_IL-15 signaling via JAK-STAT cascade 
Apoptosis and survival_Anti-apoptotic TNFs/NF-kB/IAP pathway 10 
Heme metabolism 12.5 
Prostaglandin 1 biosynthesis and metabolism 14 
Apoptosis and survival_Anti-apoptotic TNFs/NF-kB/Bcl-2 pathway 17.5 
NDPA N-Acylethanolamines. N-Acyltransferase pathway 
Peroxysomal straight-chain fatty acid β-oxidation 
Neurophysiological process_Sweet taste signaling 6.5 
Immune response_Classical complement pathway 10.5 
Glycolysis and gluconeogenesis p. 2 11 
Androstenedione and testosterone biosynthesis and metabolism p.2 17.5 
LRRK2 in neurons in Parkinson’s disease 17.5 
HZBZ Phospholipid metabolism p.3 
Immune response_Antigen presentation by MHC class II 6.5 
Immune response_Antiviral actions of Interferons 13 
LRRK2 and immune function in Parkinson’s disease 13 

Comparison of Transcriptional and Noncancer BMDs for Noncancer Risk Assessment

Given the diversity in target organs affected in this study, multiple modes of action are likely involved in the pathogenesis of the noncancer apical effects across the different chemicals. Therefore, there was no expectation of identifying meaningful pathways that were correlated with the noncancer apical responses across chemicals, as performed in a previous study (Thomas et al., 2012). Instead, the correlation between the noncancer apical BMD values and the most sensitive pathway transcriptional BMD was evaluated to examine the relationship between the initiation of pathological and molecular changes (Fig. 2). The BMD was chosen as the primary metric for comparison because the BMD values represent the central estimate of the dose expected to yield the BMR of interest, whereas the BMDL refers to a statistical lower bound estimate of a confidence interval for the BMD. The BMDL values account for uncertainty in dose-response estimates due to characteristics of experimental design and the endpoint being measured (e.g., sample size and random error). Despite the limited number of chemicals, the overall correlation at each time point between the most sensitive apical BMD value and the most sensitive pathway transcriptional BMD value was greater than 0.9. In only two cases, MDMB at 4 weeks and TTCP at 13 weeks, did the two BMD values differ by more than 10-fold. The consistency of the correlations suggests that the overall relationship between the most sensitive noncancer apical response and the most sensitive transcriptional responses at the pathway level was stable over time. It should be noted that despite the use of six chemicals in the study, only five chemicals had noncancer apical BMD values for comparison with the transcriptional responses (HZBA did not; see Table 2).

FIG. 2.

Scatter plots of the relationship between BMD values for noncancer apical endpoints and transcriptional BMD values for the most sensitive signaling pathway at (A) 5 days, (B) 2 weeks, (C) 4 weeks, and (D) 13 weeks of exposure. The data points were colored based on the target tissue. For each chemical and tissue, the lowest apical BMD value that was considered adverse was plotted. All BMD values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equal apical and transcriptional BMD values. The red lines represent 10-fold difference between the BMD values. Correlation coefficients (r) are included in the lower right hand corner of each panel.

FIG. 2.

Scatter plots of the relationship between BMD values for noncancer apical endpoints and transcriptional BMD values for the most sensitive signaling pathway at (A) 5 days, (B) 2 weeks, (C) 4 weeks, and (D) 13 weeks of exposure. The data points were colored based on the target tissue. For each chemical and tissue, the lowest apical BMD value that was considered adverse was plotted. All BMD values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equal apical and transcriptional BMD values. The red lines represent 10-fold difference between the BMD values. Correlation coefficients (r) are included in the lower right hand corner of each panel.

In order to increase the statistical power to detect a relationship between a given noncancer apical response and the most sensitive pathway transcriptional response, the data from the five chemicals in this study at the 13-week time point were combined with data from a previous study (Thomas et al., 2012). The combined data consisted of 10 chemicals, four target organs (bladder, liver, thyroid, and lung), multiple routes of exposure (gavage, feed, and inhalation), and both sexes of two species (mice and rats). The overall correlation between the noncancer apical BMD values and the most sensitive pathway transcriptional BMD values was 0.827 (p = 0.0031), with a median apical-to-transcriptional BMD ratio of 1.05 (Fig. 3). A lower correlation was observed for the associated BMDL values.

FIG. 3.

Scatter plots of the relationship between (A) BMD and (B) BMDL values for noncancer apical endpoints and transcriptional BMD and BMDL values for the most sensitive signaling pathway at 13 weeks of exposure when combined with data from a previous study (Thomas et al., 2012). The data points were colored based on the target tissue, and symbol shape represents species. For each chemical and tissue, the lowest apical BMD and BMDL value that was considered adverse was plotted. All BMD and BMDL values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equal apical and transcriptional BMD and BMDL values. The red lines represent 10-fold difference between the BMD and BMDL values. Correlation coefficients (r) and associated p value are included in the lower right hand corner of each panel. The median log2 ratio of the apical-to-transcriptional BMD values is provided in the upper left hand corner. The antilog of the transformed ratio is provided in parentheses.

FIG. 3.

Scatter plots of the relationship between (A) BMD and (B) BMDL values for noncancer apical endpoints and transcriptional BMD and BMDL values for the most sensitive signaling pathway at 13 weeks of exposure when combined with data from a previous study (Thomas et al., 2012). The data points were colored based on the target tissue, and symbol shape represents species. For each chemical and tissue, the lowest apical BMD and BMDL value that was considered adverse was plotted. All BMD and BMDL values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equal apical and transcriptional BMD and BMDL values. The red lines represent 10-fold difference between the BMD and BMDL values. Correlation coefficients (r) and associated p value are included in the lower right hand corner of each panel. The median log2 ratio of the apical-to-transcriptional BMD values is provided in the upper left hand corner. The antilog of the transformed ratio is provided in parentheses.

As noted previously (Thomas et al., 2011), a direct analysis of the relationship between noncancer apical BMD values and pathway transcriptional BMD values is confounded by the correlation between the noncancer apical BMD values and the MTD. In order to understand the extent of the confounding by the MTD, MTD values were collected for 8 of the 10 chemicals. The two chemicals without MTDs, TRBZ and propylene glycol mono-t-butyl ether, did not show a 10% decrease in body weight or the development of a limiting pathological lesion at any of the doses tested in our study or previous studies (Carlson and Tardiff, 1977; NTP, 2004). One chemical in our study, HZBZ, did not have a noncancer apical BMD value. Based on eight chemicals, the correlation between the noncancer apical BMD values and the transcriptional BMD values for the most sensitive pathway was 0.9946. By comparison, the correlation between the noncancer apical BMD values and the MTD was 0.9937. The partial correlation coefficient between the noncancer apical BMD values and the pathway transcriptional BMD values after accounting for the MTD was 0.3912 (p = 0.338). This result suggests that the overall correlation between the noncancer apical BMD values and the pathway transcriptional BMD values was equivalent to that observed with the MTD.

Comparison of Transcriptional and Tumor Incidence BMDs for Cancer Risk Assessment

Similar to the noncancer endpoints, multiple modes of action were likely involved in the resultant cancer-related apical effects across the three different chemicals that targeted three different organ systems. Therefore, no attempt was made to identify common pathways that were correlated with the cancer-related apical responses across chemicals. To examine the overall relationship between cancer-related pathological effects and the initiation of molecular changes, the cancer-related apical BMD values for each chemical were plotted as a function of the transcriptional BMD values for the most sensitive pathway (Fig. 4). Similar to the analysis of the noncancer apical endpoints, the BMD, as opposed to the BMDL, was chosen as the basis of comparison because the BMD values represent the central estimate of the dose expected to yield the BMR of interest. Because only three of the six chemicals had existing rodent cancer bioassays, no correlation coefficients were calculated. For HZBZ, no adverse noncancer pathological changes were present at any time point (Table 2); nonetheless, the transcriptional BMD values for the most sensitive pathway were similar to the tumor BMD value.

FIG. 4.

Scatter plot of the relationship between BMD values for cancer-related apical endpoints and transcriptional BMD values for the most sensitive signaling pathway. The data points were colored based on exposure duration. All BMD values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equivalent apical and transcriptional BMD values. The red lines represent 10-fold difference between the BMD values.

FIG. 4.

Scatter plot of the relationship between BMD values for cancer-related apical endpoints and transcriptional BMD values for the most sensitive signaling pathway. The data points were colored based on exposure duration. All BMD values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equivalent apical and transcriptional BMD values. The red lines represent 10-fold difference between the BMD values.

In order to increase the statistical power to detect a relationship between the cancer-related apical response and the most sensitive pathway transcriptional response, the data from the three chemicals in this study at the 13-week time point were combined with data from a previous study (Thomas et al., 2012). The combined data had eight chemicals with one chemical that produced tumors in two different tissues. The combined data set also had four target organs (bladder, liver, thyroid, and lung), multiple routes of exposure (gavage, feed, and inhalation), and both sexes of two species (mice and rats). The overall correlation between the cancer-related apical BMD values and the transcriptional BMD values for the most sensitive pathway was 0.940 (p = 0.0002), with a median apical-to-transcriptional BMD ratio of 1.44 (Fig. 5). A slightly better correlation was observed between the corresponding BMDL values.

FIG. 5.

Scatter plots of the relationship between (A) BMD and (B) BMDL values for cancer-related apical endpoints and transcriptional BMD and BMDL values for the most sensitive signaling pathway at 13 weeks of exposure when combined with data from a previous study (Thomas et al., 2012). The data points were colored based on the target tissue, and symbol shape represents species. All BMD and BMDL values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equal apical and transcriptional BMD and BMDL values. The red lines represent 10-fold difference between the BMD and BMDL values. Correlation coefficients (r) and associated p value are included in the lower right hand corner of each panel. The median log2 ratio of the apical-to-transcriptional BMD values is provided in the upper left hand corner. The antilog of the transformed ratio is provided in parentheses.

FIG. 5.

Scatter plots of the relationship between (A) BMD and (B) BMDL values for cancer-related apical endpoints and transcriptional BMD and BMDL values for the most sensitive signaling pathway at 13 weeks of exposure when combined with data from a previous study (Thomas et al., 2012). The data points were colored based on the target tissue, and symbol shape represents species. All BMD and BMDL values were converted to milligrams per kilogram-day to be on a common scale. The black line dissecting the graph indicates equal apical and transcriptional BMD and BMDL values. The red lines represent 10-fold difference between the BMD and BMDL values. Correlation coefficients (r) and associated p value are included in the lower right hand corner of each panel. The median log2 ratio of the apical-to-transcriptional BMD values is provided in the upper left hand corner. The antilog of the transformed ratio is provided in parentheses.

As noted for the noncancer effects, a direct analysis of the relationship between the cancer-related apical BMD values and the pathway transcriptional BMD values is confounded by the correlation between the tumor BMD values and the MTD (Krewski et al., 1993; Razzaghi and Gaylor, 1996; Thomas et al., 2011). MTD values were available for seven of the eight chemicals with positive rodent cancer bioassays, with propylene glycol mono-t-butyl ether being the only chemical without an MTD. Based on these seven chemicals, with one chemical, methylene chloride, having tumor responses in two different tissues, the correlation between cancer-related apical BMD values and the transcriptional BMD values for the most sensitive pathway was 0.978. By comparison, the correlation between the tumor BMD values and the MTD was 0.784. The partial correlation coefficient between cancer-related apical BMD values and the pathway transcriptional BMD values after accounting for the MTD was 0.990 (p < 0.0001). This analysis indicates that the overall correlation between the cancer-related apical BMD values and the transcriptional BMD values for the most sensitive pathway is greater than that observed with the MTD.

DISCUSSION

Despite the rapid technological advances over the past decade, incorporation of newer technologies such as trancriptomics to chemical risk assessment has been relatively slow. In some respects, the delay in applying transcriptomics to risk assessment has been due to challenges associated with linking transcriptional perturbation with adverse effects and using transcriptomic data to define the mode of action for a chemical with a level of confidence required for regulatory decision making. To compensate for this lack of confidence, proposals for incorporating transcriptomic data into risk assessment have focused on using it as part of an overall weight of evidence for a mode of action (EPA, 2009; NRC, 2007). Although the challenges associated with interpreting transcriptomic data may be overcome in the long term, leveraging other strengths of the technology, such as the comprehensive characterization of transcriptional perturbations, may be utilized in the near term to make more efficient and economical decisions for the multitude of chemicals that lack published reference values.

In this study, the application of transcriptomic data to chemical risk assessment was evaluated in two ways. First, we examined the relationship between the dose that causes pathological changes and the dose that begins to perturb molecular pathway. In other words, we did not attempt to filter the transcriptomic results based on a relationship of the gene families to pathology or mode of action. We simply assessed whether a point of departure based on transcriptional perturbations would be protective of the observed apical responses. Second, we evaluated the manner in which duration of exposure affected the relationship between apical and transcriptional BMD values. The results demonstrate that transcriptional BMD values based on the most sensitive pathway are highly correlated with BMD values based on traditional apical responses for both noncancer and cancer-related endpoints. The correlation was robust across species, sex, tissues, and routes of exposure. Further, transcriptional perturbation did not occur at doses that were orders of magnitude lower than those producing pathological changes or tumors: the transcriptional BMD values based on the most sensitive pathway were, on average, less than twofold different than the apical BMD values.

The rationale for using the transcriptional BMD values from the most sensitive pathway was twofold. First, the transcriptional BMD value from the most sensitive pathway was generally conservative relative to the apical response. For example, 60% of the chemicals had transcriptional BMD values that were at or below the noncancer apical BMD values at the 13-week time point, whereas at the 5-day time point, all chemicals had transcriptional BMD values that were at or below the noncancer apical BMD values. For cancer-related responses, 89% of the chemicals had transcriptional BMD values that were at or below the apical BMD values. Second, the use of transcriptional BMD values from the most sensitive pathway did not require a priori knowledge of the linkage between specific pathways and adverse apical endpoints nor did it require additional supporting mode-of-action studies for use in a weight-of-evidence framework. Instead, transcriptomic testing of new chemicals can be performed regardless of prior knowledge of mode of action or biological assumptions.

The correlation between the pathway transcriptional BMD values and the apical BMD values for both noncancer and cancer-related endpoints is confounded by the correlation with the MTD (Thomas et al., 2011). This confounding occurs for both technical and biological reasons. For example, animals cannot be given a dose higher than the MTD for a specific chemical without severe weight loss, limiting pathological lesions, or death. This upper bound constrains the range of the dose-response curve and the possible BMD values. In addition, cytotoxicity frequently occurs at doses near the MTD, which, by itself, can increase the risk of developing tumors. To account for this confounding, the correlation between the transcriptional BMD values and the apical BMD values was compared with the correlation between the apical BMD values and the MTD. Further, a partial correlation was calculated between the transcriptional BMD values and the apical BMD values when adjusted for the MTD.

For the noncancer endpoints, the correlation between the transcriptional BMD values for the most sensitive pathway and the apical BMD values was equivalent with that between the apical BMD values and the MTD. However, both correlation coefficients were extremely high (0.9946 vs. 0.9937). It is not possible to discriminate between the two because of the limited number of chemicals that had MTD values. More chemicals will ultimately need to be evaluated in order to determine the utility of the pathway transcriptional BMD values over the MTD for noncancer endpoints. Advantages of the transcriptomic endpoints are that they are mechanistically anchored, which may provide a bridge to other high information content data streams (e.g., metabolomics, high-throughput in vitro screening) and may ultimately lead to more biologically based methods of predicting risk. For cancer-related endpoints, the correlation between the transcriptional BMD values for the most sensitive pathway and the apical BMD values was larger than that between the apical BMD values and the MTD. The partial correlation coefficient was statistically significant. In this case, the transcriptional BMD values based on the most sensitive pathway appear to provide significant value in estimating the apical BMD values relative to the MTD.

Although the largest data set for all chemicals was available at the 13-week time point, a subset of chemicals was analyzed at multiple exposure durations. The results suggest that the correlation between the transcriptional BMD values for the most sensitive pathway and the apical BMD values was relatively stable over time for both noncancer- and cancer-related endpoints. As a result, full 13-week exposures may not be necessary to estimate apical BMD values for an unknown chemical and that shorter term, more cost-effective studies may suffice. However, even though the pathway transcriptional and apical BMD values were generally correlated at each time point, certain chemicals showed a progressive reduction in BMD values. For example, TRBZ showed a 3- and 2.4-fold reduction in the apical and pathway transcriptional BMD values, respectively, from 5 days to 13 weeks. The largest difference among the pathway transcriptional BMD values between the 5-day and 13-week time point was 2.5-fold. Although longer duration exposures may be preferred to identify these chemicals, it may not be worth the increased costs associated with performing a 13-week study. One potential solution may be the inclusion of an additional uncertainty factor of two- to threefold when estimating the point of departure using a 5-day or 2-week transcriptomic study.

For most of the chemicals analyzed in this study, the most sensitive pathway was not consistent across time points. However, just as histopathological changes evolve over time in response to treatment to reflect the temporal evolution of tissue injury, remodeling, and adaptation, so will gene expression changes. For example, a large number of the most sensitive pathways for TRBZ at the 5-day time point were involved in cell cycle regulation. This is consistent with the histopathological data that reported a minimal increase in mitotic nuclei only at this time point (Dodd et al., 2012c). At the 2-week time point, a variety of metabolic pathways (e.g., amino acid, short-chain fatty acids, nucleotide, and xenobiotic) were among the most sensitive for TRBZ. Previous studies have reported treatment-related increases in the activities of liver cytochrome P450, NADPH cytochrome c reductase, O-ethyl O-p-nitrophenyl phenylphosphonothioate detoxification, azoreductase, and benzpyrene hydroxylase (Carlson and Tardiff, 1977). At 4 weeks, some of the metabolic pathways persisted, but immune response and apoptosis-related pathways were present. Finally, at 13 weeks, the top two most sensitive pathways were involved in lipid metabolism, which is consistent with the histological appearance of hepatocyte vacuolation at this time point. Even though they were not always in the top 20 most sensitive pathways, the BMD values for the CAR and PXR pathways were consistently between 10 and 20mg/kg-day across all four time points. Activation of CAR and PXR at these doses is consistent with the increased liver weight and hepatocyte hypertrophy observed at doses ≥ 10mg/kg-day (Dodd et al., 2012c).

The combined results from the study can be used to refine a previously proposed framework for applying transcriptomic data to risk assessment (Thomas et al., 2011). In a modified framework, dose-response studies could be performed on untested chemicals at any single time point between 5 days and 13 weeks (Fig. 6). The tests would be performed in mice and rats of both sexes. A battery of eight tissues that include those most frequently positive in rodent cancer bioassays (liver, lung, mammary gland, stomach, vascular system, kidney, hematopoietic system, and urinary bladder) could then be collected. In a previous analysis, these eight tissues cover 92 and 82% of targets for all mouse and rat carcinogens, respectively (Gold et al., 2001). Additional tissues that could include brain, the developing fetus, and gonadal tissue could be added for assessing neurological, developmental, and reproductive effects. Gene expression microarray analysis on these tissues would allow the estimation of pathway transcriptional BMD and BMDL values. The lowest transcriptional BMD value across all analyzed tissues would be used to estimate a point of departure. In parallel, assays for genotoxicity could further expand the utility of the transcriptional BMD values. A weight-of-evidence analysis could then be used to determine the genotoxicity potential of the chemical and inform decisions regarding the use of the transcriptional BMD values for cancer-related slope factors or noncancer reference values.

FIG. 6.

A flow chart outlining the application of transcriptomic data to chemical risk assessment. The dark gray boxes represent the experimental studies and weight-of-evidence evaluation required to support the assessments. The blue and red boxes represent the steps involved for noncancer- and cancer-related assessments, respectively.

FIG. 6.

A flow chart outlining the application of transcriptomic data to chemical risk assessment. The dark gray boxes represent the experimental studies and weight-of-evidence evaluation required to support the assessments. The blue and red boxes represent the steps involved for noncancer- and cancer-related assessments, respectively.

The application of transcriptomics within this framework poses both significant benefits and challenges. Importantly, the framework provides a relatively rapid and cost-effective means to provide provisional points of departure for chemicals without published reference values. The provisional transcriptional points of departure could be easily incorporated into the second tier of EPA’s NexGen risk assessment program (Cote et al., 2012). The second tier was designed to use high- and medium-throughput bioassay data for decision making on chemicals with limited traditional data. Similar to the proposal by Cote and colleagues, if or when more data are collected, the provisional transcriptional points of departure could either be updated or transitioned to a higher risk assessment tier. Among the obvious challenges, the proposed framework would require a shift from the current hazard-based labeling paradigm that relies on apical responses to a “region of safety” paradigm where the most sensitive adverse effect is not known, but the regulatory dose is set based on the lack of transcriptional perturbation. Concerns that may accompany this shift in focus may be whether the pathway transcriptional BMD values are adequately protective across a wide range of chemicals and the difficulty in associating a transcriptional perturbation with actual risk. Although both concerns are valid, it is generally accepted that adverse phenotypic responses would not occur without direct or indirect alterations in transcriptional programs (Farr and Dunn, 1999). In addition, for many chemicals, a point of departure based on transcriptional perturbation is preferable to having no point of departure at all.

Additional concerns may include whether a point of departure based on transcriptional perturbation is overly protective, especially when a given perturbation may not have a direct link to an adverse effect, and whether “no transcriptional perturbation” could become the new default. Based on the chemicals we have studied, concern that transcriptional perturbation would be overly protective may be unfounded because the transcriptional BMD values for the most sensitive pathway were, on average, within a factor of two of both the noncancer- and cancer-related apical BMD values. However, this cannot be ruled out for a broader set of chemicals and mechanisms. The concern that this may potentially establish a new default is valid. However, traditional animal studies could be performed if it is believed that the “no transcriptional perturbation” level was too conservative.

The proposed framework outlines a practical, near-term application of transcriptomic data to chemical risk assessment. Over the longer term, better methods for interpreting transcriptomic data in the context of mode of action and adverse endpoints will become available and should be utilized. In addition, the transcriptomic data may be combined with other emerging technologies such as high-throughput in vitro screening to identify molecular initiating events and assemble the likely sequence of key events in the mode of action. The combined use of these technologies for identifying mode of action will ultimately lead to more rapid and economical risk assessments while also helping to extrapolate across species and understand human relevance.

SUPPLEMENTARY DATA

Supplementary data are available online at http://toxsci.oxfordjournals.org/.

FUNDING

American Chemistry Council ’s Long-Range Research Initiative.

ACKNOWLEDGMENTS

This manuscript has been reviewed in accordance with the policy of the National Center for Environmental Assessment, U.S. EPA, and approved for publication. The views expressed herein are those of the author(s) and do not necessarily reflect the views or policies of the U.S. EPA. The mention of trade names or commercial products does not constitute endorsement or recommendations for use.

REFERENCES

Akaike
H
.
1973
.
Information theory and an extension of the maximum likelihood principle
. In
2nd International Symposium on Information Theory
 , (B.N. Petrov and F. Csaki, Eds.), Budapest; Academiai kiado:
267
281
.
Carlson
G. P.
Tardiff
R. G
.
1977
.
Effect of 1,4-dibromobenzene and 1,2,4-tribromobenzene on xenobiotic metabolism
.
Toxicol. Appl. Pharmacol
 .
42
,
189
196
.
Cote
I.
Anastas
P. T.
Birnbaum
L. S.
Clark
R. M.
Dix
D. J.
Edwards
S. W.
Preuss
P. W
.
2012
.
Advancing the next generation of health risk assessment
.
Environ. Health Perspect
 .
120
,
1499
1502
.
Dodd
D. E.
Pluta
L. J.
Sochaski
M. A.
Banas
D. A.
Thomas
R. S
.
2012
.
Subchronic hepatotoxicity evaluation of 2,3,4,6-tetrachlorophenol in sprague dawley rats
.
J. Toxicol
 .
2012
,
376246
.
Dodd
D. E.
Pluta
L. J.
Sochaski
M. A.
Banas
D. A.
Thomas
R. S
.
2012
.
Subchronic hepatotoxicity evaluation of bromobenzene in Fischer 344 rats
.
J. Appl. Toxicol
 .
33
,
370
377
.
Dodd
D. E.
Pluta
L. J.
Sochaski
M. A.
Funk
K. A.
Thomas
R. S
.
2012
.
Subchronic hepatotoxicity evaluation of 1,2,4-tribromobenzene in Sprague-Dawley rats
.
Int. J. Toxicol
 .
31
,
250
256
.
Dodd
D. E.
Pluta
L. J.
Sochaski
M. A.
Funk
K. A.
Thomas
R. S
.
2012
.
Subchronic thyroid toxicity evaluation of 4,4’-methylenebis(N,N’-dimethyl)aniline in Fischer 344 rats
.
J. Toxicol. Environ. Health. A
 .
75
,
637
648
.
Dodd
D. E.
Pluta
L. J.
Sochaski
M. A.
Funk
K. A.
Thomas
R. S
.
2012
.
Subchronic urinary bladder toxicity evaluation of N-Nitrosodiphenylamine in Fischer 344 rats
.
J. Appl. Toxicol
 .
33
,
383
389
.
Dodd
D. E.
Pluta
L. J.
Sochaski
M. A.
Wall
H. G.
Thomas
R. S
.
2012
.
Subchronic hepatotoxicity evaluation of hydrazobenzene in Fischer 344 rats
.
Int. J. Toxicol
 .
31
,
564
571
.
EPA
1988
.
Recommendations for and Documentation of Biological Values for Use in Risk Assessment
 , 600/6–87/008.
U.S. Environmental Protection Agency
,
Washington, DC
.
EPA
1989
.
Risk Assessment Guidance for Superfund Volume I: Human Health Evaluation Manual (Part A)
 , EPA/540/1–89/002.
U.S. Environmental Protection Agency
,
Washington, DC
.
EPA
1997
.
Chemistry Assistance Manual for Premanufacture Notification Submitters
 , EPA 744-R-97-003.
U.S. Environmental Protection Agency
,
Washington, DC
.
EPA
2002
.
A Review of the Reference Dose and Reference Concentration Processes
 , EPA/630/P-02/002F.
U.S. Environmental Protection Agency
Washington, DC
.
EPA
2005
.
Guidelines for Carcinogen Risk Assessment
 , EPA/630/P-03/001F.
U.S. Environmental Protection Agency
Washington, DC
.
EPA
2009
.
An Approach to Using Toxicogenomic Data in U.S. EPA Human Health Risk Assessments: A Dibutyl Phthalate Case Study
 , EPA/600/R-09/028F.
National Center for Environmental Assessment, U.S. Environmental Protection Agency
,
Washington, DC
.
Farr
S.
Dunn
R. T.
2nd
1999
.
Concise review: Gene expression applied to toxicology
.
Toxicol. Sci
 .
50
,
1
9
.
Gold
L. S.
Manley
N. B.
Slone
T. H.
Ward
J. M
.
2001
.
Compendium of chemical carcinogens by target organ: Results of chronic bioassays in rats, mice, hamsters, dogs, and monkeys
.
Toxicol. Pathol
 .
29
,
639
652
.
Irizarry
R. A.
Bolstad
B. M.
Collin
F.
Cope
L. M.
Hobbs
B.
Speed
T. P
.
2003
.
Summaries of Affymetrix GeneChip probe level data
.
Nucleic Acids Res
 .
31
,
e15
.
Krewski
D.
Gaylor
D. W.
Soms
A. P.
Szyszkowicz
M
.
1993
.
An overview of the report: Correlation between carcinogenic potency and the maximum tolerated dose: Implications for risk assessment
.
Risk Anal
 .
13
,
383
398
.
NRC
2007
.
Applications of Toxicogenomic Technologies to Predictive Toxicology and Risk Assessment
 .
National Academy of Sciences Press
,
Washington, DC
.
NTP
1978
.
Bioassay of Hydrazobenzene for Possible Carcinogenicity
 , 92.
U.S. Department of Health and Human Services National Toxicology Program
,
Washington, DC
.
NTP
1979
.
Bioassay of 4,4’-Methylenebis-(N,N-Dimethyl)Benzenamine for Possible Carcinogenicity
 , 186.
U.S. Department of Health and Human Services National Toxicology Program
,
Washington, DC
.
NTP
1979
.
Bioassay of N-Nitrosodiphenylamine for Possible Carcinogenicity
 , 164.
U.S. Department of Health and Human Services National Toxicology Program
,
Washington, DC
.
NTP
2004
.
Toxicology and Carcinogenesis Studies of Propylene Glycol Mono-t-Butyl Ether in F344/N Rats and B6C3F1 Mice
 , 515.
U.S. Department of Health and Human Services National Toxicology Program
Washington, DC
.
Posada
D.
Buckley
T. R
.
2004
.
Model selection and model averaging in phylogenetics: Advantages of akaike information criterion and bayesian approaches over likelihood ratio tests
.
Syst. Biol
 .
53
,
793
808
.
Razzaghi
M.
Gaylor
D. W
.
1996
.
On the correlation coefficient between the TD50 and the MTD
.
Risk Anal
 .
16
,
107
113
.
Steel
R. G.
Torrie
J. H
.
1980
.
Principles and Procedures of Statistics: A Biometrical Approach
 .
McGraw Hill
,
New York, NY
.
Thomas
R. S.
Allen
B. C.
Nong
A.
Yang
L.
Bermudez
E.
Clewell
H. J.
3rd
Andersen
M. E
.
2007
.
A method to integrate benchmark dose estimates with genomic data to assess the functional effects of chemical exposure
.
Toxicol. Sci
 .
98
,
240
248
.
Thomas
R. S.
Clewell
H. J.
3rd
Allen
B. C.
Wesselkamper
S. C.
Wang
N. C.
Lambert
J. C.
Hess-Wilson
J. K.
Zhao
Q. J.
Andersen
M. E
.
2011
.
Application of transcriptional benchmark dose values in quantitative cancer and noncancer risk assessment
.
Toxicol. Sci
 .
120
,
194
205
.
Thomas
R. S.
Clewell
H. J.
3rd
Allen
B. C.
Yang
L.
Healy
E.
Andersen
M. E
.
2012
.
Integrating pathway-based transcriptomic data into quantitative chemical risk assessment: A five chemical case study
.
Mutat. Res
 .
746
,
135
143
.
Yang
L.
Allen
B. C.
Thomas
R. S
.
2007
.
BMDExpress: A software tool for the benchmark dose analyses of genomic data
.
BMC Genomics
 .
8
,
387
.