Neonicotinoid Insecticides Alter the Gene Expression Profile of Neuron-Enriched Cultures from Neonatal Rat Cerebellum

Neonicotinoids are considered safe because of their low affinities to mammalian nicotinic acetylcholine receptors (nAChRs) relative to insect nAChRs. However, because of importance of nAChRs in mammalian brain development, there remains a need to establish the safety of chronic neonicotinoid exposures with regards to children’s health. Here we examined the effects of long-term (14 days) and low dose (1 μM) exposure of neuron-enriched cultures from neonatal rat cerebellum to nicotine and two neonicotinoids: acetamiprid and imidacloprid. Immunocytochemistry revealed no differences in the number or morphology of immature neurons or glial cells in any group versus untreated control cultures. However, a slight disturbance in Purkinje cell dendritic arborization was observed in the exposed cultures. Next we performed transcriptome analysis on total RNAs using microarrays, and identified significant differential expression (p < 0.05, q < 0.05, ≥1.5 fold) between control cultures versus nicotine-, acetamiprid-, or imidacloprid-exposed cultures in 34, 48, and 67 genes, respectively. Common to all exposed groups were nine genes essential for neurodevelopment, suggesting that chronic neonicotinoid exposure alters the transcriptome of the developing mammalian brain in a similar way to nicotine exposure. Our results highlight the need for further careful investigations into the effects of neonicotinoids in the developing mammalian brain.


Introduction
Neonicotinoids are used worldwide as pesticides, but their toxicity also affects beneficial insects, such as honey bees, and other invertebrate species [1,2]. Neonicotinoids are agonists of nicotinic acetylcholine receptors (nAChRs) and their effectiveness has been attributed to their high affinity for insect nAChRs. Furthermore, their low affinity for several subsets of mammalian nAChRs in binding assays suggests that neonicotinoids are relatively safe for mammals including humans [3,4].
However, there have been reports of severe or fatal cases of poisoning by the neonicotinoids imidacloprid (IMI), acetamiprid (ACE), and thiacloprid [5][6][7]. Marfo et al. reported the occurrence of subacute neonicotinoid intoxication through daily consumption of contaminated food in Japan [8]. They detected high levels of neonicotinoid metabolites in the urine of patients with symptoms typical of neo-nicotinoid intoxication, including headache, general fatigue, palpitation/chest pain, abdominal pain, muscle pain/weakness/spasm, cough, finger tremor, and short-term memory loss. Furthermore, an in vitro study reported that IMI and clothianidin agonize human α4β2 nAChR and disrupt receptor responses against the endogenous ligand acetylcholine (ACh) [9]. Moreover, some metabolites of the neonicotinoids have a high affinity for mammalian nAChRs, similar to that of nicotine (NIC) [3]. Together, these findings suggest that neonicotinoids might be a risk to human health, but there remains a paucity of data on this subject [10].
In rodent experiments, exposure to IMI in utero was shown to cause impairments in sensorimotor performance and overexpression of glial fibrillary acidic protein (GFAP) in the motor cortex and hippocampus of neonatal rats [11]. Furthermore, the neonicotinoids thiamethoxam and clothianidin induced dopamine release in the rat striatum via nAChRs [12] and thiamethoxam altered behavioral and biochemical processes related to the cholinergic systems in rat [13]. Gestational administration of clothianidin induces behavioral disorders in F1 mice [14] and exposure to clothianidin and other stress causes behavioral and reproductive abnormalities in male mice [15]. Exposure to IMI alters learning performance and related gene expression in infant and adult rats [16]. Recently, Sano et al. reported that exposure to low dose ACE in utero and through lactation induces abnormalities in socio-sexual and anxiety-related behaviors in male mice [17]. Additionally, an in vivo mouse study revealed that four types of neonicotinoids (ACE, IMI, thiacloprid, and nitenpyram) readily pass through the blood-brain barrier (BBB) and are detectable in the brain at 60%-70% of the serum concentration within 15 min of intraperitoneal administration [18].
We previously reported that excitatory Ca 2+ -influx was evoked by NIC, ACE, and IMI at concentrations greater than 1 µM in small neurons in neonatal rat cerebellar cultures expressing α3, α4, and α7 nAChR subunit mRNA [19]. The proportion of excited neurons, and peak excitatory Ca 2+ -influx induced by ACE and IMI were lower than those induced by NIC. However, ACE and IMI had greater effects on mammalian neurons than would be predicted based on the results of binding assay studies. Because three nAChR antagonists significantly inhibited ACE-and IMI-induced neuronal Ca 2+ -influx, it is likely that ACE and IMI have direct agonist activity at nAChRs in cerebellar neurons.
During mammalian brain development, expression of several nAChR subtypes is higher during the fetal stage than at the mature stage, and are important for synapse formations. In the human fetal cerebellum, α4, α7, β2, and β4 nAChRs are highly expressed [20]. Recently, several lines of evidence indicated that endogenous cholinergic signaling via nAChRs is important in determining the morphological and functional maturation of neural circuit formation [21]. Glutamatergic synapse formation is promoted by α7-containing nAChRs and affected by NIC exposure in hippocampal and cortical neurons [22]. Retinal β2 nAChRs are necessary for visual circuit formation [23] and prenatal NIC exposure alters the visual cortex system in baboons [24]. Even at a dose lower than that necessary to activate the receptor, NIC causes desensitization of nAChRs [25], which results in a disturbance of normal synapse formations at the developmental stage [26].
Although neonicotinoids have a lower affinity for mammalian nAChRs than for insect nAChRs, they have similar chemical characteristics to NIC, which is known to disrupt normal brain development. Therefore, their effects on human health especially in the developing brain must be carefully examined. In 2013, the European Food Safety Authority proposed that ACE and IMI had the potential to be developmentally neurotoxic [27]. However, the impact of neonicotinoids on mammalian nAChRs and brain development has not been sufficiently studied.
In the present study, we examined whether chronic neonicotinoid exposure affects the gene expression profile in the developing mammalian brain using microarrays of neuron-enriched cultures from neonatal rat cerebellum. Transcriptome profile analysis is a useful tool to demonstrate possible effects of neonicotinoids on multiple nAChR functions. Rodent neonatal cerebellar cultures are commonly used for studying developing neurons [28]. Maturation of the rat cerebellum takes about 30 days (from the 12-day embryo to postnatal day (P) 19) [29]. The cerebellum is used to model important neurodevelopmental processes in vitro, such as neuronal cell differentiation, neurite outgrowth, and synapse formation. Therefore, we used neuron-enriched cultures derived from cerebella from 1-day-old rats, with or without exposure to NIC, ACE, or IMI, and identified significant differential expression of several genes important for brain development.

Animals and Ethics Statement
Gestational Sprague-Dawley rats were purchased from Clea Japan, Inc. (Tokyo, Japan). All experimental procedures were performed in accordance with and approved by the Animal Use and Care Committee of the Tokyo Metropolitan Institute of Medical Science (permit no. 15003), and all animals were cared for and treated humanely in accordance with our institutional animal experimentation guidelines.

Neuron-Enriched Cultures of Neonatal Rat Cerebellum
The details of the culture methods have been described previously [19]. In brief, cerebella of P1 neonatal rats were digested with papain, and dissociated cells were suspended in a synthetic culture medium containing 1% fetal bovine serum. The cells were plated at a density of 1.2 × 10 6 cells/0.8 mL in a chamber slide well (4 cm 2 /well, Permanox, Nunc Lab-Tek, Sigma-Aldrich) that was precoated with 0.1 mg/mL poly-L-lysine (Sigma-Aldrich) and 10 µg/mL laminin (BD Biosciences, Franklin Lakes, NJ, USA). Four wells (4.8 × 10 6 cells) per group (control, NIC, ACE, and IMI) were used for each experiment. The remaining cerebellar cells were plated in chamber slide wells (0.7 cm 2 /well) at the same density as the 4 cm 2 wells, cultured in the same manner, and then used for immunocytochemistry.
We performed six independent culture experiments. For each experiment, cultures for all four groups were prepared from the same lot of cerebellar cells (from 12-15 neonates of one litter).

RNA Preparation
Total RNA was isolated from 16 DIV cerebellar cells using TRIzol reagent (Invitrogen, Thermo Fisher Scientific) and the Direct-zol RNA Miniprep Kit (Zymo Research, Irvine, CA, USA) following the manufacturers' recommendations. Quality of the total RNA was strictly controlled: by several parameters. The RNA extracts were analyzed by a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to confirm RNA quality with an A260/A280 ratio close to 2 and an A260/A230 ratio greater than 2. In total, 24 samples were tested, comprising six samples for each of the four treatments: control, NIC, ACE, and IMI.

Microarray
RNA quality was further assessed using a 2100 Bioanalyzer (Agilent Technologies, Valencia, CA, USA). cDNA and cRNA were synthesized with a Low Input Quick Amp Labeling Kit (Agilent Technologies) according to the manufacturer's instructions. The labeled cRNA samples were then hybridized to the Whole Rat Gene Expression ver.3 4x44k Microarray (G4847B; Agilent Technologies). The array slides were washed and scanned using a SureScan Microarray Scanner (Agilent Technologies). Microarray data were analyzed using Feature Extraction software v11.5.1.1 (Agilent Technologies).

Data Normalization and Quality Control
Data were analyzed using GeneSpring GX 13.1 software (Agilent Technologies). Raw signals were log transformed and normalized using the quantile normalization method. For each probe, the median of the log summarized values from all the samples was calculated and subtracted from each of the samples to get a transformed baseline. Quality control for each array was checked according to GeneSpring instructions. The final data set contained an expression signal above the 20th percentile in more than half the samples within range, so as not to exclude a particularly large change caused by some treatment. These microarray data have been submitted to the Gene Expression Omnibus and are accessible through accession number GSE 80656.

Differential Expression Analysis and Statistics
To determine the specific effects of NIC, ACE, and IMI, each treatment versus control was analyzed and the results were compared. Not all of the normalized array data followed a normal distribution, so we used the nonparametric paired Mann-Whitney U test; this reduced the variation within each culture lot (each from a single litter). Storey's bootstrapping estimation of the false discovery rate (FDR; q) was used to correct for multiple comparisons [30]. Differences were considered statistically significant for p and q values <0.05, and genes with a 1.5 (log 2 0.585)-fold change (FC) were selected for analysis. For several genes, the arrays contained multiple probes on the chip. In these cases, we limited each gene to corroborate direction and magnitude of change. We also excluded the compromised data and the probes without gene symbols or descriptions.
To determine α3, α4, and α7 nAChR mRNA expression in the cultures, we performed qRT-PCR analysis of the genes encoding these subunits (Chrna3, Chrna4, and Chrna7, respectively; primers purchased from Takara) using the same RNA samples.

Gene Ontology and Bioinformatics Analysis
To identify the molecular functions and diseases associated with the differentially expressed (DE) genes after NIC, ACE, or IMI exposure, we used the Protein Analysis Through Evolutionary Relationships (PANTHER) Classification System [35] and MetaCore analysis software (Thomson Reuters, New York, NY, USA). Overviews of DE gene ontology (GO) were examined using PANTHER. With regard to MetaCore, because ontologies typically contain many terms, some of them may turn out to be significant for a particular list and given p-value. The enrichment p values vary dependent on the size of the gene list and the selection of the background. The enrichment analysis FDR was controlled using the Benjamini and Hochberg method. Each DE gene list was analyzed by MetaCore enrichment analysis including GO processes (gene functions) and related diseases to identify genes involved in neurodevelopment. Next, all three gene lists were analyzed by comparison experiment tools, which can be used for comparing experimental data by analyzing their intersections.

Immunocytochemistry
At 16 DIV, cultured cerebellar cells were fixed with 4% paraformaldehyde in 0.1 M phosphate buffer for 20 min at room temperature. The following primary antibodies were used to detect specific antigens: mouse monoclonal anti-Tuj1 (Covance, Princeton, NJ, USA) as a marker for immature neurons, mouse monoclonal anti-calbindin D28k (Sigma-Aldrich) as a Purkinje cell marker, rabbit anti-GFAP (Dako, Glostrup, Denmark) as an astrocyte marker, mouse anti-oligodendrocyte marker O4 (Chemicon, Merck Millipore, Billerica, MA, USA), and mouse anti-CD11b OX42 IgM (Serotec, Bio-Rad), as a microglial marker. For double-labeling, the cells were incubated with primary antibodies against TuJ1, O4, and OX42 or GFAP, then Alexa 488-labeled anti-mouse IgG or IgM (Molecular Probes, Thermo Fisher Scientific), and Alexa 594 labeled anti-rabbit IgG secondary antibodies (Molecular Probes). The nucleus was stained with Hoechst 33342 (Molecular Probes). The cells were observed by confocal laser microscopy (Zeiss, LSM780, Oberkochen, Germany). Percentages of neural cells were calculated from the number of immunostained cells per Hoechst-positive-nucleus (n = approximately 1000) across three to four experiments, using MetaMorph imaging software (Molecular Device, Sunnyvale, CA, USA). To measure the dendritic arborization of Purkinje cells, the calbindin D28k-immunoreactive area was quantified from 10 cells per sample, using MetaMorph, in three independent experiments, as described previously [36] (see Figure 2 for representative results).

Morphological Characterics of 16 DIV Cerebellar Neuron-Enriched Cultures
To determine the effects of NIC, ACE, or IMI exposure on survival rate and morphological features of cells in the cerebellar cultures, we immunostained cerebellar cells for specific antibodies (Figure 1). We observed no notable morphological differences in Tuj1-positive neurons (large Purkinje cells, small granule cells, and other interneurons such as basket cells, Golgi cells, stellate cells), and glial cells (astrocytes, oligodendrocytes, and microglia) between the four groups in any of the six culture lots. However, staining for calbindin D28k to show the detailed morphology of Purkinje cells revealed a slight disturbance in the dendritic arborization of cells exposed to NIC, ACE, and IMI, and the reduction of dendritic area was 38%~40% in the treatment groups compared to a control ( Figure 2). In each of the four groups in each experiment, the majority (approximately 61%) of cells were Tuj1-positive neurons. Other cells included approximately 28% GFAP-positive astrocytes, 5% O4-positive oligodendrocytes, and 3% OX42-positive microglia (for details see supplemental Figure S1). Other neural cells, which were not identified by staining, were likely to be NG2 glial cells or neural stem cells [37,38]. Hoechst staining revealed approximately 1% abnormal nuclei (indicating cell death). . From the first to the fifth row, the same sample is shown. In the OX42 and O4 rows, red represents GFAP and blue is Hoechst stain. The photos of microglia and oligodendrocytes were obtained from areas of sparse cell density, because their morphology was easy to identify. Bars = 50 μm. No notable morphological differences were observed between the groups. . From the first to the fifth row, the same sample is shown. In the OX42 and O4 rows, red represents GFAP and blue is Hoechst stain. The photos of microglia and oligodendrocytes were obtained from areas of sparse cell density, because their morphology was easy to identify. Bars = 50 µm. No notable morphological differences were observed between the groups.

Figure 2.
Dendritic arborization of Purkinje cells in cerebellar cultures exposed to nicotine (NIC), acetamiprid (ACE), and imidacloprid (IMI) for 14 days (16DIV). A slight disturbance in the dendritic arborization exposed to NIC, ACE, and IMI was observed. Upper panel shows a representative cell from each group stained for anti-calbindin D28k. Bar = 50 μm. Graph shows MetaMorph quantification of calbindin D28k-reactive dendritic area (without cell soma) (n = 10 cells per group). Error bar represent standard errors. T-tests were conducted for each treatment versus control. * p < 0.05.

Differential Gene Expression after Exposure to NIC, ACE and IMI
For control versus NIC (CvN), 4550 of 23,011 filtered probes were identified as significant (p < 0.05, q < 0.05), and were used for differential expression analysis (fold change (FC) ≥ 1.5). The 92 probes that passed were checked for multiple probes, gene symbol, descriptions, and of these 34 DE genes were selected (Tables 1 and 2).
Next, we constructed Venn diagrams to visually assess the overlap and separation between the three DE gene lists, CvN (34 genes), CvA (48 genes), and CvI (67 genes). As shown in Figure 3 and Table 1, nine genes (four upregulated and five downregulated versus control) were common to three gene lists. Common to at least two lists were three genes between CvN and CvA, five genes between CvN and CvI, and four genes between CvA and CvI.
Other DE genes unique to CvN, CvA, or CvI are shown in Tables 2-4. The differential expression levels of the three groups are presented as heat maps and standard deviations that display expression differences between each of the treatments (NIC, ACE, and IMI) and control ( Figure 4, Table S2).

Figure 2.
Dendritic arborization of Purkinje cells in cerebellar cultures exposed to nicotine (NIC), acetamiprid (ACE), and imidacloprid (IMI) for 14 days (16DIV). A slight disturbance in the dendritic arborization exposed to NIC, ACE, and IMI was observed. Upper panel shows a representative cell from each group stained for anti-calbindin D28k. Bar = 50 µm. Graph shows MetaMorph quantification of calbindin D28k-reactive dendritic area (without cell soma) (n = 10 cells per group). Error bar represent standard errors. T-tests were conducted for each treatment versus control. * p < 0.05.

Differential Gene Expression after Exposure to NIC, ACE and IMI
For control versus NIC (CvN), 4550 of 23,011 filtered probes were identified as significant (p < 0.05, q < 0.05), and were used for differential expression analysis (fold change (FC) ≥ 1.5). The 92 probes that passed were checked for multiple probes, gene symbol, descriptions, and of these 34 DE genes were selected (Tables 1 and 2).
Next, we constructed Venn diagrams to visually assess the overlap and separation between the three DE gene lists, CvN (34 genes), CvA (48 genes), and CvI (67 genes). As shown in Figure 3 and Table 1, nine genes (four upregulated and five downregulated versus control) were common to three gene lists. Common to at least two lists were three genes between CvN and CvA, five genes between CvN and CvI, and four genes between CvA and CvI.
Other DE genes unique to CvN, CvA, or CvI are shown in Tables 2-4. The differential expression levels of the three groups are presented as heat maps and standard deviations that display expression differences between each of the treatments (NIC, ACE, and IMI) and control ( Figure 4, Table S2).

Quantitative qRT-PCR Validation
qRT-PCR was used to validate microarray data for 20 randomly selected DE genes and two control genes. Almost all genes showed consistent trends between qRT-PCR and microarray results ( Figure 5), with the Pearson correction coefficient of 0.80 between the qRT-PCR and microarray data. In Tada2b, log2 FC of PCR was noticeably lower than those of microarray, especially CvN. Of the 22 genes selected, changes in only three (Lyn, Pcdhgb7, and Dcdc2) were not confirmed by qRT-PCR. This may be because of differences in primer positions, resulting in transcript variants or degradation of mRNA.

Quantitative qRT-PCR Validation
qRT-PCR was used to validate microarray data for 20 randomly selected DE genes and two control genes. Almost all genes showed consistent trends between qRT-PCR and microarray results ( Figure 5), with the Pearson correction coefficient of 0.80 between the qRT-PCR and microarray data. In Tada2b, log 2 FC of PCR was noticeably lower than those of microarray, especially CvN. Of the 22 genes selected, changes in only three (Lyn, Pcdhgb7, and Dcdc2) were not confirmed by qRT-PCR. This may be because of differences in primer positions, resulting in transcript variants or degradation of mRNA. Figure 5. Confirmation of the microarray data using quantitative real-time PCR (qRT-PCR). Twenty DE genes and two control genes were selected at random for validation. Red bars, microarray data; blue bars, qRT-PCR data. The qRT-PCR data were normalized against the reference gene Actb; similar results were obtained using Gapdh as reference (not shown). The similarity of the expression patterns (up-and downregulation) between the microarray and qRT-PCR analyses confirmed the results of the microarray. Error bars represent standard errors from six experiments. Control genes' descriptions are follows; Atp5f1, ATP synthase, H + transporting, mitochondrial Fo complex, subunit B1; Bche, butyrylcholinesterase.

nAChRs Expressions
To confirm the expression of nAChRs in the treated cultures at 16 DIV, mRNA expression of Chrna3, Chrna4, and Chrna7 (coding for α3, α4, and α7 nAChRs, respectively) was determined using qRT-PCR. Figure 6 shows their relative expression as measured by PCR and microarray. Expression of these mRNAs after exposure to NIC, ACE or IMI was very similar to that in the control group and their FCs were log2 −0.44 to log2 0.29. Figure 5. Confirmation of the microarray data using quantitative real-time PCR (qRT-PCR). Twenty DE genes and two control genes were selected at random for validation. Red bars, microarray data; blue bars, qRT-PCR data. The qRT-PCR data were normalized against the reference gene Actb; similar results were obtained using Gapdh as reference (not shown). The similarity of the expression patterns (up-and downregulation) between the microarray and qRT-PCR analyses confirmed the results of the microarray. Error bars represent standard errors from six experiments. Control genes' descriptions are follows; Atp5f1, ATP synthase, H + transporting, mitochondrial Fo complex, subunit B1; Bche, butyrylcholinesterase.

nAChRs Expressions
To confirm the expression of nAChRs in the treated cultures at 16 DIV, mRNA expression of Chrna3, Chrna4, and Chrna7 (coding for α3, α4, and α7 nAChRs, respectively) was determined using qRT-PCR. Figure 6 shows their relative expression as measured by PCR and microarray. Expression of these mRNAs after exposure to NIC, ACE or IMI was very similar to that in the control group and their FCs were log 2 −0.44 to log 2 0.29. Figure 6. Relative mRNA expression of three nAChRs in cerebellar cultures exposed to nicotine, acetamiprid or imidacloprid for 14 days, at 16 days in vitro, versus control, measured using qRT-PCR and microarray. Error bars represent standard error of six independent experiments. T-tests showed no significant differences in expression between PCR and microarray, or between control versus nicotine (CvN), acetamiprid (CvA), or imidacloprid (CvI).

Bioinformatics Analysis of DE Genes for CvN, CvA, and CvI
First we applied the PANTHER Classification System to the three group lists of DE genes to gain the overview of the functions of these genes (Figure 7). Most patterns were similar between the three groups and comprised 10 categories: "cellular process", "metabolic process", "developmental process", "biological regulation", "localization", "multicellular organismal process", "response to stimulus", "biological adhesion", "cellular component organization or biogenesis", and "immune system process". The three GO categories of biological processes ("apoptotic process", "locomotion", and "reproduction") were common to CvA and CvI. Only "growth" was unique to CvI (Figure 7). Details of PANTHER classifications are shown in supplemental Table S3.
We next applied the DE gene lists (CvN, CvA, and CvI), with their FCs and p values to the MetaCore system, and performed bioinformatics enrichment analysis on these data, which revealed 10 lists of significant-scoring GO processes and biomarkers of diseases (Tables 5-7). Genes for two types of calcium channels, Cacna1h (voltage-dependent T-type calcium channel subunit α-1H) and Cacng1 (voltage-dependent calcium channel γ-1 subunit) showed differential expression in all three groups and frequently included various categories in related genes, such as "calcium or divalent ion transport" (Tables 5-7, GO processes) and "epilepsy", "brain diseases", and "seizures" (Tables 5-7, diseases). Also, F2rl2 (coagulation factor II (thrombin) receptor-like 2), alias Par3 (proteinase-activated receptor 3) and belonging to a multifunctional subfamily of G-protein coupled receptors (GPCRs), was a DE gene common to all three groups and related to multiple categories, such as "ion transport", "brain diseases", "autistic disorder", and "child development disorders" (Table 5-7). DE genes unique to each group were also related to important categories. These included Ihh (Indian hedgehog; CvN), related to "cell-cell adhesion" and "brain diseases" (Table 5), Hrh2 (histamine receptor H2; CvA), and Htr2c (serotonin receptor 2C; CvI), both related to important categories of GO processes and diseases (Tables 6 and 7). . Relative mRNA expression of three nAChRs in cerebellar cultures exposed to nicotine, acetamiprid or imidacloprid for 14 days, at 16 days in vitro, versus control, measured using qRT-PCR and microarray. Error bars represent standard error of six independent experiments. T-tests showed no significant differences in expression between PCR and microarray, or between control versus nicotine (CvN), acetamiprid (CvA), or imidacloprid (CvI).

Bioinformatics Analysis of DE Genes for CvN, CvA, and CvI
First we applied the PANTHER Classification System to the three group lists of DE genes to gain the overview of the functions of these genes (Figure 7). Most patterns were similar between the three groups and comprised 10 categories: "cellular process", "metabolic process", "developmental process", "biological regulation", "localization", "multicellular organismal process", "response to stimulus", "biological adhesion", "cellular component organization or biogenesis", and "immune system process". The three GO categories of biological processes ("apoptotic process", "locomotion", and "reproduction") were common to CvA and CvI. Only "growth" was unique to CvI (Figure 7). Details of PANTHER classifications are shown in supplemental Table S3.
We next applied the DE gene lists (CvN, CvA, and CvI), with their FCs and p values to the MetaCore system, and performed bioinformatics enrichment analysis on these data, which revealed 10 lists of significant-scoring GO processes and biomarkers of diseases (Tables 5-7). Genes for two types of calcium channels, Cacna1h (voltage-dependent T-type calcium channel subunit α-1H) and Cacng1 (voltage-dependent calcium channel γ-1 subunit) showed differential expression in all three groups and frequently included various categories in related genes, such as "calcium or divalent ion transport" (Tables 5-7, GO processes) and "epilepsy", "brain diseases", and "seizures" (Tables 5-7, diseases). Also, F2rl2 (coagulation factor II (thrombin) receptor-like 2), alias Par3 (proteinase-activated receptor 3) and belonging to a multifunctional subfamily of G-protein coupled receptors (GPCRs), was a DE gene common to all three groups and related to multiple categories, such as "ion transport", "brain diseases", "autistic disorder", and "child development disorders" (Tables 5-7). DE genes unique to each group were also related to important categories. These included Ihh (Indian hedgehog; CvN), related to "cell-cell adhesion" and "brain diseases" (Table 5), Hrh2 (histamine receptor H2; CvA), and Htr2c (serotonin receptor 2C; CvI), both related to important categories of GO processes and diseases (Tables 6 and 7).   pie charts A, B, and C, respectively. Each gene list was applied to PANTHER software. GO biological processes are shown in pie charts, CvN (A), CvA (B), and CvI (C). A number of genes were not recognized by the PANTHER system, including Ndufaf2 and Sdr42e2 (common to all three lists); Car3 and LOC317356 (CvN); LOC100363332 (CvA); Cyp4a1, Klk1c7, LOC100364862, LOC100912563, and LOC294497 (CvI).   Next, we analyzed three groups of DE gene lists by MetaCore's Compare Experiments Workflow, which is a tool that can be used for comparing experimental data by analyzing their intersections in terms of their mappings onto MetaCore's various ontologies, including GO processes and diseases. The comparison analysis results of three DE gene lists in CvN, CvA, and CvI revealed significant GO processes and diseases categories ( Table 8). The GO processes were associated with "calcium ion transport", "divalent metal ion transport", "divalent inorganic cation transport", etc., most of which were correlated to common DE genes including calcium channels (Cacna1h and Cacng1) and GPCRs (F2rl2). The listed diseases are "epilepsy", "seizures", "autistic disorder", "pervasive child development disorders", and "mental disorders diagnosed in childhood", most of which were also correlated with calcium channels and GPCRs, same as GO processes.

Discussion
This is the first in vitro investigation of transcriptome changes in rat neonatal cerebellar neuron-enriched cultures following long-term exposure (14 days) to NIC and the neonicotinoids ACE and IMI, which are widely used in pesticides. Our results indicate that alterations in transcriptome profiles occur after exposure to these agents.
The DE gene lists of the CvN, CvA, and CvI comparisons included genes unique to each comparison and those genes common to two or more comparisons and support our previous data regarding NIC-like effects of ACE and IMI [19]. Neonicotinoids and NIC all excite small cerebellar neurons, but the firing patterns, firing peaks, and excited cell numbers differ between the different agents.
In the DE genes common to the three groups, those for two types of calcium channels, Cacna1h and Cacng1, are important for neuronal activity. The Cacna1h gene is registered as a strong candidate in the autism-related gene database, SFARI [39]. Mutations in this gene are also found in patients with epilepsy [40,41]. In cerebellar Purkinje cells, parallel fiber excitatory postsynaptic potentials are regulated by K channels and the Cav3.2 (Cacna1h) calcium channel [42]. The Cacng1 gene expressed in human fetal and adult brain [43] and developing mouse cerebellum [44]. A mutation of this gene was also proposed as a candidate for epilepsy [45] and its expression in hippocampal CA1 is altered in drug-induced epilepsy [46]. Another common DE gene, F2rl2 (Par3) belongs to an important subfamily of GPCRs, that are correlated with various phases of neurodevelopment [47][48][49][50], including axon-dendrite differentiation [51]. The remaining DE genes common to the three experimental groups were Cramp1l (alias: hematological and neurological expressed 1-like protein), which is a DNA and chromatin binding protein; Tada2b, a transcriptional adaptor protein; Lmod3, a fetally expressed member of the leiomodin family, which increases actin polymerization rate; Ndufaf2, which is NADH:ubiquinone oxidoreductase complex assembly factor 2; Sdr42e2, one of the short-chain dehydrogenases/reductases; and Unc45b, which acts as a co-chaperone for HSP90 and is required for the proper folding of the myosin motor domain [52]. These genes all have basic proliferation or differentiation functions, and expression of each was confirmed in the developing mouse cerebellum database [44] or the human brain database, GTEx Portal [53].
NIC, ACE, and IMI appeared to cause a slight disturbance in the dendritic arborization of Purkinje cells (Figure 2). This disturbance might be correlated to altered expression of the nine genes. Indeed, F2rl2 (Par3) regulates axon and dendrite differentiation [51], and its disruption might affect this important developmental stage. Further pharmacological and mechanistic investigations are needed to identify the role of these DE genes in the development of Purkinje (and other) cells. Expression of Calb1 (calbindin D28k) was very similar between control cultures and each experimental group in the microarrays.
In the SFARI autism database, the following DE genes were registered: Celf6, which codes for a transcription associated protein and was common to CvN and CvI; Magel2 (melanoma antigen family L2), unique to CvI; and Ampd1 (adenosine monophosphate deaminase), also unique to CvI. Three types of CvN unique genes were registered in the autism KB database [54,55]: Cadm3 (cell adhesion molecule), Gpr83 (G protein-coupled receptor), and Acta1 (actin α1). Four CvA-unique genes were registered in the same database: Napb (related to amyotrophy), Rasl10b, Iqcf1, and Myog (transcription factors). Finally, three CvI-unique are also registered in the Autism KB database: Htr2c, Rbfox2 (an RNA binding protein), and Txk (a tyrosine kinase). The functions and correlations of these disorder-related genes are not yet clear, but given the developmental nature of autism, they may be involved in brain development.
The PANTHER ontology analysis of the three DE gene lists (CvN, CvA, CvI) showed similar and fundamental gene categories for brain development (Figure 7 and supplemental Table S3). Independent analysis of each DE gene list using MetaCore indicated similarities and differences in significant categories of GO processes and diseases (Tables 5-7). In the process category, all three experimental groups included "calcium ion transport", which is essential for neuronal activity. Within the diseases category, "epilepsy", "autistic disorders", and "child development disorders, pervasive" were common to the three groups. Moreover, comparing analyses of all three groups confirmed these results (Table 8), and other essential GO process categories ("cellular response to potassium ion", "metal ion transport", thrombin receptor signaling", "positive regulation of exocytosis", "membrane depolarization") and disease categories ("mental disorders diagnosed in childhood", "mitochondrial complex 1 deficiency", "delirium", "neurologic manifestations") are listed with significant p values.
In each enrichment analysis, two GO processes ("glucocorticoid biosynthetic process", "cellular component morphogenesis") were common to CvN and CvA, and two diseases ("Tourette syndrome", "neurotoxicity syndromes") were common to CvA and CvI. Other disease categories, unique to each group, included "movement disorders" (CvN), "communication disorders" (CvA), and "dyskinesia" (CvI), etc. These analyses suggest that chronic exposure to NIC, ACE or IMI might alter gene expression profiles, which are important for normal brain development. It will be necessary to confirm this in future studies.
A limitation of the present study was that we could not determine whether the altered gene expression profiles were caused by the direct action of NIC, ACE and IMI on nAChRs. However, we did confirm the mRNA expression of α3, α4, and α7 nAChR subunits, both here ( Figure 6) and in our previous study [19]. Our previous report indicated that NIC-, ACE-, or IMI-induced Ca 2+ -influx in cerebellar neurons was inhibited by the specific nAChR antagonists dihydro-β-erythroidine (α4β2 specific) and α-bungarotoxin (α7 specific), which suggests that NIC and neonicotinoids act directly at these receptors. Furthermore, we cannot rule out the possibility that some toxic metabolites of ACE or IMI are produced during cultivation, and both the agents and their metabolites might act at these receptors.
Other reports have suggested that neonicotinoids have direct actions on mammalian nAChRs. IMI and clothianidin directly activate human α4β2 nAChRs [9], and modify the receptor's responses to ACh, even at a concentration below the activation threshold of the receptors.
ACE and IMI have some potency at mammalian α3 and α4β2 nAChRs [3]. Although the authors concluded that IMI and ACE were inactive at α7 nAChRs, their data indicate that ACE and IMI actually inhibited the binding of α-bungarotoxin to rat or human α7 nAChRs. Moreover, α-bungarotoxin antagonized IMI-induced activation of stellate cells of the mouse cochlear nucleus, which suggests that IMI binds to mouse α7 or α9α10 nAChRs [56] and IMI activated mammalian α7 nAChRs expressed in Xenopus oocytes [57]. In excitation assays of chicken α7 nAChRs, two neonicotinoids IMI and thiacloprid actually activated the receptors, even though the excitation peaks were lower than those by ACh or NIC [58]. Those authors expected similar results with human α7, because of its similarity to the avian α7 (93.8% in the extracellular agonist binding regions). Together, these data support the notion that IMI and ACE bind to mammalian α7 nAChRs.
ACh and other ligands often cause desensitization of nAChRs, including α4β2 and α7, even concentrations below the receptor's activation threshold [25,68]. In addition, many factors change the kinetics of receptor desensitization, including subunit composition, temperature, post-translational modification, types of ligand, duration of ligand exposure, and interactions with modulator proteins [69][70][71]. Together, this suggests that the transcriptome profile alterations induced by NIC, ACE, and IMI might be caused by activation and/or desensitization of nAChRs.
In the mammalian nervous system, nAChRs are expressed in glia as well as neurons. For example, α7 subtypes are expressed in astrocytes [72] and microglia [73], and α3, α4, α5, α7, β2, and β4 subtypes are found in oligodendrocyte precursor cells [74]. Such reports, together with the present results, suggest that these glial nAChRs have important regulatory functions that can be disrupted by NIC, ACE, and IMI.
Recently, Meijer et al. reported that several pesticides disturbed KCl-induced Ca 2+ -influx via voltage-gated calcium channels in PC12 or primary cortical cells, but that IMI and carbaryl did not [75,76]. In contrast, our previous report indicated that IMI disturbed KCl-induced Ca 2+ -influx after nAChR activation in 33-40% of rat cerebellar neurons. This discrepancy may be a factor of the time course examined, or differences in the number, subclass, or localization of nAChRs in PC12 and primary cortical cells. Alternatively, the difference may be due to the fact that IMI induced Ca 2+ -influx in small proportion of rat cerebellar neurons, whereas KCl induced significant Ca 2+ -influx in the majority of PC12 cells and also our cerebellar neurons.
To date, there has been no definitive evidence for mammalian developmental neurotoxicity caused by neonicotinoids [77], because of the widely held view that they have a low affinity for mammalian nAChRs, penetrate poorly into the mammalian brain through the BBB, and have a lack of effects in common with NIC. However, increasing evidence to the contrary [11][12][13][14][15][16][17][18][19], and consideration of some neonicotinoid metabolites showing high affinities for mammalian nAChRs [3], suggests that neonicotinoids have potential developmental neurotoxicity.
In utero exposure to NIC in combination with an organophosphate produces significant sensorimotor deficits in rat offspring [78]. In 2012, the American Academy of Pediatrics stated that "epidemiologic evidence demonstrates associations between early life exposure to pesticides and pediatric cancers, decreased cognitive function, and behavioral problems [79]," and "children's exposures to pesticides should be limited as much as possible [80]". This statement affirms the neurotoxicity of organophosphate pesticides, which disturb the acetylcholinesterase and cholinergic systems. Actions of NIC and neonicotinoids on nAChRs may increase the cholinergic neurotoxicity of organophosphates.
In a recent report, three major classes of pesticides were shown to be commonly present in urine samples from 3-year-old children in Japan (n = 223) [81]. In total, 79.8% of the samples were positive for at least one neonicotinoid (median 3.97 nM and maximum 308.18 nM), and 100% were positive for at least one organophosphate metabolite (median 211.45 nM and maximum 2232.78 nM) and pyrethroid metabolites, 3-PBA (median 1.01 µg/L and maximum 25.29 µg/L). Studies such as this highlight the urgency for investigations into the effects of these neonicotinoids and other pesticides on human health. Although the pesticide concentrations were within the acceptable daily intake, the pesticide registration system in Japan does not include developmental neurotoxicity tests, and combined effects of various pesticides have not yet been examined.
There are several limitations to this study; for example, we examined only one brain region. In several preliminary experiments, we exposed cultured cells from the cerebellum, hippocampus, and cerebral cortex to 1~10 µM NIC, ACE, or IMI at 10 DIV, which resulted in variable transcriptome disturbances, not all of which were statistically significant (not shown). Moreover, we did not examine the protein expressions of DE genes. Despite these limitations, our study suggests that neonicotinoids might affect mammalian brain development. Further research, including dose-dependency and time course studies, functional assays (e.g., Ca2 + imaging), protein expression measurements, and in vivo experiments, are needed to define the effects of neonicotinoids on the developing mammalian brain.

Conclusions
Chronic low-dose exposure of the neonicotinoids ACE, IMI can induce certain alterations in the transcriptome of the developing rodent brain. Several genes known to be essential for brain development were up or downregulated after exposure of cerebellar cultures to NIC, ACE, and IMI, including two types of calcium channels (Cacna1h and Cacng1) and a subfamily of GPCRs (F2rl2/Par3). Bioinformatic analysis suggested a correlation between these transcriptome alterations and developmental disorders such as epilepsy and autistic disorders. In addition, immunocytochemistry in Purkinje cells revealed that chronic exposure to NIC, ACE, or IMI induced small disturbances in dendritic arborization that might correlate to the observed transcriptome alterations. Taken together, our results indicate that neonicotinoids used as pesticides may have some adverse effects on the developing mammalian brain. Even though the affinities of IMI and ACE for mammalian nAChRs are lower than to those for insect nAChRs, the effects of neonicotinoids on the developing brain should be carefully investigated.