Cortical neuron culture

Primary mouse cortical neuron cultures were prepared as previously described from E14.5 pregnant C57BL/6J (Cat. #000664, Jackson) dams crossed to CAST/EiJ (Cat. #000928, Jackson) males18. Hybrid cultures are thought to better model genetic variation associated with human populations58. Dissociated cells were placed in multiwell plates coated with poly-D-lysine (0.1 mg ml−1) in Neurobasal medium (Life Technologies) containing 5% fetal bovine serum (Gibco), B27 (17504-044, Invitrogen), Antibiotic-Antimycotic (15240-062, Invitrogen) and GlutaMAX (35050-061, Invitrogen). At DIV 3, a half medium change was performed with feeding medium identical to the plating medium except that we omitted fetal bovine serum and included 4.84 μg ml−1 uridine 5′-triphosphate (U6625, Sigma-Aldrich) and 2.46 μg ml−1 5-fluoro-2′-deoxyuridine (F0503, Sigma-Aldrich) to inhibit mitosis in dividing cells.

Chemicals

The chemical library was donated by the EPA ToxCast Program. For verification and replication experiments, the following chemicals were purchased from Sigma-Aldrich: DL-α tocopherol acetate (T3376), DL-sulforophane (S4441), vincristine sulfate salt (V8879), oxyfluorfen (35031), rotenone (45656), fenamidone (33965), pyraclostrobin (33696), trifloxystrobin (46447), myxothiazol (T5580), pyridaben (46047), azoxystrobin (31697), fluoxastrobin (33797), fenpyroximate (31684) and kresoxim-methyl (37899). Famoxadone was purchased from Chem Service, Inc (N-11943). Topotecan hydrochloride was purchased from Tocris (4562). Paclitaxel was purchased from Fisher Scientific (AC32842). All chemical stocks were prepared in DMSO unless otherwise noted. Vehicle samples were prepared with an equivalent DMSO concentration of ≤0.5% in feeding medium.

Live/dead assay

Cultures were grown on poly-D-lysine coated 96-well plates at a density of 5 × 104 cells per well. Culture plates were screened for quality by qualitative assessment of density and by confirming low percentage (<10%) of dead cells in control wells as established by the method described below. Chemicals were applied in quadruplicate by performing a half medium change with a 2 × concentration in prewarmed feeding medium on culture plates meeting these quality criteria, at DIV 7. After 24 h, another half medium change was performed with medium containing NucBlue (R37605, Life Technologies) and SYTOX Green (S7020, Life Technologies) dyes according to the manufacturer’s directions. Cells were then rinsed twice with PBS (137 mM NaCl, 10 mM Na 2 HPO 4 , 1.8 mM KH 2 PO 4 ) and fixed with 4% paraformaldehyde in 0.1 M phosphate buffer for 10 min at room temperature. After rinsing an additional two times with PBS, the cells were imaged on a Nikon Eclipse Ti epifluorescent microscope at × 10 magnification with ultraviolet excitation to reveal the NucBlue dye marking all cells, and with fluorescein isothiocyanate filter set to identify dead SYTOX-positive cells. Channel-merged colour images were imported into ImageJ software (NIH) and converted into RGB stacks. A manual threshold was performed on the NucBlue channel to yield an image amenable to binary particle analysis. The threshold set for the NucBlue channel was utilized for the SYTOX channel to reduce any potential bias. Particle analysis was performed on each channel for each image to obtain the total number of nuclei and the total number of nuclei from dead cells, respectively. The percentage of dead cells was calculated and averaged across four replicates. All chemicals from the library were initially tested at 10 μM. A concentration causing 10% or greater cell death relative to vehicle was considered toxic. Lower doses were tested by reducing an order of magnitude until a non-toxic dose was identified.

RNA-seq

Cultures were treated with the non-cytotoxic dose of each chemical, by exchanging half of the original medium with chemicals at 2 × concentration in prewarmed feeding medium, on DIV 7 for 24 h in 12-well plates at a density of 5 × 105 cells per well. RNA was isolated using RNeasy plus mini kit (Cat. #74134, Qiagen). RNA yield and quality were determined using a Nanodrop 1000 Spectrophotometer (Thermo Scientific). Samples were further assessed for quality using either an Agilent Bioanalyzer 2100 or TapeStation 2200 to obtain a RIN. RIN values exceeding 7 were used for sequencing. RNA samples were used to generate and barcode complementary DNA libraries using the TruSeq RNA Library Preparation Kit at the UNC High Throughput Sequencing Facility. Pools of 24 multiplexed samples were sequenced per lane on an Illumina HiSeq 2500 using 50-bp paired-end reads.

RNA-seq data processing

RNA-seq reads were filtered using TagDust and aligned to the reference mouse genome (mm9) with TopHat using default parameters59,60. Reads aligning to ribosomal RNA genes were removed. Transcript abundance was estimated by computing Reads per kilobase per million mapped reads (RPKM) using RefSeq gene models aggregated by gene symbol61. For differential expression analyses, raw counts over RefSeq exons were used, and then were compared across samples using DESeq62.

Hierarchical clustering and cluster membership

RPKMs for all RefSeq genes and all samples were filtered such that a given gene had >90% of samples with RPKM >0 and gene length exceeded 500 bp. RPKMs were then log 10 -transformed, quantile-normalized, median-centred, and the effects due to culture or sequencing batch were removed using a mixed model analysis of variance. This batch correction step was performed on individual experiments rather than chemical–vehicle ratios to reduce the effects due to variations in culture composition as well as batch-to-batch variations at the level of sequencing. Genes were then median-centred again, and filtered such that their s.d. exceeded 0.1. Chemical replicates were then combined using the median. Classes of chemicals were then discovered in two ways. First, we performed hierarchical clustering using average linkage and Pearson correlation distance measures. The boundaries of the discovered classes from clustering were set by computing pairwise Spearman correlations for all samples. A given cluster had to have at least three members with the minimum pairwise Spearman correlation coefficient exceeding 0.2, though cluster members often exceeded a correlation of 0.6 (Supplementary Fig. 5c,d). Seven chemical clusters were discovered, however one was removed due to failing to have a positive silhouette width, a metric of cluster robustness63. These definitions enabled the three positive control topoisomerase inhibitors to form a cohesive group of chemicals with known structural and functional similarity (cluster 5). None of the culture or sequencing batches contributed to the formation of any one of the final six clusters (Supplementary Fig. 4).

Hierarchical clustering of EPA toxicological assays

Concentration at 50% maximum activity (AC 50 ) values for all assays and all ToxCast Phase I chemicals were retrieved from the EPA Dashboard. We retained only those assays with data available for 99% of the chemicals and required that a given assay must have at least five active chemicals (AC 50 ≤10 μM); 199 assays were retained. We then transformed the data values as described in Sipes et al.64 and performed hierarchical clustering to find the chemical groups with concordant activity on certain assays.

Pathway analysis

To detect differentially regulated pathways, chemicals within each of the six clusters were separately compared with all other chemicals tested using GSA65, supplying a modified gene pathway file that contained all MSigDB C2 annotations as well as gene sets from published gene expression studies in human diseases and mouse models. Nearly all added gene sets were derived directly from the associated publication. Raw data derived from Durrenberger et al.66 were analysed by applying quantile normalization, and differential expression was detected using SAM using 500 permutations67. Some brain disease gene expression studies had no significantly differentially expressed genes (q<0.05), and hence could not be included in our analysis. For ‘BLALOCK_ALZHEIMERS_UP’, a gene set included within MSigDB C2, we removed gene symbols that were not included in the RefSeq annotation. The Mod1 pathway (Gupta et al.10) was filtered to include genes with a KM1 score of 0.4 or greater, to reduce the number of genes in this pathway below 1,500 so it could be analysed by GSA. A given pathway (size between 10 and 1,500 genes) was considered if it achieved statistical significance (false discovery rate<0.1) for at least one cluster compared with all other chemicals (500 permutations). Expression of a single pathway was then summarized by taking the median expression value across all genes in that pathway for a given chemical. Once the entire data matrix was assembled, pathways were median-centred and hierarchically clustered, while keeping chemical ordering consistent with Fig. 2. Disease-relevant ontologies were then extracted and plotted separately; colour was applied based on the pathway enrichment score (blue indicates pathway downregulation and red indicates pathway upregulation), but only for the pathways achieving statistical significance (false discovery rate<0.1).

Quantitative real-time PCR

Total RNA was extracted and purified using the RNeasy plus mini kit (Qiagen) following the manufacturer’s instructions. RNA was quantified using Nanodrop 1000 (Thermo Scientific) and reverse transcribed using iScript Reverse Transcription Supermix (Bio-Rad). Quantitative real-time PCR analysis was performed in technical duplicates using SYBR Green PCR master mix (Applied Biosystems) in a 7500 Real-Time PCR instrument (Applied Biosystems). An amount of 10 ng of complementary DNA was used in each reaction. Fold change of expression of the target RNA in each treatment condition relative to the untreated sample were calculated by raising 2 to the power of −ΔΔCt. ΔΔCt was calculated by subtracting the ΔCt of the untreated sample to the ΔCt of each treatment condition. The ΔCt of each sample was calculated by subtracting the average Ct of Gapdh to the average Ct of the target RNA. The following primers were used at a final concentration of 0.5 μM:

Rest (F: 5′-GTGCGAACTCACACAGGAGA-3′, R: 5′-AAGAGGTTTAGGCCCGTTGT-3′); Gsta4 (F: 5′-CGGCTGGAGTGGAGTTTGAG-3′, R: 5′-CCAAGGGTACTTGGCCGAAA-3′); Gstm1 (F: 5′-CCGTGCAGACATTGTGGAGA-3′, R: 5′-CTGCTTCTCAAAGTCAGGGTTG-3′); Hmox1 (F: 5′-AGGCTTTAAGCTGGTGATGGC-3′, R: 5′-GGGGCATAGACTGGGTTCTG-3′); and Gapdh (F: 5′-TATGACTCCACTCACGGCAAAT-3′, R: 5′-GGGTCTCGCTCCTGGAAGAT-3′).

Fluorescence Immunohistochemistry

Mouse cortical neuron cultures were prepared identically as for live/dead and RNA-Seq experiments but at an equivalent density of 2.5 × 105 cells per well in a 24-well plate containing precoated glass coverslips. At DIV 7, cells were treated with fenamidone at a concentration of 10 μM for 24 h. Coverslip-adherent cells were rinsed, paraformaldehyde fixed, blocked with 10% donkey serum and probed with the following cell-type-specific primary antibodies: guinea pig anti-NeuN (1:400, ABN90P, Millipore), goat anti-Gfap (1:1,200, ab53554, Abcam) and rabbit anti-Iba1 (1:400, 019-19741, Wako Pure Chemical Industries). Species-specific AlexaFluor-conjugated donkey secondary antibodies (1:200, Life Technologies) were then applied along with Nucblue (Life Technologies) to counterstain all nuclei. Coverslips were inverted into mounting substance (Fluorogel, 17985-10, Electron Microscopy Sciences) and imaged on a Zeiss LSM 710 laser scanning confocal microscope using a × 10 objective tiling scan of entire coverslips including three z planes to account for uneven coverslips. A researcher blind to the experimental conditions of images utilized maximal projected, spliced images to obtain particle counts for the total Nucblue-positive nuclei and the NeuN-positive neurons using particle analysis in ImageJ. Gfap- and Iba1-positive features containing clear Nucblue-positive nuclei were manually counted using the count function in Adobe Photoshop.

Mitochondrial superoxide detection and live-cell imaging

Cultures were prepared on commercial poly-D-lysine-coated glass coverslips (GG-12-pdl, neuVitro) in 24-well plates at a density of 1.5 × 105 cells per well. At DIV 7, cells were treated with chemical or an equivalent DMSO concentration in feeding medium for 2 h. Cells were then treated with MitoSOX Red (M36008, Life Technologies) mitochondrial superoxide indicator for 10 min at 37 °C. Medium containing the dye was aspirated, and the cells rinsed twice with 37 °C feeding medium and replaced with 37 °C artificial cerebral spinal fluid composed of 150 mM NaCl, 5 mM KCl, 1 mM MgCl 2 , 2 mM CaCl 2 , 10 mM HEPES and 10 mM dextrose (pH 7.3). Individual coverslips were transferred to a stage-top perfusion system mounted on an inverted Nikon Ti Eclipse microscope with constant flow of warmed artificial cerebral spinal fluid. Images (× 20) were collected using an Andora Clara charge-coupled device camera for brightfield and the Texas Red channel to ascertain the mitochondrial superoxide levels with the same exposure across all experiments. ImageJ software (NIH) was used to trace the soma of the cells in bright-field images as regions of interest (ROIs) using the polygon selection tool. ROIs were then superimposed on the fluorescent image and the total cell fluorescence was calculated per ROI to account for the size of the ROI and normalize to background fluorescent levels, as described68. Fluorescence intensity was then normalized to average corrected vehicle intensity. At least two non-overlapping fields from each coverslip were collected for each well. ROI tracing was performed by a researcher blind to the experimental condition. For select conditions, the proportion of cells showing altered soma morphology was manually calculated per image by the same experimenter performing ROI tracing.

Tubulin polymerization assay

Chemicals were tested in a commercially available cell-free biochemical assay (BK011P, Cytoskeleton, Inc). Briefly, control and experimental solutions were prepared and loaded in duplicate into a warm (37 °C) 96-well plate. A solution of purified tubulin and a fluorescent reporter were added to the wells containing the warm substances and placed into a prewarmed (37 °C) Biotek Synergy HT microplate reader where fluorescence was measured every minute for 60 min. Equimolar doses of paclitaxel (included in kit) and vincristine (V8879, Sigma-Aldrich) were included as positive and negative modulators of tubulin polymerization, respectively.

Chemical usage and approval data

Usage data were acquired from United States Geological Survey. Data are reported as total kilograms applied for each chemical across all US counties sampled per year (2000–2012) and a linear trend line was fit beginning in year 2000 or with the registration year (if after 2000). Registration dates for specific pesticide products were obtained from the National Pesticide Information Retrieval System (NPIRS; http://ppis.ceris.purdue.edu/). The NPIRS is a collection of pesticide-related databases, applications and websites under the administration of the Center for Environmental and Regulatory Information Systems at Purdue University, West Lafayette, Indiana. NPIRS obtains product information on a weekly basis from the EPA Pesticide Product Information System website. This Pesticide Product Information System information, including registered use sites and pests listed on the EPA stamped approved label, is disseminated through the NPIRS member and public websites, and is provided for informational purposes only. To ascertain which foodstuffs had the greatest amount of these chemicals, we ranked the five food commodities in descending order of maximum residue level detected using 2008–2012 data from the USDA Pesticide Data Program and the FDA Pesticide Program Residue Monitoring Program.