Gene expression profiling to identify potentially relevant disease outcomes and support human health risk assessment for carbon black nanoparticle exposure.

New approaches are urgently needed to evaluate potential hazards posed by exposure to nanomaterials. Gene expression profiling provides information on potential modes of action and human relevance, and tools have recently become available for pathway-based quantitative risk assessment. The objective of this study was to use toxicogenomics in the context of human health risk assessment. We explore the utility of toxicogenomics in risk assessment, using published gene expression data from C57BL/6 mice exposed to 18, 54 and 162 μg Printex 90 carbon black nanoparticles (CBNP). Analysis of CBNP-perturbed pathways, networks and transcription factors revealed concomitant changes in predicted phenotypes (e.g., pulmonary inflammation and genotoxicity), that correlated with dose and time. Benchmark doses (BMDs) for apical endpoints were comparable to minimum BMDs for relevant pathway-specific expression changes. Comparison to inflammatory lung disease models (i.e., allergic airway inflammation, bacterial infection and tissue injury and fibrosis) and human disease profiles revealed that induced gene expression changes in Printex 90 exposed mice were similar to those typical for pulmonary injury and fibrosis. Very similar fibrotic pathways were perturbed in CBNP-exposed mice and human fibrosis disease models. Our synthesis demonstrates how toxicogenomic profiles may be used in human health risk assessment of nanoparticles and constitutes an important step forward in the ultimate recognition of toxicogenomic endpoints in human health risk. As our knowledge of molecular pathways, dose-response characteristics and relevance to human disease continues to grow, we anticipate that toxicogenomics will become increasingly useful in assessing chemical toxicities and in human health risk assessment. Crown

Introduction

Chronic inhalation of fine and ultrafine particulate matter has been associated with adverse pulmonary effects including fibrosis and cancer, as well as exacerbation of existing conditions such as asthma, bronchitis and chronic obstructive pulmonary disorder (Bonner, 2007, Knaapen et al., 2004), in addition to cardiovascular disease (Dockery et al., 1993, Pope et al., 2004). Human exposure to manufactured nanomaterials (NMs), which have at least one size dimension that is less than 100 nm, may constitute an increased risk of adverse effects especially following inhalation exposure, and their potential to induce toxic effects is poorly understood (Handy and Shaw, 2007). Moreover, the human health risks associated with inhalation exposure have not been adequately investigated. Methods that can be effective in screening for NM toxicities are paramount, due to the countless variations in physical and chemical properties of NMs in terms of size, shape, agglomeration and surface coatings.

Traditional assays used in human health risk assessment (HHRA) generally involve chronic and subchronic rodent exposures with concomitant analyses of tumour induction (e.g., two-year rodent cancer bioassay), in addition to various non-cancer endpoints, the most sensitive of which is used for regulatory decision-making (Meek et al., 1994). These approaches form the foundation of the chemical regulatory system and have been invaluable for HHRA. However, some of these assays, such as those based on chronic animal exposures at the maximum tolerated dose, are time and resource intensive, thus limiting broad application (Suter et al., 2004). Recent discussions have identified gene expression profiling as a potentially rapid and cost-effective approach for identifying and assessing prospective hazard, characterizing chemical (or particle) mode of action, and assessing human relevance in support of HHRA (National Academy of Sciences, 2007). In order for gene expression data to become accepted for routine use in HHRA, it is necessary to demonstrate that mRNA/protein expression profiles can effectively predict the modes of action and biological outcomes of exposure at relevant doses, and to confirm that these data can be used to strengthen the foundation for HHRA and regulatory decisions. In this regard, it has been hypothesized that gene expression profiling will be extremely useful in identifying effects at low doses, and moreover, useful for distinguishing between doses that elicit an adaptive response vs. those that yield adverse effects (Boverhof and Zacharewski, 2006). To date, the application of gene expression profiling in regulatory toxicology has largely focused on qualitative identification of chemical modes of action and transcription biomarkers that can predict specific toxicities. However, the utility of gene expression profiling in quantitative determination of threshold values (e.g., benchmark doses) has not yet been rigorously explored (Thomas et al., 2012).

In the present study we investigate the utility of gene expression profiles derived from mice exposed to Printex 90 carbon black nanoparticles (CBNPs) by intratracheal installation to identify potential hazards, modes of action, and doses above which adverse effects may be expected for specific toxicological outcomes. In addition, we quantitatively compare benchmark doses for pathways to those of apical endpoints derived from the same experimental animals. We employ Printex 90 as a model NM due to the rich database of traditional toxicity information on which our findings can be anchored. Briefly, Printex 90 consists almost entirely of carbon, with very low levels of impurities in terms of polycyclic aromatic hydrocarbons and endotoxins (Bourdon et al., 2012b, Jacobsen et al., 2008, Saber et al., 2011) They generate reactive oxygen species (Jacobsen et al., 2008), induce DNA strand breaks in vitro and in vivo (Jacobsen et al., 2009, Saber et al., 2005) and mutations in vitro (Jacobsen et al., 2007) that are associated with oxidative stress (Jacobsen et al., 2011). The data in this study are from previously published experiments investigating Printex 90 CBNP exposure in C57BL/6 mice at various doses (i.e., vehicle, 18, 54 and 162 μg) collected at several time-points (1, 3 and 28 days) following a single acute instillation (Bourdon et al., 2012a). We previously characterized widespread changes in gene expression involving acute phase response and inflammation, supported by concomitant influxes of pulmonary bronchoalveolar lavage cells (BAL) and increases in tissue-specific DNA strand breaks (Bourdon et al., 2012a, Bourdon et al., 2012b). In addition to the examination of BMDs and BMDLs, we compare CBNP-modified gene expression profiles to various models of lung disease in mice and humans reported in the literature, in order to explore the utility of our data in predicting the potential risk of adverse health outcomes and the human relevance of expression changes. The work demonstrates one approach by which gene expression profiling may be integrated into HHRA to support or predict apical toxicological endpoints, dose–response, and relevance to human diseases.

Materials and methods

Summary of experimental model and published results

Details of the mouse exposures, particle characterization and pulmonary phenotype were previously published in Bourdon et al., 2012a, Bourdon et al., 2012b. Briefly, female C57BL/6 mice were exposed to a single installation of vehicle or Printex 90 (18, 54 or 162 μg) and euthanized 1, 3 and 28 days post-exposure (n  = 6/group). The intratracheal instillation route of exposure allows for deposition of known doses directly in the lungs of the mice, and controls for potential dermal- and ingestion-related CBNP exposure that can occur during whole body inhalation exposures. The doses were selected to represent 1, 3 and 9 working days of exposure at the occupational inhalation exposure limit of 3.5 mg/m3 of CB (as established by the US Occupational Safety and Health Administration (OSHA) and the US National Institute for Occupational Safety and Health (NIOSH))) for a mouse (assuming 1.8 L/h inhalation rate and 33.8% particle deposition in mouse, for an 8 h working day) (Dybing et al., 1997, Jacobsen et al., 2009). Very limited filtration of CBNPs from the nose is expected during human exposure. Printex 90 CBNPs were characterized and displayed the following properties: 14 nm primary particle size, 295–338 m2/g Brunauer Emmett and Teller (BET) surface area, 74.2 μg/g PAHs, 142 EU/g endotoxin, polydispersity index of 1, −10.7 mV zeta potential, 2.6 μm peak hydrodynamic number and 3.1 μm peak volume-size-distribution (Bourdon et al., 2012b).

Analysis of pulmonary inflammatory cellular influx in bronchoalveolar lavage (BAL) revealed neutrophilic inflammation that was sustained to day 28 at all doses. Tissue-specific genotoxicity, as observed by DNA strand breaks, persisted up to day 28 at the two highest doses and FPG-sensitive sites at all doses on day 1 and the highest dose on day 3 (Bourdon et al., 2012b). Whole mouse genome DNA microarray revealed 487 and 81 differentially expressed genes (FDR adjusted p-value ≤ 0.1 and fold changes ≥ 1.5) overall in lung and liver, respectively (Bourdon et al., 2012a). The complete microarray dataset is available through the Gene Expression Omnibus at NCBI (http://www.ncbi.nlm.nih.gov/geo/, Superseries GSE35284, SubSeries GSE35193). This dataset was previously used to examine molecular interactions between lung and liver upon CBNP exposure (Bourdon et al., 2012a).

Rank-based pathway analysis

To determine the most affected processes of CBNP exposure, pathway analysis of gene expression data was conducted using a rank based test in R (R Development Core Team, 2011) as described in Alvo et al. (2010). The relative expression for the genes in a pathway was first aligned by subtracting the median expression value for the combined treatment and control groups. These values were then ranked within each subject and the vector of average ranks was calculated for each treatment group. The distance between the two treatments was calculated and a permutation analysis was used to obtain a p-value for each pathway. Pathways with p  < 0.05 were considered significant.

Benchmark dose (BMD)/lower confidence limit benchmark dose (BMDL) calculation for apical endpoints and RT-PCR data

BMD10 (BMD representing an excess risk of 10% in exposed animals vs. controls) and BMDLs (95% confidence limit) were calculated for apical endpoint data (inflammation and genotoxicity) and for RT-PCR using EPA BMDS 2.2 (Davis et al., 2011). Only data that were statistically above control levels (p  < 0.05) for at least two of the doses were included. Prior to running the analysis, the data were screened for homogeneity of variance, and then fit against five continuous dose–response models (i.e., hill, polynomial, linear, power and exponential). Goodness of fit >0.05 and scaled residuals within ±2.0 was applied as a cut off for selection of the appropriate model, and curves were also inspected visually. When more than one model was suitable, the one with the lowest Akaike's information criterion (AIC) was selected.

BMD/BMDL calculation for genomic data

In order to determine BMDs and BMDLs for gene expression data, BMDExpress was employed (Yang et al., 2007). Briefly, microarray probes with more than one representation on each array were averaged. Analyses were performed on genes that were identified as statistically significant by one-way ANOVA (p  < 0.05) using the four following models: Hill, Power, Linear 1° and Polynomial 2°. The Power model had a power restriction of ≥1. Selection on Linear and Polynomial 2° was based on choosing a model which describes the data with the least complexity. A nested Chi-square test, with cut-off of 0.05, first selects among linear and polynomial models, followed by comparing AIC, which measures the relative goodness of fit. A Hill model was excluded if the “k” parameter of the model was less than 1/3 of the lowest positive dose (18 μg) (Black et al., 2012). Other settings included maximum iterations of 250, confidence level of 0.95, benchmark response (BMR) of 1.349 (number of standard deviation defining BMD) (Yang et al., 2007). For functional classifications and analyses, the resulting BMD datasets were mapped to KEGG pathways with promiscuous probes removed (probes that mapped to multiple annotated genes). BMDs that exceeded the highest exposure dose (162 μg) and that exceeded a goodness-of-fit p-value of 0.1 were removed from the analysis.

Prediction analysis for microarrays

To determine the correlation between gene expression profiles of mice exposed to CBNPs with those of mouse pulmonary disease models, a prediction analysis for microarrays (PAM) (Tibshirani et al., 2002) was conducted in R (R Development Core Team, 2011) using the PAMR library (Hastie et al., 2011). Data for this analysis encompassed 13 mouse lung disease models, and were obtained from the National Centre for Biotechnology Information Gene Expression Omnibus (accession #GSE4231 and #GSE11037). The samples were labelled as belonging to one of three models of lung inflammation: bacterial infection, lung injury and fibrosis, or Th2 response (allergic airway inflammation). Probes with common GENBANK accessions were collapsed to a single measurement for each sample using the mean. Using the common accession numbers, a prediction model using shrunken centroids was estimated. Cross-validation of the nearest shrunken centroid classifier was conducted to identify an appropriate threshold. PAMR implements 10-fold cross-validation. This involves dividing the samples into ten approximately equal-size parts ensuring that the classes are distributed proportionally. Ten-fold cross-validation works by fitting a model on 90% of the samples and then predicting the class labels of the remaining 10%. This procedure is repeated ten times, with each part playing the role of the test samples and the errors on all ten parts added together to compute the overall error. A threshold of 2 was selected, yielding a classifier with 753 GENBANK accessions. The means of the nine CBNP treatment conditions were then classified using the estimated prediction model.

Functional analysis using inflammatory lung disease models

Functional analysis was conducted to establish molecular perturbations that were in common or discrepant between CBNP exposed mice and inflammatory lung disease models. The analysis was conducted on genes that were common between CBNP and each lung disease model, then again for genes that were unique to CBNP, using a cut-off of FDR-adjusted p  < 0.1 and a fold-change > 1.5 for all datasets. The less stringent cut-off was employed for disease models because of the low power in several of the datasets. DAVID Bioinformatics Resources 6.7 was used to identify enriched biological functions from terms with similar genes and biological meaning (Huang et al., 2009a, Huang et al., 2009b). DAVID Biological functions with enrichment scores > 1.3 were considered significant, in accordance with DAVID recommendations (Huang et al., 2009a). Clusters with enrichment scores > 1.3 in our analysis contained at least one gene ontology term or pathway for which the Benjamini-corrected p-value was ≤0.05.

Mining of public data repositories for disease prediction

In order to predict potential disease outcomes of relevance to humans, gene expression profiles were mined against genomic data repositories. Disease prediction analysis was done in NextBio (http://nextbio.com) using the high dose exposure profiles as differentially expressed genes were identified at all time-points for this dose. Data from CBNP exposed mice were compared to curated datasets to identify disease studies with similar gene profiles, gene ranking and consistency. Pairwise gene signature correlations and rank-based enrichment statistics were employed in the calculation of NextBio scores for each disease. The disease that ranked highest in comparison with CBNP exposure was given a score of 100, and the rest of the results were normalized accordingly (Kupershmidt et al., 2010). Meta-analysis was performed using select disease models for mice, as well as for human studies representative of disease state. The analysis identified, ranked and scored all genes and biogroups that were common between the studies according to the scoring method described above for disease prediction (Kupershmidt et al., 2010). Biogroups were filtered for canonical pathways.

Results

Pathway analysis

The rank-based pathway analysis revealed a total of 151, 150 and 106 differentially expressed KEGG pathways on days 1, 3 and 28, respectively. The most affected pathways according to statistical significance were primarily related to inflammation on day 1, to steroid biosynthesis and DNA repair on day 3 and to apoptosis and inflammation on day 28. Significant pathways (p  < 0.05) pertaining to genotoxicity (DNA damage and repair) and inflammatory and immune responses are summarized in Table 1 , along with previously established phenotypes. All significant pathways are presented in Supplemental Table 1. Analysis of the number of common pathways between doses for each time-point revealed that most pathways occurring at lower doses also occur at higher doses. However, the number of significant pathways increased with dose (Fig. 1 ).

Table 1

Summary of significant KEGG pathways (p ≤ 0.05) relating to phenotypes established in Bourdon et al., 2012a, Bourdon et al., 2012b.

	Day 1			Day 3			Day 28
	Genotoxicity-related pathways	Inflammation-related pathways	Phenotype	Genotoxicity-related pathways	Inflammation-related pathways	Phenotype	Genotoxicity-related pathways	Inflammation-related pathways	Phenotype
18 μg	Apoptosis	NOD-like receptor signalling	DNA strand breaks; FPG sensitive sites; neutrophil influx			Neutrophil influx; eosinophil influx; lymphocyte influx; macrophage influx	Apoptosis		Neutrophil influx
54 μg	Nucleotide excision repair; glutathione metabolism	NOD-like receptor signalling; Toll-like receptor signalling; cytokine–cytokine receptor signalling; antigen processing and presentation	DNA strand breaks; FPG sensitive sites; neutrophil influx	Homologous recombination; base excision repair; cell cycle; p53 signalling; nucleotide excision repair; non-homologous end-joining	T Cell receptor signalling; hematopoietic cell lineage; MAPK signalling; antigen processing and presentation; Toll-like receptor signalling; NOD-like receptor signalling; Jak-STAT signalling; B cell receptor signalling	DNA strand breaks; Eosinophil influx; Lymphocyte influx	Apoptosis		DNA strand breaks; eosinophil influx; lymphocyte influx
162 μg	Glutathione metabolism; mismatch repair; cell cycle; nucleaotide excision repair; homologous recombination; non-homologous end-joining; base excision repair; p53 signalling	NOD-like receptor signalling; antigen processing and presentation; cytokine–cytokine receptor interaction; Toll-like receptor signalling; T cell receptor signalling; Hematopoietic cell lineage; Jak-STAT signalling; MAPK signalling; chemokine signalling; B cell receptor signalling	DNA strand breaks; FPG sensitive sites; neutrophil influx	Homologous recombination; glutathione metabolism; cell cycle; mismatch repair; nucleotide excision repair; base excision repair; non-homologous end-joining; p53 signalling	Antigen processing and presentation; NOD-like receptor signalling; hematopoietic cell lineage; cytokine–cytokine receptor interaction; Jak-STAT signalling; Toll-like receptor signalling; chemokine signalling pathway; complement and coagulation cascades	DNA strand breaks; FPG sensitive sites; neutrophil influx; eosinophil influx	Glutathione metabolism; non-homologous end-joining	NOD-like receptor signalling; cytokine–cytokine receptor interaction; Jak-STAT signalling; antigen processing and presentation; B cell receptor signalling; Toll-like receptor signalling; hematopoietic cell lineage	DNA strand breaks; neutrophil influx; lymphocyte influx

Fig. 1

Venn diagrams illustrating overlap of significant pathways (p < 0.05) according to dose, for each post-exposure day (1, 3 and 28 days) in C57BL/6 mice exposed to CBNPs.

BMD and BMDL

EPA BMDS 2.2 BMDs and BMDLs were generated for apical endpoints and RT-PCR data (BMD values for each endpoint and gene presented in Supplementary Table 2; curves are presented in Supplemental Fig. S1). Although many of the apical endpoints and RT-PCR data were not suitable for modelling, BMD and BMDL values generally increased over post-exposure time as expected. The mean BMDs for inflammatory apical endpoints were 0.9, 1.2 and 9.6 μg, and BMDLs were 0.6, 0.9 and 6.5 μg on days 1, 3 and 28, respectively. BMD values for RT-PCR data of genes involved in inflammation tended to be higher than for apical endpoints. Mean BMDs of inflammatory genes were 14.5, 16.7 and 29.0 μg, and mean BMDLs were 10.4, 9.1 and 20.1 μg, on days 1, 3 and 28, respectively.

BMDs and BMDLs were also generated for microarray gene expression profiles using BMDExpress. Minimum BMDs for KEGG pathways relevant to inflammation, KEGG pathways relevant to genotoxicity, for the most sensitive KEGG pathways as well as for apical endpoint data are presented in Table 2 . Minimum BMDs were calculated according to the median of all significant genes for each pathway and the 5th percentile of significant genes of all pathways, in order to increase sensitivity. Even the 5th percentile BMDs tended to be higher than BMDs generated for apical endpoints (Table 2). However, minimum BMDs, representing the most sensitive gene for each relevant pathway, were much more comparable to BMDs of apical endpoints (Table 2).

Table 2

BMD and BMDL for apical endpoints and KEGG pathways relevant to inflammation and genotoxicity.

	Day 1	Endpoint/pathway	Day 3	Endpoint/pathway	Day 28	Endpoint/pathway
Apical endpoints
Min inflammation related BMD	0.5	Neutrophil influx	0.6	Eosinophil influx	3.7	Neutrophil influx
Min inflammation related BMDL	0.3	Neutrophil influx	0.6	Eosinophil influx	2.2	Neutrophil influx
Min genotoxicity related BMD	n/a	n/a	4.5	DNA strand breaks	2.4	DNA strand breaks
Min genotoxicity related BMDL	n/a	n/a	3.0	DNA strand breaks	1.6	DNA strand breaks
KEGG pathways relevant to inflammation
Min median-based BMDa	29.4	Complement and coagulation cascades	32.8	Cytokine–cytokine receptor interaction	77.7	MAPK signalling pathway
Min 5th percentile-based BMDb	18.6	Toll-like receptor signalling pathway	20.2	Hematopoietic cell lineage	70.9	T cell receptor signalling pathway
Mean minimum BMDc	1.8	n/a	3.3	n/a	12.2	n/a
Min median-based BMDLd	17.2	Antigen processing and presentation	19.3	Cytokine–cytokine receptor interaction	54.6	Leucocyte transendothelial migration
KEGG pathways relevant to genotoxicity
Min median-based BMDa	31.9	Apoptosis	29.4	Base excision repair	29.4	Cell cycle
min 5th percentile-based BMDb	23.6	p53 signalling pathway	6.7E-13	Non-homologous end-joining	30.2	p53 signalling pathway
Mean minimum BMDc	18.8	n/a	11.6	n/a	1.88E-06	n/a
Min median-based BMDLd	18.5	Apoptosis	7.8E-06	Non-homologous end-joining	16.5	Cell cycle
All pathways
Min median-based BMDa	17.2	Heparan sulfate biosynthesis	0.7	Keratan sulfate biosynthesis	6.21E-09	Glycine, serine and threonine metabolism
Min 5th percentile-based BMDb	1.9E-13	Amyotrophic lateral sclerosis (ALS)	2.98E-13	N-Glycan degradation	1.00E-01	Alkaloid biosynthesis II
Min median-based BMDLd	2.9	Biotin metabolism	9.11E-10	Keratan sulfate biosynthesis	6.21E-09	Glycine, serine and threonine metabolism

Median-based BMD calculated for each pathway, based on individual BMDs of significant genes.

5th percentile-based BMD calculated for each pathway, based on individual BMDs of significant genes.

Minimum BMD correspond to one gene with the lowest BMD, for each pathway.

Median-based BMDL calculated for each pathway, based on individual BMDLs of significant genes.

Prediction analysis for microarray

PAM was used to compare the Printex 90 gene expression dataset to 13 pulmonary gene expression profiles that represent a range of murine pulmonary disease models (e.g., transgene over-expression, treatments with infectious agents, toxic chemicals and allergens) as described in Lewis et al. (2008), Thomson et al. (2012). The models were classified according to the three following sub-groups: (1) bacterial infection, (2) lung injury and fibrosis, and (3) Th2 response (allergic airway inflammation).

Clustering of the models using PAM is shown in Fig. 2A. Two CBNP exposure conditions (day 28 low and medium doses) did not cluster with other CBNP exposure condition or other disease models, likely due to lack of response. Models of bacterial infection did not cluster with other disease models or CBNP exposure. PAM analysis revealed an association between CBNP exposure, Th2 responses and lung injury/fibrotic responses. Although Th2 response and lung injury/fibrotic responses were more closely associated with one another than with CBNP exposure, PAM analysis revealed that CBNP exposure was more closely related to lung injury/fibrotic responses than to Th2 responses, which is also supported by probability statistics comparing CBNP exposure with each disease sub-group (Fig. 2B).

Fig. 2

PAM using mouse gene expression profiles from 13 lung disease models and CBNP exposed mice. (A) Represents clustering of 753 genes as selected by a PAM threshold cut-off of 2. Sample are classified according to CBNP exposure (purple), Th2 response (yellow), injury and fibrosis (blue) and bacterial infection (red). (B) Represents the probability statistics comparing CBNP exposure for each disease sub-group.

Functional analysis

In order to examine commonalities and discrepancies between disease models and CBNP exposure in more detail, functional analysis was conducted on (1) genes that were in common between CBNP and each disease model and (2) genes that were unique to CBNP. The number of significant genes used for each analysis is presented in Supplemental Table 3. The DAVID biological functions are summarized in Table 3 . This analysis demonstrates that inflammation was common between most models at all time-points (excluding Aspergillus extract). On day 1, commonalities for CBNP exposure were observed with bacterial infection models (i.e., due to the acute phase response) and with injury and fibrosis models (i.e., due to changes in tissue morphogenesis related genes). Day 3 revealed inflammation and cell cycle disturbances in most of the models. However, CBNP responses were more similar to bleomycin-induced lung injury as shown by the high degree of overlapping biological functions on day 3 (Table 3). CBNPs triggered an adaptive immune response on day 28 that was also only apparent in lung injury and fibrosis models.

Table 3

Comparison of CBNP profiles with lung disease models using functional analysis for genes in common (grey) and genes unique to CBNP (black).

Disease prediction

Gene expression profiles from the high dose CBNP-exposed mice vs. control were analysed in NextBio to identify closely related respiratory disease profiles in humans. On all post-exposure days, severe acute respiratory syndrome (SARS), congenital cystic adenomatoid malformation, and injury of lung, were identified as the top three respiratory diseases associated with CBNP exposure. Interestingly, fibrosis was identified as a predicted disease outcome of CBNP exposure that increased considerably with time (e.g., score of 14 on day 1, 35 on day 3 and 45 on day 28). In order to examine the molecular mechanisms that may be involved in fibrosis in more detail, a meta-analysis was completed using curated studies within NextBio that identified fibrosis as a phenotype. Meta-analysis in mouse employed 36 models that included fibrosis-induced by injury with naphthalene, bleomycin and ganciclovir, doxycycline-induced over-expression, and TGF-β over-expression in a variety of mouse models (wild type, inflammation resistant/susceptible). Meta-analysis using CBNP gene expression profiles in mouse ranked 473 canonical pathways and 21,277 genes present in at least one of the studies on select models of pulmonary fibrosis and lung injury (identified in NextBio disease correlation profiles). In order to establish human-relevance, the analysis was repeated using human studies curated in NextBio. Meta-analysis encompassed 4 studies from lung biopsies of patients affected with fibrosis, with intermediate to severe pulmonary hypertension, pneumonia and exacerbation of idiopathic pulmonary fibrosis. Overall, 472 canonical pathways and 15,795 genes were ranked as present in at least one of the studies.

The top ranked pathways and genes for the mouse and human meta-analyses are presented in Table 4 . Interestingly, comparison of fold-ranks between the mouse and human analysis revealed that the most affected pathways were the same in both species. However, the genes that were most perturbed during fibrotic responses were considerably different in CBNP-exposed mice compared to human diseases, with the exception of glycerol-3-phosphate dehydrogenase (GDP1), kruppel-like factor 4 (KLF4), secreted phosphoprotein 1 (SPP1) and ceruloplasmin (CP).

Table 4

Meta-analyses in NextBio using mouse and human profiles in which fibrosis was a phenotype. Values in parentheses represent rank in the opposite species (mouse or human).

Rank 1 (rank 2)	Pathway	Rank 1 (rank 2)	Gene (symbol)
Mouse
1 (1)	Cytokine–cytokine receptor interaction	1 (698)	Cyclin-dependent kinase inhibitor 1a (Cdkn1a)
2 (13)	p53 signalling pathway	2 (158)	Matrix metalloproteinase 14 (Mmp14)
3 (28)	Focal adhesion	3 (2353)	Tissue inhibitor of metalloporteinase 1 (Timp1)
4 (22)	Cell cycle	4 (535)	Metallothionein 2 (Mt2)
5 (2)	Toll-like receptor signalling pathway	5 (786)	Thymidine kinase 1 (Tk1)
6 (15)	Cell communication	6 (2936)	Osteoglycin (Ogn)
7 (29)	ECM receptor interaction	7 (–)	elastin (Eln)
8 (16)	Cell cycle G1 to S reactome	8 (594)	Suppressor of cytokine signalling 3 (Socs3)
9 (19)	Cell cycle G1/S checkpoint	9 (273)	Lysloxidase (Lox)
10 (3)	Hematopoitic cell lineage	10 (491)	Serum amyloid a3 (Saa3)
Human
1 (1)	Cytokine–cytokine receptor interaction	1 (1001)	Mucin 4 (MUC4)
2 (5)	Toll-like receptor signalling pathway	2 (652)	Glutathione S-transferase alpha 1 (GSTA1)
3 (10)	Hematopoitic cell lineage	3 (56)	Kruppel-like factor 4 (KLF4)
4 (17)	Erythrocyte differentiation pathway	4 (4161)	Golgi membrane protein 1 (GOLM1)
5 (14)	IL-1R pathway	5 (48)	Glycerol-3-phosphate dehydrogenase 1 (GPD1)
6 (12)	Calcineurin Nf at signalling	6 (1366)	Spermatogenesis associated, serine-rich 1 (SPATS2L)
7 (26)	Local acute inflammatory response pathway	7 (76)	Secreted phosphoprotein 1 (SPP1)
8 (27)	Cytokines and inflammatory pathway	8 (620)	Chemokine (C-X-C motif) ligand 10 (CXCL10)
9 (21)	MMP cytokine connection	9 (89)	Ceruloplasmin (CP)
10 (30)	Msp/RON receptor signalling pathway	10 (–)	Vanin 2 (VNN2)

Discussion

It is now widely accepted that toxicity is preceded by, and accompanied by, transcriptional changes, thus providing molecular signatures of direct and indirect toxic effects (Auerbach et al., 2010, Fielden et al., 2011, Gatzidou et al., 2007). It is hypothesized that toxicogenomic profiling can be used as a screening tool to prioritize the specific assays that should be conducted from the standard battery of tests, thus minimizing animal use, cost and time (Dix et al., 2007). Moreover, global analyses of transcriptional changes provide a wealth of information that can be used to identify putative modes of action and to query relevance to human adverse health outcomes (Currie, 2012). This type of approach is the general premise of the widely supported paradigm outlined in ‘Toxicity Testing in the 21st Century’ (National Academy of Sciences, 2007). However, substantive work demonstrating the ability of gene expression profiles to identify hazards, to assess risk of exposure via quantitative dose–response analysis, and to identify adverse outcomes associated with specific modes of action is required before these endpoints can be used in HHRA. The present study applies pathway- and network-based approaches, BMD modelling, and disease prediction tools to gene expression data to explore the relationship between apical endpoints and transcriptional profiles. The work investigates the potential utility of gene expression profiling in determining hazard and mode of action of NPs, in characterizing dose–response relationships and in predicting the relevance of these findings to potential disease-outcomes and human health effects for HHRA.

Hazard identification, mode of action and prediction of apical endpoints

The utility of gene expression profiling in hazard identification has been examined for a limited number of chemicals, including dibutyl phthalate and acetaminophen (Euling et al., 2011, Kienhuis et al., 2011, Makris et al., 2010). Toxicogenomic profiles of alachlor exposure in rat olfactory mucosa (Genter et al., 2002) and dimethylarsenic (DMA) exposure in human cultured bladder cells and rat bladder epithelium (Sen et al., 2005, US EPA, 2005) have also provided useful information for two final assessments of acetochlor and arsenicals (US EPA, 2004, US EPA, 2006). Our data demonstrate that gene expression profiles can also be viewed as effective predictors of the biological effects of CBNP exposure. For example, inflammatory responses manifested at the gene expression level and detected using DNA microarrays and classified in this work using KEGG pathway analyses and previously in the same mice using ingenuity pathway analysis (Bourdon et al., 2012a) are entirely consistent with the observed pulmonary influx of inflammatory markers (e.g., neutrophils, eosinophils and lymphocytes). The number of genes perturbed and the magnitude of expression changes in these pathways correlates with dose and time. In addition, observed transcriptomic changes associated with perturbations of cell cycle networks, alterations of non-homologous end-joining, and p53 signalling support the sustained genotoxicity observed in the mice, although dose and time correlations were not as apparent (e.g., levels of DNA strand breaks remained relatively constant at the two highest exposure doses (Bourdon et al., 2012b) whereas induction of DNA repair genes decreased with dose and time). The transcriptomic changes associated with alterations in glutathione metabolism and free radical scavenging correlate with induction of DNA formamidopyrimidine DNA glycoslase (FPG) sensitive sites (an indicator of oxidative DNA damage) early after the exposure. The persistence of this response is an indication of an adaptive response to oxidative stress in the lungs of the mice. Interestingly, CBNP-induced alterations in gene expression profiles also revealed a pulmonary acute phase response and unexpected changes in lipid homeostasis, which were subsequently supported by measured decreases in plasma high density lipoprotein (HDL) (Bourdon et al., 2012a). The strong association between CBNP-induced gene expression profiles and apical endpoints collectively support the use of toxicogenomics for hazard identification of NMs, and perhaps more importantly, for highlighting unexpected adverse outcomes. Moreover, ongoing work within the Organization for Economic Co-operation and Development (OECD) is actively developing adverse outcome pathways (AOP) approaches that are expected to provide tangible methods by which systems biology endpoints can be used in human health risk assessment. Toxicogenomics data that examine responses over dose and time in a variety of tissues can be very useful for such applications, as illustrated for CBNP exposure in Fig. 3 . Overall, our data suggest that gene expression profiles can be effectively used to identify putative mode(s) of action and hazards of NP exposure, in the absence of phenotypic data.

Fig. 3

Use of toxicogenomic data and observed phenotypes in an adverse outcome pathway.

Dose–response characterization

In addition to identification of hazard, it has been suggested that gene expression profiles may be useful for quantitative assessment (e.g., establishment of reference doses) of responses related to both cancer and non-cancer endpoints (Thomas et al., 2007). Benchmark doses are generally considered more informative than the no observable adverse effect level (NOAEL) in deriving reference doses as they are based on the entire dose–response relationship (Crump et al., 1995). Because alterations in gene expression can be initiated in the absence of biological effects (e.g., adaptive or stress response pathways effective in mitigating toxic effects), it is expected that reference doses for genomics endpoints may be too sensitive for use in HHRA. However, previous analyses of 5 chemicals (i.e., 1,4-dichlorobenzene, propylene glycol mono-t-butyl ether, 1,2,3-trichloropropane, methylene chloride and naphthalene) showed that median BMD and BMDLs for the most sensitive pathways and GO categories were highly correlated with BMD and BMDLs of cancer and non-cancer endpoints (Thomas et al., 2011, Thomas et al., 2012).

In the current study, rather than choosing the most sensitive (i.e., lowest) BMDs, we focussed on the analysis of pathways that were specific to biological outcomes observed in the mice (i.e., phenotypically anchored), and calculated BMDs for these relevant genes and pathways. The pathway-based BMDs and BMDLs calculated here for relevant pathways were actually less sensitive (i.e., higher BMDs) than those of the observed apical endpoints. However, the mean of the minimum BMDs and BMDLs across all the pathways that we assigned as relevant to the apical endpoints (i.e., corresponding to the most sensitive genes within the relevant pathways) were similar to those of relevant apical endpoints. Median BMDs and BMDLs for the most sensitive pathways also correlate more closely with apical endpoints even though the pathways were not necessarily relevant to these endpoints. This finding supports previous examples demonstrating a 1:1 correlation between BMDs for gene expression and apical endpoints (Thomas et al., 2011, Thomas et al., 2012). These data indicate the potential utility of using gene expression profiles in determining acceptable exposure limits for NPs. In order to determine the specific utility of pathway derived BMDs in HHRA, it will be necessary to establish a comprehensive catalogue of pathways that are actually perturbed in the event of specific adverse effects.

Disease prediction and relevance to humans

Perhaps the principal motivation for including gene expression profiling in HHRA is the wealth of information that can be used to identify key events that are correlated with adverse outcomes that are relevant to human disease, and moreover can be used to predict the likelihood of a human disease. Identification of key events at the transcriptional level can facilitate the identification of processes that are critical for disease initiation and progression, thus allowing information from animal experiments to be queried and used for extrapolation to human scenarios (Edwards and Preston, 2008). Comparison of our data with specific models of lung disease, including bacterial infection, airway hypersensitivity and lung injury revealed that CBNPs induced responses that were more closely related to lung injury and fibrosis than to other models. This finding was further supported by comparison of the expression profiles of CBNP exposed mice to those of curated studies of animals and humans exhibiting a myriad of pulmonary disease phenotypes. This analysis demonstrates that CBNP exposure perturbs genes that are known to be involved in tissue injury and fibrosis in mice. Although it is unclear if CBNP exposure would result in the same gene expression profile in humans, similar pathways including many involved in fibrotic responses were found in both mice and humans (52% of the top 50 pathways found were common between mouse and human). Despite concordance of pathways, the top ranked genes differed considerably between both species. However, many of the genes found in mice and humans had similar functions, including inflammatory and acute phase responses (e.g., Saa3, Socs3 and Mt2 in mice and CP, VNN2 and CXCL10 in humans), cell cycle progression (Cdkn1a in mice and KLF4 in humans) and bone and tissue modelling (Mmp14, Timp1, Eln and Ogn in mice and SPP1 in humans). Thus, despite discordance in the gene expression profiles between species, the similar functions of top ranked genes and concordance between pathways supports the likelihood of similar responses in the event of CBNP exposure in humans. In addition, fibrosis has been identified as an outcome of exposure to various particles and NPs in animals (Bermudez et al., 2004, Shvedova et al., 2008), including Printex 90 (e.g., 28-day nose only inhalation in Wistar WU rats) (Bellmann et al., 2009), as well as in humans (Lkhasuren et al., 2007, Wang and Christiani, 2003). The process of pulmonary fibrosis is closely related to progression of carcinogenic outcome (Hubbard et al., 2000). These data demonstrating very similar fibrotic pathways in mice and humans and a significant overlap with CBNP-induced gene expression changes thus support the use of pathway-based approaches in identifying molecular mechanisms of disease onset and progression, and using gene expression profiles to support HHRA.

Conclusions and recommendations

This study confirms several key elements that are necessary for the application of gene expression profiling for HHRA of toxicant exposures in general. First, transcriptional profiles can effectively predict the biological effects of chemical exposures. Specifically, in the absence of data for any apical endpoints, our data would have suggested that mice exposed to CBNPs exhibit an inflammatory response, oxidative stress, DNA damage and perturbations in cholesterol metabolism. Second, a comparison of BMDs and BMDLs of relevant pathways and apical endpoints confirms that minimum pathway BMDs and BMDLs are in the same range as those of apical endpoints. Third, that expression profiles can be fairly easily mined to identify potential adverse outcomes (i.e., diseases) that are relevant to humans, and might reasonably be expected to occur in humans exposed to substances that elicit specific gene expression patterns in experimental animals. We believe that our work constitutes a significant step towards the ultimate recognition of toxicogenomic endpoints for routine assessment of human health risk.

Gene expression profiling offers a promising approach to decipher the largely unknown hazards of NP exposure. Due to the unique properties of NPs, powerful technologies that can assess a multitude of adverse outcome possibilities will be required to elucidate their modes of action and potential impacts on human health within a time-frame that is suitable for prompt regulatory decision making. This same premise should hold true for any new chemical products, for which toxicity is largely or completely unknown. In order to establish a strong foundation for the integration of gene expression profiling into HHRA, it will be necessary for the approach employed here to be applied to a variety of additional chemicals/particles that span a wide range of toxicological potencies and modes of action, and using a variety of experimental designs (e.g., multiple doses and time-points). As our knowledge of molecular pathways, and of the diverse tools used to decipher their biological significance, dose–response characteristics and relevance to human disease continues to grow, we anticipate that toxicogenomics will become increasingly useful in assessing the toxicological hazards of a wide range of test articles, and by extension, for HHRA.

Conflict of interest statement

None.