Integration of metabolic activation with a predictive toxicogenomics signature to classify genotoxic versus nongenotoxic chemicals in human TK6 cells.

The use of integrated approaches in genetic toxicology, including the incorporation of gene expression data to determine the molecular pathways involved in the response, is becoming more common. In a companion article, a genomic biomarker was developed in human TK6 cells to classify chemicals as genotoxic or nongenotoxic. Because TK6 cells are not metabolically competent, we set out to broaden the utility of the biomarker for use with chemicals requiring metabolic activation. Specifically, chemical exposures were conducted in the presence of rat liver S9. The ability of the biomarker to classify genotoxic (benzo[a]pyrene, BaP; aflatoxin B1, AFB1) and nongenotoxic (dexamethasone, DEX; phenobarbital, PB) agents correctly was evaluated. Cells were exposed to increasing chemical concentrations for 4 hr and collected 0 hr, 4 hr, and 20 hr postexposure. Relative survival, apoptosis, and micronucleus frequency were measured at 24 hr. Transcriptome profiles were measured with Agilent microarrays. Statistical modeling and bioinformatics tools were applied to classify each chemical using the genomic biomarker. BaP and AFB1 were correctly classified as genotoxic at the mid- and high concentrations at all three time points, whereas DEX was correctly classified as nongenotoxic at all concentrations and time points. The high concentration of PB was misclassified at 24 hr, suggesting that cytotoxicity at later time points may cause misclassification. The data suggest that the use of S9 does not impair the ability of the biomarker to classify genotoxicity in TK6 cells. Finally, we demonstrate that the biomarker is also able to accurately classify genotoxicity using a publicly available dataset derived from human HepaRG cells.

INTRODUCTION

The Health and Environmental Sciences Institute's (HESI) Technical Committee for the Application of Genomics to Mechanism‐Based Risk Assessment is coordinating an international effort to develop a genomic biomarker to classify chemicals as either genotoxic or nongenotoxic for the purposes of human health risk assessment and pharmaceutical drug development [Goodsaid et al., 2010; Li et al., 2015]. The foundation of this initiative, which is presented in an accompanying article, is a reference database that consists of gene expression profiles from human TK6 cells exposed to twenty‐eight model compounds [Goodsaid et al., 2010; Li et al., 2015]. The reference database represents chemicals exemplifying a wide spectrum of well‐characterized mechanisms of action, both genotoxic and nongenotoxic in nature, at one concentration and one time point (4 hr; Supporting Information, Table 1). Gene expression profiles were generated using Agilent microarrays to build the reference database [Goodsaid et al., 2010; Li et al., 2015]. Ultimately, a 65‐gene classifier was developed that included various genes associated with a DNA damage response and that was derived by applying the nearest shrunken centroids (NSC) method [Li et al., 2015; Tibshirani et al., 2002]. This classifier was named the TGx‐28.65 classifier, indicating that 28 reference compounds were used in the training set to derive a 65‐gene signature (Supporting Information, Table 2). Our companion article demonstrates that the TGx‐28.65 classifier, in which the P53‐signaling pathway is significantly overrepresented, can differentiate genotoxic from nongenotoxic compounds with 100% accuracy (refer to Li et al., 2015, for working definition of genotoxicity). The 28 chemicals used in the development of the HESI genomic biomarker include compounds that do not require metabolic activation in order to produce DNA damage. This limits the suitability of the assay for environmental chemicals, a larger proportion of which may require metabolic activation to exert their genotoxic effects [Shimada et al., 2013]. Therefore, to broaden the utility of this toxicogenomics (TGx) biomarker, we asked whether the TGx‐28.65 gene signature could be used in the presence of a metabolic activation system (MAS) for chemicals requiring metabolic activation.

In this study, we investigate the performance of the TGx classifier to accurately classify compounds requiring metabolic activation as genotoxic or nongenotoxic in the presence of rat liver S9 microsomal fraction. We also investigate the concentration‐response relationship of the classifier, and its correlation with induction of genetic damage (micronuclei: MN), apoptosis, and cytotoxicity for chemicals requiring metabolic activation. To do this, we exposed TK6 cells, a human lymphoblastoid cell line, to increasing concentrations of two genotoxic chemicals, benzo[a]pyrene (BaP) and aflatoxin B1 (AFB1), and two nongenotoxic chemicals, dexamethasone (DEX) and phenobarbital (PB), in the presence of rat liver S9 for 4 hr. Agilent gene expression profiles immediately postexposure, in addition to 4 hr and 20 hr postexposure, were aligned against the existing training set, and the TGx‐28.65 classifier was used to assign a genotoxicity score for each chemical concentration and time point. The classifier was interpreted in the context of relative survival (RS), apoptosis, and MN frequency collected at the 24 hr time point by flow cytometry. The data clearly demonstrate the potential of the classifier to predict genotoxicity in the presence of S9 and suggest that the strength of the classifier is correlated with genotoxicity for BaP and AFB1. The classifier was further validated using a published dataset from more metabolically competent human HepaRG cells exposed to genotoxic and nongenotoxic chemicals [Doktorova et al., 2013]. In this external validation, the classifier was 100% accurate in classifying genotoxic and nongenotoxic chemicals, indicating that the signature is suitable for use with other cell lines and microarray platforms.

MATERIALS AND METHODS

Chemicals

Benzo[a]pyrene (CAS No. 50‐32‐8), aflatoxin B1 (CAS No. 1162‐65‐8), dexamethasone (CAS No. 50‐02‐2), and phenobarbital (CAS No. 50‐06‐6) were purchased from Sigma‐Aldrich (St. Louis, MO). All chemicals were dissolved and diluted in dimethylsulfoxide (DMSO), which served as the vehicle control for cellular exposures.

Cell Culture

TK6 cells, a human lymphoblastoid cell line, were obtained from American Type Culture Collection (ATCC# CRL‐8015; ATCC, Manassas, VA). TK6 cells express wild type P53, demonstrate rapid growth in suspension and are often used in standard genetic toxicology studies; they are considered to be a suitable cell line for use in the in vitro mammalian MN test [Islaih et al., 2005; OECD, 2014; Godderis et al., 2012]. TK6 cells were maintained as described elsewhere [Ellinger‐Ziegelbauer et al., 2009]. Briefly, cells were cultured and maintained in RPMI 1640 medium containing 10% heat inactivated horse serum, in addition to 0.1% pluronics, sodium pyruvate, and antibiotics (penicillin at 20 units/ml and streptomycin at 20 µg/ml) at 37 ± 1°C and 6 ± 1% CO2 in air. Immediately before chemical exposure, cells were seeded at a density of 4 (± 0.5) × 105 cells/ml in 12‐well plates with a final volume of 3 ml per well. Chemical exposures were conducted in the presence of 1% 5,6 benzoflavone‐/phenobarbital‐induced rat liver S9 (BF/PB‐induced S9) (Moltox, Boone, NC) with NADPH generating system cofactors. The in vitro cellular exposures and associated cytotoxicity/genotoxicity measures were performed by Integrated Laboratory Systems, Inc. (ILS; Research Triangle Park, NC).

Range Finder Studies

Range finder studies were conducted for each of the four test chemicals using a 4 ± 0.5 hr exposure, followed by a 20 ± 0.5 hr recovery (i.e., 24 hr time point). Exponentially growing TK6 cells were exposed to four to eight different concentrations of each chemical as follows: BaP (0.125–10 μg/ml), AFB1 (0.025–0.2 μM), DEX (7.813 μg/ml–10 mM (3,924.6 μg/ml)), and PB (0.2–21 nM), in addition to the appropriate negative and positive controls. Negative controls for the chemical exposures consisted of media only (negative control (NC)) and DMSO vehicle controls (VC), in the presence of S9 (NC (+S9) and VC (+S9)). In addition, TK6 cells exposed to 24 µg/ml cisplatin served as the positive control. As cisplatin is a direct‐acting chemical and does not require metabolic activation, cisplatin treatment, and its corresponding negative controls were performed in the absence of S9 (NC (‐S9) and VC (‐S9)). Exposures were done in duplicate for BaP and as single replicates for AFB1, DEX, and PB for the range‐finder studies. All exposures were done in the presence of 1% BF/PB‐induced rat liver S9 with NADPH generating system cofactors. Following the 4 hr exposure, cells were collected by centrifugation and rinsed with PBS before the replacement of fresh media. Cells were then placed back in the incubator for an additional 20 ± 0.5 hr. At the end of the 24 hr time period, cells were harvested for cytotoxicity and MN analysis.

Cytotoxicity, Apoptosis, and MN Frequency

The In Vitro MicroFlow kit (Litron Laboratories, Rochester, NY), a high content flow‐cytometry based assay, was used to determine the cytotoxic and genotoxic effects of the test chemicals in human TK6 cells. This assay was used to determine the percent relative survival (% RS) compared with control cells, and the percentage of apoptotic/necrotic cells, while simultaneously scoring MN in the same cell population. Reagent preparation, cell harvest, sequential staining and flow cytometric analysis were performed according to the In Vitro MicroFlow Instructional Manual (also described in [Recio et al., 2012]). The data were collected using a Becton‐Dickinson FACSCalibur 2 laser 4‐color instrument (Becton Dickinson, San Jose, CA). Briefly, RS was determined using an absolute counting technique in which 6 μm latex microspheres were added to the cell preparation as “counting” beads, which serve as internal standards for flow cytometry. RS was then calculated using the intact viable nuclei to bead ratios in treated versus control cells. Apoptotic and necrotic cells were also identified by flow cytometry, as ethidium monoazide (Nucleic Acid Dye A contained in the kit) crosses the outer membrane of necrotic and apoptotic cells. Finally, MN frequency was simultaneously scored in treated and control cells by analyzing 20,000 (±2,000) cells from each sample using the double staining procedure outlined in the instruction manual for the In Vitro MicroFlow kit.

Concentration Selection Criteria for Definitive Studies

Range finder studies were conducted to select three concentrations per chemical (low, medium, and high) for subsequent definitive studies. A low concentration was selected that induced expression of at least one of three selected genes: Activating Transcription Factor 3 (ATF3; a stress response gene), Cyclin‐Dependent Kinase Inhibitor 1a (CDKN1A; a p53 target gene), and Growth Arrest and DNA Damage Inducible, alpha (GADD45A; a stress response gene), compared with controls using real‐time quantitative RT‐PCR in the absence of cytotoxicity and minimal induction of MN at 24 hr. The high concentration was selected to cause approximately 55 ± 5% cytotoxicity (i.e., 40–50% relative survival), as recommended in OECD test guideline 487 [OECD, 2014; Kirkland et al., 2011]. Finally, the medium concentration was selected to induce MN (more than twofold induction; P < 0.05) at the 24 hr time point and was between the low and the high concentrations. We point the reader to our companion article for a more complete description of the utility of these genes in establishing the appropriate concentration for subsequent testing [Li et al., 2015]. In the absence of observed cytotoxicity or genotoxicity, a second range finder was conducted to select the top concentration using the following selection criteria: 1 mM or 0.5 mg/ml, and 10 mM or 5 mg/ml, whichever is lower in both cases, whenever solubility in the vehicle or culture medium or observed cytotoxicity was not a limiting factor. This selection criterion was based on the revised and former ICH Guidelines for the Guidance on Genotoxicity Testing and Data Interpretation for Pharmaceuticals Intended for Human Use [ICH, 1996; ICH, 2011]. These guidelines were used for DEX dose selection only, as there was no observed cytotoxicity or robust target gene expression changes following two range finder studies. As such, four additional concentrations of DEX were tested, including 1 mM, 5 mM, 7.5 mM, and 10 mM concentrations; however, the 10 mM dose heavily precipitated out of solution and thus 7.5 mM was used as the top concentration for this chemical.

Definitive Studies

TK6 cells were exposed to three concentrations of each chemical as follows: BaP (0.45 μg/ml, 1.4 μg/ml, 10 μg/ml), AFB1 (0.025 μM, 0.075 μM, 0.1 μM), DEX (0.63 mM, 1 mM, 7.5 mM) and PB (1 mM, 3.2 mM, 10 mM), with a minimum of three technical replicates. The highest concentration of BaP, AFB1, and DEX were also tested in the absence of S9 alongside VC (–S9) at 4 hr, 8 hr, and 24 hr. All experiments were run in parallel with cisplatin‐treated positive controls. Cisplatin was also tested in the presence of S9. Because the 4 hr time point was less effective in predicting genotoxicity than the 8 hr and 24 hr time points (for BaP and AFB1, which were run first), this time point was not performed for PB. Separate plates were used for exposures to negative and positive controls (NC (±S9), VC (±S9), and PC‐24 (−S9)). Time points were named as follows: 4 hr = 4 hr of exposure followed by immediate sample collection; 8 hr = 4 hr of exposure, media replaced and sampled 4 hr later; 24 hr = 4 hr of exposure, media replaced, and sampled 20 hr later. The definitive exposures were conducted as described in the range finder studies described above. At the end of each 4 hr exposure, the exposed and control TK6 cells for each time point were removed by centrifugation and the cells were rinsed in PBS. The 4 hr samples were immediately harvested and flash frozen for RNA collection and a portion was used for measurement of cytotoxicity. The 8 hr and 24 hr samples were resuspended in 3 ml of media and re‐incubated for an additional 4 ± 0.5 and 20 ± 0.5 hr, respectively, before being harvested for RNA and cytotoxicity/MN assays. Cells dedicated for the purposes of RNA extraction were collected from each well (4 ± 0.5 × 106 cells), pelleted by centrifugation and flash frozen in liquid nitrogen. Cells were stored at −80oC until RNA extraction was performed. At the end of each time point, 100 μl of each cellular sample was transferred to a 96‐well plate to assess cellular viability using the CellTiter96® AQueous One Solution Cell Proliferation Assay (Promega, Madison WI), which measures cell viability and proliferation by means of a colorimetric assay with 3‐(4,5‐dimethylthiazol‐2‐yl)‐5‐(3‐carboxymethoxyphenyl)‐2‐(4‐sulfophenyl)‐2H‐tetrazolium (MTS), an MTT‐equivalent assay, following the manufacturer's instructions. This cytotoxicity endpoint was measured in exposed and control cells at the 4 hr, 8 hr, and 24 hr time points. The remaining volume of cells from the 24 hr time point was used to assess cytotoxic and genotoxic effects using a flow cytometry‐based high content cytotoxicity and MN assay described above. The cytotoxicity results from the MTS assays can be found in the Supporting Information (Supporting Information Figs. S1–S4 for BaP, AFB1, DEX and PB, respectively). Any further reference to cytotoxicity will signify RS data.

Total RNA Extraction

An RNeasy Mini Kit (Qiagen, Toronto, ON, Canada) was used to isolate and purify total RNA using the spin technology protocol with an on‐column DNase I digestion, following the manufacturer's instructions. Total RNA was quantified and assessed for purity and integrity using a NanoDrop® ND‐1000 spectrophotometer and an Agilent 2100 Bioanalyzer. Only high quality total RNA samples, with RNA Integrity Numbers (RINs) ranging from 8.6 to 10 and A260/280 absorbance ratios ≥2.0 were used for gene expression analysis.

Quantitative Real‐Time RT‐qPCR

First‐strand cDNA was synthesized from 325 to 650 ng of total RNA using SuperScript™ III First‐Strand Synthesis System with the oligo‐dT primer (Invitrogen, Burlington, ON, Canada). Quantitative real‐time RT‐qPCR was performed using TaqMan® Gene Expression Assays (TaqMan® MGB probes, FAM™ dye‐labeled) on a Bio‐Rad CFX96™ Real‐Time PCR Detection system. Gene expression for ATF3, CDKN1A, and GADD45A was measured using commercially available TaqMan® gene expression assays (Life Technologies Inc., Burlington, ON, Canada), according to the manufacturer's instructions. TaqMan® Gene Expression Assay IDs were as follows: Hs00231069_m1 (ATF3); Hs00355782_m1 (CDKN1A); Hs00169255_m1 (GADD45A) [Ellinger‐Ziegelbauer et al., 2009]. A human GAPDH endogenous control was used as the reference gene for these gene expression studies. GAPDH demonstrated stable expression for each chemical across treated and control samples at each time point. RT‐qPCR analysis was performed for all four genes, including the three genes of interest (GOI) and the reference gene (REF), at all 3 time points (4 hr, 8 hr, and 24 hr; n = 3–6). Normalized expression (NE) was calculated using the following equation, where E is equal to the PCR reaction efficiency and C(t) is the crossing threshold [Simon, 2003]: Normalized expression values were analysed using a one‐way ANOVA independently for each time point. Pairwise comparisons were conducted and P values were adjusted using the Tukey method [Miller, 1981]. Results were back‐transformed and the fold changes linearized. The delta method was used to obtain the standard errors for the back transformed results. RT‐qPCR gene expression data were used to confirm appropriate selection of concentration for microarray analysis (described in [Li et al., 2015]). Specifically, the highest concentration of each genotoxic chemical had to increase the expression of at least one of these genes across all biological replicates (P < 0.05 and ≥1.5 times up regulation of gene expression), while exhibiting 55 ± 5% cytotoxicity and induction of MN at the 24 hr time point.

Microarray Analysis

Transcriptome measurements were performed using a two‐color dye swap design [Kerr and Churchill, 2001]. Microarray analysis was conducted on all time points from the definitive studies (n = 3 per concentration and time point), along with the corresponding pooled vehicle and positive controls. Agilent Low‐Input Quick Amp Labeling kits were used to ultimately generate cyanine‐3 and cyanine‐5 labeled cRNA from 100 ng of total RNA, according to the manufacturer's instructions (Agilent Technologies, Mississauga, ON, Canada). Labeled cRNA samples (325 ng of cyanine‐3 labeled cRNA and 325 ng of cyanine‐5 labeled cRNA) were co‐hybridized to Agilent SurePrint G3 Human GE 8x60K oligonucleotide microarrays (Agilent Catalogue No: G4851A, Agilent Microarray Design ID: 028004; Agilent Technologies, Mississauga, ON, Canada) at 65°C for 17 hr. Slides were washed according to the manufacturer's specifications and scanned using an Agilent DNA Microarray Scanner at 3 μm resolution. Agilent Feature Extraction (version 11.0.1.1) was used for data extraction from the generated image files and to generate QC reports for each array.

Statistical and Bioinformatic Analyses

The gene expression results, both normalized and raw intensity values, for all microarray data sets are compliant with the minimal information about a microarray experiment (MIAME) standards [Brazma et al., 2001] and have been deposited in the National Centre for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database under accession number GSE51175. Detailed methodology regarding statistical and bioinformatics analyses are outlined in Li et al. (2015; companion article). Briefly, the intensity values for distinct Agilent probe IDs, which represent the same Genbank accession number (gene symbol), were merged to generate the error‐weighted average for each gene. Hierarchical clustering was performed using the hclust function in R (www.r-project.org). Agglomerative clustering was based on complete linkage with 1 − ρ° as the dissimilarity metric, where ρ° is the error‐weighted Pearson's correlation coefficient for two treatments or two genes [Li et al., 2015]. The Nearest Shrunken Centroids (NSC) method, implemented in the pamr function of R (www.bioconductor.org), was employed to generate the classifier (TGx‐28.65), which now consists of 64 genes due to an annotation update (Supporting Information Table 2; Li et al., 2015 (companion article)). The NSC method was also used to classify the gene expression profiles of the test agents in this study by examining them for similarities with the transcriptome profiles of the reference chemicals in the database using statistical and bioinformatics tools [Tibshirani et al., 2002]. In brief, the standardized centroid (SC) is computed by applying the NSC method for each class of the training set and is the mean expression level for each gene in a class divided by its within‐class standard deviation. For each class, the SC is shrunken in the direction of the overall centroid to create the NSC. Experimental samples are then classified through comparison of their gene expression profile to the class of NSCs. Sample classification is accomplished by assigning it to a class that is closest to it in squared distance such that the probability of class membership is greater than 0.90 [Li et al., 2015].

Analysis of the TGx‐28.65 Classifier in Human Liver HepaRG Cells

In a study by Doktorova et al. [2013], human hepatoma‐dervied HepaRG cells were exposed to 15 prototypical compounds belonging to three toxic classes: (i) genotoxic carcinogens [AFB1; 4‐(methylnitrosamino)‐1‐(3‐pyridyl)‐1‐butanone (NNK); 2‐nitrofluorene (2NF); BaP; cyclophosphamide (CYCLO)], (ii) nongenotoxic carcinogens [methapyrilene hydrochloride (MPH); piperonylbutoxide (PIPB), Wy‐14643 (WYE), phenobarbital sodium (SPB), 12‐O‐tetradecanoylphorbol‐13‐acetate (TPA)]; and (iii) noncarcinogens [nifedipine (NIF); clonidine (CND); d‐mannitol (MAN); tolbutamide (TOL); diclofenac sodium (SDF)] for 72 hr to study the transcriptomic responses. Raw CEL files from this study generated using Affymetrix U133 Plus 2.0 GeneChips (GEO accession number GSE40117) were downloaded to test the TGx‐28.65 biomaker in a human cell line with inherent metabolic capabilities [Doktorova et al., 2013]. The human HepaRG microarray data were RMA normalized in R using the ReadAffy function in the Affy package. Using 53 genes from the TGx‐28.65 classifier that were present on both microarray platforms, the shrunken centroid classifier developed using the TK6 data was applied to predict genotoxicity following exposure to the aforementioned chemicals in human HepaRG cells.

RESULTS

Benzo[a]pyrene

A range finder study was conducted in TK6 cells exposed to five concentrations of BaP (0.125–10 μg/ml), alongside positive, negative, and vehicle controls (data not shown). The following concentrations were chosen for subsequent gene expression microarray studies: 0.45, 1.4, and 10 μg/ml. Relative survival, apoptosis and MN for the definitive studies are shown in Figure 1. At the low concentration, there were slight indications of cytotoxicity (70.9% RS) and genotoxicity (1.6‐fold increase in MN). Although there was no induction of ATF3 at the lowest concentration of BaP, CDKN1A was up‐regulated at all three time points, in addition to GADD45A at the 24 hr time point (Fig. 2). The mid concentration resulted in a significant decline in RS to 56.5% and a 3.5‐fold induction of MN (Fig. 1). The high concentration caused a significant decline in RS to 47.7% and a 6.1‐fold increase in the formation of MN (Fig. 1). Significant induction of ATF3, CDKN1A, and GADD45A at 4 hr, 8 hr, and 24 hr was noted for the mid and high concentrations (Fig. 2). A heat map for cells treated with BaP using the TGx‐28.65 classifier showed both concentration‐ and time‐dependent effects of BaP and the positive cisplatin control (Figs. 3A and 3E, respectively). Analysis of the TGx‐28.65 gene set in BaP‐treated TK6 cells against HESI's database of genomic profiles of known direct‐acting genotoxic and nongenotoxic compounds [Li et al., 2015] predicts that BaP is genotoxic at the mid and high concentrations in the presence of S9 metabolic activation for all time points. The NSC classification probabilities using the TGx‐28.65 classifier for BaP treatment are shown by the prediction bar in Figure 3A. The data suggest that the classifier is quantitative in nature, as the low concentration of BaP was classified as nongenotoxic at 4 hr and was unclassified at 8 hr and 24 hr, consistent with our observation that this concentration did not cause overt signs of cytotoxicity or genotoxicity. Hierarchical cluster analysis would lead us to classify all concentrations of BaP as a genotoxic compound by comparison of its transcriptome profile in human TK6 cells using the TGx‐28.65 classifier when aligned against the database of genotoxic and nongenotoxic reference chemicals, except for the low concentration at 4 hr, which clusters with the nongenotoxic agents (Fig. 4). Finally, the 4 hr and 24 hr BaP data were published in part in a case study on BaP to demonstrate the potential usefulness of these data in human health risk assessment [Moffat et al., 2015].

Figure 1

Cytotoxicity and genotoxicity measurements in human TK6 cells (n = 3) following exposure to BaP using the In Vitro Microflow kit (Litron Laboratories). Percent relative survival (▪), percentage of apoptotic/necrotic cells (●), and percentage of MN (▴) are depicted at the 24 hr time point (4 hr exposure + 20 hr recovery; n = 3). *P < 0.001 compared with the vehicle control. Error bars depict the standard error. % RS is a measure of the relative survival, % apoptosis is a measure of apoptotic cells, and % MN is a measure of micronucleus induction.

Figure 2

Normalized gene expression of targeted genes, including ATF3, CDKN1A, and GADD45A, in human TK6 cells by TaqMan RT‐qPCR following BaP exposures at 4 hr, 8 hr, and 24 hr. Normalized expression was calculated using human GAPDH as the reference gene (n = 3–6). Vehicle control +S9 (VC (+S9)) served as the negative control for BaP treatment *P < 0.05, **P < 0.01, and ***P < 0.001 compared with the vehicle control. Error bars depict standard error.

Figure 3

Heat maps of the TGx‐28.65 gene signature following exposure to: (A) BaP, (B) AFB1, (C) DEX, (D) PB, and (E) cisplatin in human TK6 cells at 4 hr, 8 hr, and 24 hr. The TGx‐28.65 classifier was used to distinguish genotoxic from nongenotoxic compounds in the presence of rat liver S9. The heat map on the left hand side depicts the 28 reference compounds used to generate the TGx‐28.65 biomarker. The labels on the right hand side of the heat map are the GenBank accession numbers for the classifier genes contained within the predictive gene signature. The color scale indicates gene expression fold‐changes relative to control: up‐regulated genes are shown in red, down‐regulated genes are shown in green, and genes that are not regulated are shown in black. Predictions of genotoxicity and NSC classification probabilities using the TGx‐28.65 classifier for all four treatments are shown using red (genotoxic) and blue (nongenotoxic) bars above each heat map.

Figure 4

Comparison of the transcriptome profiles of AFB1, BaP, DEX, and PB in human TK6 cells to the training set of genotoxic and nongenotoxic chemicals used to develop the toxicogenomic classifier (TGx‐28.65). Hierarchical clustering of the expression levels of genes in this predictive genotoxicity signature is depicted.

The classification of the lowest concentration of BaP as nongenotoxic or unclassified suggests that the results are not simply due to the effects of the S9 alone. However, to examine this in more detail, we also studied the effects of BaP in the absence of S9. The highest concentration of BaP used above (10 μg/ml) was assayed in the absence of S9 against a matched VC (‐S9) control at the 4 hr, 8 hr, and 24 hr time points. In the absence of S9, the high concentration of BaP decreased RS to 75.4%, but did not result in increased MN frequencies relative to vehicle controls (data not shown). Applying the TGx‐28.65 biomarker revealed a nongenotoxic classification for BaP in the absence of metabolic activation at all three time points (Supporting Information Fig. S5). This confirms that BaP needs to be metabolically activated to cause the signature and that the signature is not the result of the parent compound.

In addition, to ensure that the genotoxic signature was not produced simply because of the MAS, we compared VC (−S9) to VC (+S9). There were no differences in relative survival or MN frequencies when these controls were compared (data not shown). The expression profile for this contrast was classified as nongenotoxic (data not shown). These results suggest that the genotoxic signature of BaP is due to the formation of genotoxic reactive metabolites in TK6 cells produced by enzymes within the S9 and not due to nontarget S9 effects on gene expression patterns. Please also see work on the nongenotoxic chemicals dexamethasone and phenobarbital below for further confirmation of this finding.

Aflatoxin B1

AFB1 concentrations ranging from 0.025 to 0.2 μM, in addition to positive, negative, and vehicle controls, were tested in a range‐finding experiment (data not shown). The concentrations selected for the subsequent microarray analysis were: 0.025, 0.075, and 0.1 μM. A decline in RS (70.3%) and a marginal induction of MN (2.0‐fold) was observed at the low concentration in the definitive study, indicating that this concentration has some degree of cytotoxicity and is weakly genotoxic (Fig. 5). There was a significant induction of ATF3 and CDKN1A at 8 hr and 24 hr (Fig. 6), but no effect on the RT‐qPCR panel of genes at the 4 hr time point at the low concentration. The mid concentration exhibited a significant decline to 48.1% RS and an 8.1‐fold induction of MN (Fig. 5). The highest concentration selected for the definitive study caused a significant decline in RS to 40.2% and a 13.7‐fold increase in MN frequency (Fig. 5). Both the mid and the high concentrations caused significant up‐regulation of ATF3 and GADD45A at 8 hr and 24 hr, in addition to induction of CDKN1A at all three time points (Fig. 6). The heat map of AFB1‐treated TK6 cells showed concentration‐ and time‐dependent effects using the TGx‐28.65 classifier (Fig. 3B). NSC classification probabilities using the 65‐gene classifier for AFB1 treatment are shown by the prediction bar in Figure 3B. The toxicogenomic signature classified AFB1 as a genotoxic compound at the mid and high concentrations at all time points (Fig. 3B). In the absence of S9, the high concentration of AFB1 did not result in a decline in RS, nor an induction of MN relative to vehicle controls (data not shown). Furthermore, the high concentration of AFB1 was classified as nongenotoxic in the absence of S9 at 4 hr, 8 hr, and 24 hr (Supporting Information Fig. S5). Finally, AFB1 was identified as a genotoxic compound by comparison of its transcriptome profile in human TK6 cells to the TGx‐28.65 classifier by hierarchical clustering with the set of genotoxic and nongenotoxic reference chemicals for all concentrations and time points (Fig. 4).

Figure 5

Cytotoxicity and genotoxicity measurements in human TK6 cells following exposure to AFB1 using the In Vitro Microflow kit (Litron Laboratories). Percent relative survival (▪), percentage of apoptotic/necrotic cells (●), and percentage of MN (▴) are depicted following 24 hr time point (4 hr exposure + 20 hr recovery; n = 3). *P < 0.001 compared with the vehicle control. Error bars depict the standard error. % RS is a measure of the relative survival, % apoptosis is a measure of apoptotic cells, and % MN is a measure of micronucleus induction.

Figure 6

Normalized gene expression of targeted genes, including ATF3, CDKN1A, and GADD45A, in human TK6 cells by TaqMan RT‐qPCR following AFB1 exposures at 4 hr, 8 hr, and 24 hr. Normalized expression was calculated using human GAPDH as a reference gene (n = 3). Vehicle control +S9 (VC (+S9)) served as the negative control for AFB1 treatment. *P < 0.05, **P < 0.01, and ***P < 0.001 compared with the vehicle control. Error bars depict standard error.

Dexamethasone

Range finder studies tested concentrations from 7.8125 to 250 μg/ml, in addition to positive, negative, and vehicle controls (data not shown). Since DEX is nongenotoxic, we did not expect to see the induction of MN; however, we were aiming to achieve 55 ± 5% cytotoxicity as the top concentration. At 250 µg/ml (0.63 mM), a decline in RS to 82.4% was observed. Thus, we based our selection of maximum concentration for mammalian cell assays on the ICH Guidelines [ICH, 1996, 2011]. The old ICH Guideline criteria for the selection of top concentration for mammalian cell assays stated 5 mg/ml or 10 mM, whichever is lower; the revised guidelines have reduced the top concentration by 10‐fold to 500 µg/ml or 1 mM, when solubility in the culture medium and/or cytotoxicity are not limiting factors [ICH, 1996, 2011]. We prepared both of the recommended top concentrations: 1 mM (392.46 μg/ml) and 10 mM (3924.6 µg/ml), in addition to 5 mM and 7.5 mM concentrations. The 1 mM, 5 mM, and 7.5 mM concentrations resulted in a slight precipitate, but not enough to impede cellular exposures. The 10 mM concentration heavily precipitated out of solution and was not suitable for the subsequent cellular exposures. Thus, the 0.63 mM, 1 mM, and 7.5 mM concentrations of DEX were selected for gene expression analysis. These concentrations caused a concentration‐dependent decline in relative survival to 82.4%, 71.3%, and 66.4%, respectively (Fig. 7). There was minimal induction of MN at all concentrations; MN frequencies were 1.3, 1.2, and 1.1‐fold for the 0.63 mM, 1 mM, and 7.5 mM concentrations relative to controls, respectively (Fig. 7). ATF3, CDKN1A, and GADD45A were either unchanged or down‐regulated at 4 hr, 8 hr, and 24 hr following treatment of TK6 cells with DEX at all concentrations (Fig. 8). The heat map of DEX‐exposed TK6 cells revealed concentration‐ and time‐dependent effects using the TGx‐28.65 classifier and a very distinct profile from the positive controls (Fig. 3C). NSC classification probabilities using the classifier for DEX treatment revealed nongenotoxic classification at all concentrations and time points, which is shown by the prediction bar in Figure 3C. In this instance, only cisplatin, the positive control, was classified as genotoxic. In the absence of S9, the high concentration of DEX also classified as nongenotoxic at 4 hr, 8 hr, and 24 hr (Supporting Information Fig. S5). The genomic signature obtained for DEX‐treated TK6 cells classified this chemical as a nongenotoxic compound at all concentrations for all of the time points tested by comparison of its transcriptome profile to the reference chemicals used to generate the TGx‐28.65 classifier (Fig. 4).

Figure 7

Cytotoxicity and genotoxicity measurements in human TK6 cells following exposure to DEX using the In Vitro Microflow kit (Litron Laboratories). The experiments were performed separately for each concentration and thus are shown separately. Percent relative survival (▪), percentage of apoptotic/necrotic cells (●), and percentage of MN (▴) are depicted following 24 hr time point (4 hr exposure + 20 hr recovery; n = 3). *P < 0.001 compared with relevant vehicle control (i.e., VC1 (+S9) vs. 0.63 mM DEX, VC2 (+S9) vs. 1 mM DEX, and VC3 (+S9) vs. 7.5 mM DEX). Error bars depict the standard error. % RS is a measure of the relative survival, % apoptosis is a measure of apoptotic cells and % MN is a measure of micronucleus induction.

Figure 8

Normalized gene expression of targeted genes, including ATF3, CDKN1A, and GADD45A, in human TK6 cells by TaqMan RT‐qPCR following DEX exposures at 4 hr, 8 hr, and 24 hr. Normalized expression was calculated using human GAPDH as a reference gene (n = 3). Relevant vehicle control +S9 (VC (+S9)) served as the negative control for DEX treatments (i.e., VC1 (+S9) for 0.63 mM, VC2 (+S9) for 1 mM, and VC3 (+S9) for 7.5 mM). *P < 0.05, **P < 0.01, and ***P < 0.001 compared with the appropriate vehicle control. Error bars depict standard error.

Phenobarbital

Following a range finder study (data not shown), the concentrations selected for the definitive study were: 1, 3.2, and 10 mM. Relative survival, apoptosis and MN induction are shown in Figure 9. We remind the reader, that PB was tested at the 8 hr and 24 hr time points exclusively. At the low concentration, there was no indication of cytotoxicity or genotoxicity (Fig. 9). The mid concentration resulted in a decline in RS to 85.1% without induction of MN, while the high concentration caused a decline in RS to 53.8% without significant induction of MN (Fig. 9). The expression of ATF3 and CDKN1A increased by more than twofold for the mid and high concentrations compared with the vehicle controls at 8 hr, but were otherwise unchanged (Fig. 10). The heat map of PB‐exposed TK6 cells revealed concentration‐ and time‐dependent effects using the TGx‐28.65 classifier and a very distinct profile from the positive controls (Fig. 3D). NSC classification probabilities using the classifier for PB treatment are shown by means of the prediction bar in Figure 3D. PB was classified as nongenotoxic at all concentrations and time points, except for the high concentration at the 24 hr time point, in the presence of S9 metabolic activation. The high concentration of PB was incorrectly classified as genotoxic at 24 hr, indicating confounding effects of cytotoxicity at late time points that may lead to misclassification. Hierarchical cluster analysis would lead us to classify all other concentrations of PB as a nongenotoxic compound by comparison of its transcriptome profile in human TK6 cells using the TGx‐28.65 classifier when aligned against the database of genotoxic and nongenotoxic reference chemicals (Fig. 4).

Figure 9

Cytotoxicity and genotoxicity measurements in human TK6 cells (n = 3) following exposure to PB using the In Vitro Microflow kit (Litron Laboratories). Percent relative survival (▪), percentage of apoptotic/necrotic cells (●) and percentage of MN (▴) are depicted at the 24 hr time point (4 hr exposure + 20 hr recovery; n = 3). *P < 0.001 compared with the vehicle control. Error bars depict the standard error. % RS is a measure of the relative survival, % apoptosis is a measure of apoptotic cells, and % MN is a measure of micronucleus induction.

Figure 10

Normalized gene expression of targeted genes, including ATF3, CDKN1A, and GADD45A, in human TK6 cells by TaqMan RT‐qPCR following PB exposures at 8 hr and 24 hr. Normalized expression was calculated using human GAPDH as a reference gene (n = 3). Vehicle control +S9 (VC (+S9)) served as the negative control for PB treatment. Error bars depict standard error.

Cisplatin

Cisplatin is a direct‐acting DNA‐alkylating agent and was one of the 28 training chemicals used in the development of the TGx‐28.65 genotoxicity classifier. As such, we included this direct‐acting positive control in our experiments to ensure the genotoxicity classifier was working in the absence of S9 while we were testing the performance of the classifier in the presence of a MAS. Microarray experiments were conducted using the same concentration of cisplatin used in the development of the classifier (80 μM or 24 μg/ml). Cisplatin caused a significant decline in relative survival (43.3%) and a significant induction in MN (17.4‐fold increase compared with the VC (−S9) controls; data not shown). Significant induction of ATF3, CDKN1A, and GADD45A at 4 hr, 8 hr, and 24 hr was noted (Supporting Information Fig. S6). The heat map of cisplatin‐treated TK6 cells showed time‐dependent effects using the TGx‐28.65 classifier (Fig. 3E). NSC classification probabilities using the 65‐gene classifier for cisplatin treatment are shown by the prediction bar in Figure 3. The toxicogenomic signature classified cisplatin as a genotoxic compound for all time points, in both the presence and absence of S9 (Fig. 3E, Supporting Information Fig. S5). Hierarchical cluster analysis would lead us to classify cisplatin as a genotoxic compound by comparison of its transcriptome profile in human TK6 cells using the TGx‐28.65 classifier when aligned against the database of genotoxic and nongenotoxic reference chemicals (Fig. 4).

External Validation of the TGx‐28.65 Biomarker in Human HepaRG Cells

As a means to externally validate the TGx‐28.65 genotoxicity classifier, a publically available microarray dataset was accessed to test the performance of the genomic biomarker in another human cell line, hepatoma‐derived HepaRG cells. HepaRG cells have inherent metabolic capabilities, thus eliminating the requirement for the addition of an external metabolic activation system [Doktorova et al., 2013]. In a study by Doktorova et al. [2013], HepaRG cells were exposed to 15 well‐characterized chemicals belonging to three different classes: (i) genotoxic carcinogens, (ii) nongenotoxic carcinogens; (iii) noncarcinogens for 72 hr. The TGx‐28.65 classifier correctly classified all 15 chemicals using 53 of the 65 transcripts that were present on both the Affymetrix arrays used in the aforementioned study and the Agilent arrays used in the current study (Fig. 11).

Figure 11

Analysis of an externally derived dataset of human HepaRG transcriptomic profiles using the TGx‐28.65 classifier to predict genotoxicity following exposure to 15 compounds belonging to three toxic classes: (i) genotoxic carcinogens [AFB1; 4‐(methylnitrosamino)‐1‐(3‐pyridyl)‐1‐butanone (NNK); 2‐nitrofluorene (2NF); BaP; cyclophosphamide (CYCLO)], (ii) nongenotoxic carcinogens [methapyrilene hydrochloride (MPH); piperonylbutoxide (PIPB), Wy‐14643 (WYE), phenobarbital sodium (SPB), 12‐O‐tetradecanoylphorbol‐13‐acetate (TPA)]; (iii) non‐carcinogens [nifedipine (NIF); clonidine (CND); d‐mannitol (MAN); tolbutamide (TOL); diclofenac sodium (SDF)] [Doktorova et al., 2013].

DISCUSSION

Genotoxicity testing is a critical component of the chemical testing paradigm. The current test battery includes in vitro assays in bacterial and mammalian cells to test a substance's inherent ability to compromise the integrity of genetic material [Muller et al., 1999; FDA, 2007 (revised); ICH, 2011]. Increasing emphasis is being placed on obtaining mechanistic information to better inform human health risk assessment [Thybaud et al., 2007]. In our companion article, Li et al. (2015) developed a mechanism‐based genomics biomarker (TGx‐28.65) that is enriched in P53‐regulated DNA damage response genes. The TGx‐28.65 classifier was used to test three chemicals to demonstrate its ability to classify a genotoxic compound, a nongenotoxic compound and an irrelevant positive (i.e., positive genotoxicity results in vitro and negative genotoxicity results in vivo) using isopropyl methanesulfonate (iPMS), 3‐nitropropionic acid (3‐NP), and tri‐methylxanthine (caffeine), respectively. The present study demonstrates the ability of the TGx‐28.65 classifier to predict genotoxicity in the presence of rat liver S9 in human TK6 cells (to expand its application potential), as well as predict genotoxicity in a more metabolically competent human liver hepatocyte cell line (HepaRG) using a different microarray technology and experimental design. Thus, in addition to confirming the utility of the TGx‐28.65 biomarker in the presence of S9 and in a different cell line, the work also provides further validation for the TGx‐28.65 classifier overall by demonstrating its efficacy in an independent data set produced in two different laboratories using different technologies.

It has already been demonstrated that transcriptomic‐based analysis can correctly classify chemical genotoxicity in vitro in HepG2 cells in both the presence of β‐naphthoflavone‐/phenobarbital‐induced rat liver S9 [Boehme et al., 2011] and in the absence of S9 [Magkoufopoulou et al., 2012], as well as in mouse fibroblast cells [Rohrbeck et al., 2010] and in mouse lymphoma cells, both with and without metabolic activation (i.e., S9) [Dickinson et al., 2004; Hu et al., 2004]. The ability to differentiate genotoxic mechanisms has also been demonstrated in human TK6 cells with direct acting chemicals [Ellinger‐Ziegelbauer et al., 2009]. Thus, proof of principle has been published, but internationally coordinated validation exercises have not been conducted for any of these genomic biomarkers. Moreover, characterization of the concentration‐response of the signature anchored against traditional measures (e.g., MN frequency) has not been explored in detail.

In our study, TK6 cells were treated with increasing concentrations of BaP, AFB1, DEX, and PB in the presence of rat liver S9, and expression profiles were derived at 4 hr, 8 hr, and 24 hr (note that the 4 hr time point was not done for PB). We anchored the microarray gene expression results to relevant cytotoxicity (RS and apoptosis) and genotoxicity (MN) data, as well as to targeted expression of three key stress response genes (ATF3, CDKN1A, and GADD45A) by TaqMan RT‐qPCR (as recommended in [Ellinger‐Ziegelbauer et al., 2009]) in order to be consistent with the dose selection strategy employed by Li et al. (2015; companion article). NSC classification probabilities and hierarchical clustering for each chemical were performed using the TGx‐28.65 biomarker to predict whether BaP, AFB1, DEX, and PB are genotoxic or nongenotoxic. The BaP and AFB1 mid and high concentrations were potent inducers of MN, apoptosis, and cytotoxicity, and all were classified with >90% probability as genotoxic, as expected. Although DEX treatment caused robust changes in gene expression (False Discovery Rate corrected P value <0.05 and a fold change ≥±2) at both 0.63 mM and 1 mM concentrations at 4 hr (116 genes and 254 genes deregulated, respectively), at 8 hr (185 and 359, respectively) and at 24 hr (246 and 291, respectively), these gene expression changes were not associated with a genotoxic signature (data not shown). All concentrations of DEX tested (including 7.5 mM tested at 24 hr, the highest soluble concentration sampled at a time point exhibiting cytotoxicity) produced signatures that were greater than 90% predictive of being nongenotoxic at all times points. PB was classified as nongenotoxic at all concentrations and time points, except the highly cytotoxic concentration at 24 hr. This may result from the confounding effects of cytotoxicity that influence expression profiles at 24 hr [Ellinger‐Ziegelbauer et al., 2009]. We also note that while the S9 MAS may cause mild cytotoxicity [Boehme et al., 2011] and induce changes in gene expression, these changes do not mask the ability to detect a nongenotoxic signature or produce a genotoxic classification on their own. Indeed, the TGx‐28.65 signature produced by comparing vehicle control (+S9) to vehicle control (–S9) was clearly nongenotoxic. Thus, the misclassification of PB is not due to interference by S9. Overall, based on the profiles obtained with these test chemicals, the HESI TGx‐28.65 classifier appears to be effective at classifying genotoxic and nongenotoxic chemicals in the presence of this particular S9 MAS at 8 hr.

In addition to demonstrating that the classifier works in the presence of metabolic activation, the data suggest that the classifier is also partially quantitative in nature; the low concentrations of BaP and AFB1, which were either nongenotoxic or very weakly genotoxic, yielded equivocal results (i.e., no clear evidence to support genotoxicity or nongenotoxicity, with the exception of AFB1 at 4 hr, which classified as genotoxic). However, hierarchical cluster analysis with the training set compounds using the TGx‐28.65 biomarker revealed clustering of these low doses with known genotoxic agents (all time points except for the low concentration of BaP at 4 hr, which clustered with the nongenotoxic chemicals). Therefore, the predictive accuracy of the marker may provide important information relating to probability of genotoxicity even at low concentrations. We are presently extending this analysis to an additional 10 genotoxic and nongenotoxic compounds in the presence of different S9 MASs to confirm this finding in TK6 cells.

As microarray analysis gives us a fixed snapshot of gene expression changes at a specific point in time, it is critical to choose the appropriate time point in moving forward with further testing and validation of this genotoxicity classifier in the presence of metabolic activation. Our initial work focused on three time points (4 hr, 8 hr, and 24 hr). The 4 hr time point was chosen as the HESI database of direct‐acting reference chemicals was constructed using a single 4 hr time point [Goodsaid et al., 2010; Li et al., 2015]. We also included an 8 hr and a 24 hr time point to allow for adequate metabolic activation of these chemicals to their reactive metabolites. It has been suggested that early time points (i.e., 4 hr and 7 hr) may be more likely to generate gene expression data relevant to the mechanisms of genotoxic damage in the absence of appreciable cytotoxicity, which lends additional support to the inclusion of the 4 hr and 8 hr time points [Ellinger‐Ziegelbauer et al., 2009]. In an additional study conducted in HepG2 cells, in which cells were treated with 34 compounds for 12 hr, 24 hr, and 48 hr with the intent of developing and testing a transcriptomics‐based in vitro assay for predicting chemical genotoxicity in vivo, it was determined that the 24 hr exposure duration led to the most accurate predictions of genotoxicity [Magkoufopoulou et al., 2012]. Taken together, we decided to include the 4 hr, 8 hr, and 24 hr time points in the initial development and optimization of this genotoxicity classifier in the presence of metabolic activation. We demonstrate that PB was misclassified at the highest concentration at the 24 hr time point. The findings support that misclassification may occur at the 24 hr time point under conditions of high levels of cytotoxicity as has been suggested previously [Amundson et al., 2005; Ellinger‐Ziegelbauer et al., 2009]. Given that the TGx‐28.65 biomarker correctly classified all chemicals at the 8 hr time point (i.e., BaP and AFB1 were classified as genotoxic and DEX and PB were classified as nongenotoxic), we propose that this is the optimal sampling time for chemicals requiring metabolic activation. However, at this time we recommend the use of multiple time points for testing in the presence of S9 until the TGx‐28.65 genotoxicity classifier is further validated and optimized. We also advise that all chemical testing applying the TGx‐28.65 genomic biomarker be done both in the presence and absence of a MAS, as is standard practice for in vitro genetic toxicology assays.

The TGx‐28.65 genomic biomarker also correctly classified 15 chemicals belonging to three chemical classes (genotoxic carcinogens, nongenotoxic carcinogens, and noncarcinogens) in a different human cell line (HepaRG cells) using an Affymetrix microarray platform. Thus, the biomarker appears to provide a robust measure of genotoxic potential because it performs well in a different experimental model using a different microarray technology. This is very promising given the rapidly evolving methodologies used for gene expression profiling and the increasingly sophisticated cell culture models being developed to meet challenges associated with reduced animal use. The data suggest that the biomarker may be of potential value to alternative models and approaches, such as integration with high‐throughput screening.

In summary, we have demonstrated the utility of the TGx‐28.65 genomic biomarker in classifying genotoxic and nongenotoxic compounds in human TK6 cells in the presence of metabolic activation for BaP, AFB1, DEX, and PB. This model, which was anchored to relevant cytotoxicity, genotoxicity, and targeted gene expression data to ensure appropriate concentration selection, is of moderate throughput, is relatively cost‐effective, and was generated using a validated, commercial microarray platform. The use of the TGx‐28.65 classifier is meant to complement the current genotoxicity test battery to offer mechanistic information on genomic responses associated with genotoxicity. However, formal and rigorous genomic biomarker validation programs may enable qualified gene classifiers to eventually replace existing costly, time‐consuming and animal‐intensive tests with more insightful and biologically‐relevant tools to better inform the human health risk assessment process. Formal validation programs, such as the FDA Biomarker Qualification Program, are effective facilitators of new biomarker development to expedite their eventual use in toxicity testing for human health risk assessment [Goodsaid et al., 2010]. These aforementioned activities are currently ongoing for the HESI TGx‐28.65 classifier.

Acknowledgments

The authors thank Amanda Green, Kim Shepard, and Phillip Garibaldi for their valued technical support in conducting all of the cellular exposures, in addition to the cytotoxicity and genotoxicity testing, at Integrated Laboratory Systems, Inc. The authors also thank Alexandra Long and Dr. Errol Thomson for their insightful comments during the revision of this article.

AUTHOR CONTRIBUTIONS

A.F., J.A., and C.Y. designed the study and made important intellectual contributions to the manuscript. C.Y. obtained funding to support the project. L.R. and C.S. supervised the cellular exposures and flow cytometry‐based methods. J.B. and I.M. conducted all of the microarray and RT‐qPCR experiments. A.W. and D.H. conducted all the statistical analysis and prepared some of the figures. J.B. prepared the manuscript with important intellectual input from C.Y. All authors approved the final manuscript. C.Y., A.W., I.M., and J.B. had complete access to the study data.

Additional Supporting Information may be found in the online version of this article

Supporting Information Figure S1

Click here for additional data file.

Supporting Information Figure S2

Click here for additional data file.

Supporting Information Figure S3

Click here for additional data file.

Supporting Information Figure S4

Click here for additional data file.

Supporting Information Figure S5

Click here for additional data file.

Supporting Information Figure S6

Click here for additional data file.

Supporting Information Table S1

Click here for additional data file.

Supporting Information Table S2

Click here for additional data file.