Organophosphorus flame retardants are developmental neurotoxicants in a rat primary brainsphere in vitro model

Due to regulatory bans and voluntary substitutions, halogenated polybrominated diphenyl ether (PBDE) flame retardants (FR) are increasingly substituted by mainly organophosphorus FR (OPFR). Leveraging a 3D rat primary neural organotypic in vitro model (rat brainsphere), we compare developmental neurotoxic effects of BDE-47—the most abundant PBDE congener—with four OPFR (isopropylated phenyl phosphate—IPP, triphenyl phosphate—TPHP, isodecyl diphenyl phosphate—IDDP, and tricresyl phosphate (also known as trimethyl phenyl phosphate)—TMPP). Employing mass spectroscopy-based metabolomics and transcriptomics, we observe at similar human-relevant non-cytotoxic concentrations (0.1–5 µM) stronger developmental neurotoxic effects by OPFR. This includes toxicity to neurons in the low µM range; all FR decrease the neurotransmitters glutamate and GABA (except BDE-47 and TPHP). Furthermore, n-acetyl aspartate (NAA), considered a neurologic diagnostic molecule, was decreased by all OPFR. At similar concentrations, the FR currently in use decreased plasma membrane dopamine active transporter expression, while BDE-47 did not. Several findings suggest astrogliosis induced by the OPFR, but not BDE-47. At the 5 µM concentrations, the OPFR more than BDE-47 interfered with myelination. An increase of cytokine gene and receptor expressions suggests that exposure to OPFR may induce an inflammatory response. Pathway/category overrepresentation shows disruption in 1) transmission of action potentials, cell–cell signaling, synaptic transmission, receptor signaling, (2) immune response, inflammation, defense response, (3) cell cycle and (4) lipids metabolism and transportation. Taken together, this appears to be a case of regretful substitution with substances not less developmentally neurotoxic in a primary rat 3D model. Electronic supplementary material The online version of this article (10.1007/s00204-020-02903-2) contains supplementary material, which is available to authorized users.


Introduction
Flame retardants (FR) are a group of compounds, which are added to consumer products, including upholstered furniture, electrical devices, baby products, textiles and plastics, to restrain or delay flame propagation to prevent fire spreading (Dishaw et al. 2014b;EPA US 2005;Jarema et al. 2015). The global FR consumption has surpassed 2 million tons and is yet expected to increase due to international flammability standards (Ceresana 2018). However, FR exhibit characteristics similar to environmental toxicants, such as heavy metals, air pollutants and pesticides that are well recognized as hazardous to human health and inducers of neurodevelopmental damage. Prior to 2005, halogenated polybrominated diphenyl ethers (PBDEs) were the primary FR used in the USA. However, the halogenated FR have been linked to the development of cancer, endocrine disruption, immunotoxicity, reproductive toxicity, and fetal and child development perturbation (Birnbaum and Staskal 2004;Costa and Giordano 2007;Roze et al. 2009;Shaw et al. 2010). In particular, 2,2′,4,4′-tetrabromodiphenyl ether (BDE-47)-the most abundant PBDE congener in the environment and human serum (Birnbaum and Staskal 2004;EPA US 2008)-has been shown to affect the adult and developing nervous system (Dingemans et al. 2007;Eriksson et al. 2002).
Due to these human health concerns, PBDEs have been banned in Europe and phased out in the USA (Birnbaum and Staskal 2004;Feo et al. 2013), and mainly been replaced by organophosphorus FR (OPFR) (Blum et al. 2019;Dishaw et al. 2014aDishaw et al. , 2014bStapleton et al. 2014). Following the phase out of PBDEs, there is increasing evidence showing higher exposure to OPFRs compared to PBDEs from hand wipes in toddlers, and house dust suggesting that the magnitude of exposure via hand-to-mouth and dermal transfer pathways is potentially greater for OPFRs than for PBDEs (Blum et al. 2019).
The use of OPFRs is further anticipated to increase following the recent proposal by the European Commission to prohibit the class of organohalogen chemicals in electronic display enclosures and stands effective since April 2021 (Commission 2018). Hence, it is likely that the manufacturers will explore the potential for use of OPFRs in televisions and other electronics as alternative methods to meet flammability codes.
Although there has been an increase in the use of OPFRs, there is still relatively limited information on their potential health effects. This is of concern, as OPFRs bear some structural similarities to organophosphate pesticides that are well known to induce (developmental) neurotoxicity (Burke et al. 2017;Grandjean and Landrigan 2014;Mie et al. 2018). Especially, the use of FR in baby products and the exposure to children are of concern as the developing brain is much more vulnerable to environmental perturbation than the brain of adults (O'Rahilly and Muller 2008;Smirnova et al. 2014).
The exposure to industrial chemicals, including FR, and drugs during early development has been associated with the occurrence of neurodevelopmental disorders in children such as autism, mental retardation, dyslexia, epilepsy or mental deficit Landrigan 2006, 2014;Rice and Barone 2000;Zhong et al. 2020).
Current DNT testing guidelines for chemicals and pesticides (EPA U 1998;OECD 2007OECD , 2011 are based on traditional in vivo animal studies. These guidelines require elaborated, lengthy and costly study protocols that are incompatible with the assessment of large numbers of chemicals in addition to elicit uncertain predictivity for human risks (Smirnova et al. 2014). Hence, the screening of chemicals requires the employment of more reliable, cheap and fast tools capable of determining DNT potential to ensure the safety of children's health (Bal-Price et al. 2015a, 2018Bjorling-Poulsen et al. 2008;Fritsche et al. 2017).
The use of in vitro, in silico and non-mammalian species-based methodologies and models has been proposed to enhance the DNT assessment in terms of cost, time and mechanistic understanding (Bal-Price et al. 2015a, 2018Coecke et al. 2007;Smirnova et al. 2014). Several promising in vitro cell models have been developed; however, most of them still have difficulties in simulating the complex structure of the developing brain. Three-dimensional (3D) primary organotypic cell cultures have the advantage of reproducing the complex multicellular environment that closely resembles in vivo conditions (Alepee et al. 2014;Pamies et al. 2017;Sundstrom et al. 2005). The 3D rat primary neural organotypic in vitro model, used in this study, is a model able to recreate the CNS in vivo structural characteristics and biochemical signaling (Forsby et al. 2009;Honegger and Monnet-Tschudi 2001;van Vliet et al. 2007). It consists of most of the relevant cell types in the brain such as neurons, astrocytes, oligodendrocytes and microglia . The model has been extensively characterized by immunohistochemistry, electrophysiology, pharmacological behavior and expression of neurodevelopment marker genes (van Vliet et al. 2007). The model is considered mature after 21-28 days in vitro (DIV), as electrical activity, synaptogenesis and myelination are robust at this time. The DNT consensus process identified the rat aggregating cell model among the most representative models for DNT studies .
Taking into consideration the importance of identifying and characterizing potential DNT chemicals by applying suitable testing methods, this study aimed to investigate the DNT potential of four OPFR (isopropylated phenyl phosphate-IPP, triphenyl phosphate-TPHP, isodecyl diphenyl phosphate-IDDP, and tricresyl phosphate (also known as trimethyl phenyl phosphate)-TMPP) relative to one replaced PBDE (BDE-47) using metabolomics and transcriptomics approaches. The exposure to the FR induced significant alterations in gene expression and metabolite levels, with stronger effect after exposure to OPFRs. The major effects were observed for genes and metabolites associated with the neurotransmitter glutamate and its receptors, followed by general neuronal markers and reactivated glial cells. Many of the alterations could be linked to existing DNT adverse outcome pathways (AOPs) (Sachana et al. 2018;Spinu et al. 2019;Wang et al. 2018), which increase the concern of OPFRs usage as replacements of PBDE and imply the need for additional assessment of these compounds.

Animals
Pregnant female Sprague-Dawley rats were used as the source of embryonic brain tissue for dissociation and reaggregation in vitro. The animals were kept in the Johns Hopkins Bloomberg School of Public Health Animals Resources Facility for 48 h. They were housed individually under standard laboratory conditions with controlled temperature of 22 °C and a 12-h light/dark photoperiod. Food and water were supplied ad libitum. Housing and experimental protocols were approved by the Institutional Animal Care and Use Committee (IACUC) (Protocol RA15H122).
Sixteen-day pregnant females were anesthetized by inhalation of tribromoethanol (Sigma-Aldrich) following immediate decapitation to minimize any pain or distress. Aseptic conditions were maintained throughout the procedure to prevent contamination of tissue cell cultures. The rats were sterilized with 70% ethanol and an incision was made through the skin over the midline of the abdomen for the removal of the uterus containing the fetuses. The fetuses were then excised from the uterus and the whole brain was dissected out for the preparation of re-aggregating brain cell cultures.

Chemical treatments
Five FR (Table 1) and corresponding vehicle control (DMSO) were studied. FR were supplied by the Division of National Toxicology Program, U.S. National Institute of Environmental Health Sciences (NTP, NIEHS, Research Triangle Park, NC). Stock solutions were prepared in DMSO (Sigma-Aldrich). The final concentration of DMSO did not exceed 0.1% (v/v), which is non-toxic in the rat brainsphere model.
At DIV7, aggregating cell cultures were pooled and distributed over six-well plates, each containing approximately 100 aggregates in 2 ml of serum-free media. The cultures were then exposed to the chemicals up to DIV14 and DIV21, a period that covers crucial developmental processes and different stages of cell maturation (Fig. 1). The concentrations of FR were chosen based on preliminary range-finding experiments, where wide ranges of concentrations were tested using cell viability assay (resazurin). In final experiments, three non-cytotoxic concentrations at DIV14 (0.1, 1, and 5 µM) were selected for metabolomics and gene expression analysis. Importantly, using in vitro to in vivo extrapolation suggests that these are relevant human exposures (Blum et al. 2019). Half of the culture media were exchanged with new chemical treatment twice per week. One non-cytotoxic concentration at DIV21 (1 µM) of IPP was selected for transcriptomics analyses.

Assessment of cell viability
Cell viability was determined after exposure to the selected chemicals (0.1-20 µM) at DIV14 and DIV21 using the resazurin reduction assay (O'Brien et al. 2000). The blue colored dye resazurin is reduced to fluorescent resorufin by redox reactions in viable cells with active metabolism. Resazurin (0.01 mg/ml, Sigma-Aldrich) in PBS was added directly to the six-well plates, without removing the medium at the end of the period of exposure to the tested compounds. The plates were incubated for 2 h at 37 °C, 10% CO 2 . After incubation, 100 µl of medium from each sample was transferred to a 96-well plate and the fluorescence of the resazurin metabolite, resorufin, was measured at 530 nm/590 nm (excitation/emission) in a multi-well fluorometric reader (Cytofluor Multi-well Plate Reader Series 4000, Perseptive Biosystems). The differentiation medium was incubated with resazurin in parallel as a blank control. Cell viability was calculated as % of fluorescence intensity relative to solvent-treated controls after subtracting average blanks. Data are presented as mean ± SD of at least three independent experiments performed in two to four replicates. Differences between treated and DMSO control groups were assessed by one-way ANOVA (GraphPad Prism 8.4.3), followed by Dunnett's multiple comparison post hoc test including correction for multiple testing. Post hoc test was only performed vs. controls. Statistical significance is indicated as follows: *p < 0.05 (treated vs. control).

RNA purification, reverse transcription, and quantitative real-time PCR (RT-qPCR)
Cell samples were lysed for mRNA expression analysis and total RNA extraction was carried out according to the manufacturer's protocol of RNeasy Mini Kit (Qiagen). RNA integrity was assessed with the Nanodrop 2000 (Thermo Scientific) UV-Vis spectrophotometer at 260 nm. Reverse transcription was performed as follows: 500 ng RNA was incubated with 2.5 mM PCR Nucleotide Mix (Promega) and 12.5 µg/ml random primers (Promega) for 5 min at 65 °C using a Techne PCR system. Subsequently, 2 units/µl RNa-seOut inhibitor (Thermo Fisher Scientific), 10 units/µl Moloney murine leukemia virus (M-MLV) reverse transcriptase (Promega) and the samples were incubated for 10 min at 25 °C for annealing, 60 min at 37 °C for cDNA synthesis and 15 min at 70 °C for inactivation of enzymes. RNasefree DNase set (Qiagen) was used to avoid contamination with DNA. cDNA was diluted 1:5 and quantitative RT-PCR was performed using the Fast Applied Biosystems 7500 System (Life Technologies). Genes were selected (Table S1) based on previous studies (Hogberg et al. 2009 and transcriptomics data for IPP. The expression was measured using TaqMan gene expression assay (Life Technologies)

Fig. 1
Schematic experimental design of rat brainsphere model exposed during development to FR at DIV7. Samples were collected at DIV14 and 21 for metabolomics and transcriptomics analysis and FastStart Universal Probe Master mix (ROX) (Roche) according to the manufactures protocol. Relative RNA quantification was performed using the comparative CT method, normalizing the data to a standard calibrator (a mixture of samples from the different time points of the cell proliferation and differentiation), and to the 18S rRNA content (Schmittgen and Livak 2008). Data were calculated and presented as average log 2 -fold change in each independent experiment ± SD of at least three independent experiments performed in two to three replicates. Differences between treated and non-treated groups were assessed by two-way ANOVA (GraphPad Prism 8.4.3), followed by Bonferroni's comparison post hoc test including correction for multiple testing. Post hoc test was only performed vs. controls. Statistical significance is indicated as follows: *p < 0.05.

Transcriptomics sample preparation and analyses
Perturbations in transcriptome were analyzed by microarray after exposure to 1 µM of IPP or solvent control (DMSO) from DIV7 to DIV21. Transcriptomics was performed in triplicate from one experiment. 100 ng of total RNA from treated and control cells was converted into cDNA and then into labeled cRNA using Agilent LowInput QuickAmp Labeling Kit (Agilent). The resulting cRNA was labeled with Cy3. Labeled cRNAs were then purified, and RNA concentration and dye incorporation were measured using Nanodrop 2000 spectrophotometer (Thermo Scientific). Hybridization to SurePrint G3 Rat Gene Expression 8 × 60 k (Agilent, Product No. G4853A, Grid No. 028279) was conducted following the manufacturer's protocol. Microarrays were scanned with an Agilent DNA microarray scanner. Feature Extraction (12.0.0.7 version, Agilent) was used to calculate the signal intensity and ratios. All arrays met each of the ten quality evaluation metrics parameters of Feature Extraction ("Good"). After deleting non-detected probes and quantile-normalization, statistical and pathway overrepresentation analyses were performed with activated mean centering and scaling option (GeneSpring V13.1, Agilent).

Metabolomics sample preparation and analyses
Cells were collected in 1.5 ml Eppendorf tubes and washed three times with ice-cold PBS. After removal of PBS, icecold high-purity methanol (MeOH) (Sigma-Aldrich) was added. Cells were stored at -80 ℃ until use. For metabolite extraction, 75 µl of HPLC-grade water was added to the 300 µl MeOH to allow 80:20 v/v mixture of high purity MeOH:water. The cells were disrupted using an ultrasound sonicator (Qsonica, CT, USA) for 10-20 s until no more intact cell could be detected. The total protein content of the homogenates was quantified according to the manufacturer's protocol of the Bradford assay (Bio Rad) to control for potential differences in tissue quantities. After being stored at − 20 ℃ for at least 2 h to precipitate the proteins, the tubes were centrifuged at 14,000 rcf for 10 min at 4 ℃. The supernatant was transferred to a new Eppendorf tube and evaporated to dryness at room temperature in a speedvac concentrator (Savant, Thermo Fisher Scientific). The dried samples were reconstituted in 100 µl of 60% MeOH with 0.1% formic acid (FA) (Sigma-Aldrich) and transferred to plastic vials for LC-MS measurements.
Metabolite separation was achieved using an Agilent 1260 series HPLC system (Agilent, Santa Clara, CA). The injection volume of each sample was 5 µL and the column was maintained at 35 °C. QCs (pool of all samples within the experiment) and standards were run at the beginning and the end of each sequence and every four sample runs to monitor shift in the retention time on the column. For negative mode, a Cogent Diamond Hydride TM (MicroSolv, Eatontown, NJ, USA, Cat# 70,000-15P-2) aqueous normal phase (ANP) column (150 × 2.1 mm i.d., 4 µm particle size, 100 Å pore size) was used. The run time was 25 min at a flow rate of 0.4 ml/min. Chromatography was performed using solvent A (50% MeOH/50% water/0.05% FA) and solvent B (90% acetonitrile with 5 mM ammonium acetate). The gradient was: 0 min, 100% B; 20-25 min, 40% B; post-run time for equilibration, 10 min in 100% B. For positive mode, a Targa C18 reverse phase column (50 × 2.1 mm i.d., 3 µm particle size, 120 Å pore size, Higgens Analytical Mountain View, CA, USA, Cat# TS-0521-C183) was used. The run time was 25 min at a flow rate of 0.3 ml/min. Chromatography was performed using solvent A (Water with 0.1% formic acid) and solvent B (98% acetonitrile/2% water with 0.1% formic acid). The gradient was: 0 min, 2% B; 20-25 min, 100% B; post-run time for equilibration, 5 min in 2% B.
A 6520B Q-TOF LC-MS (Agilent, Santa Clara, CA) was operated in both positive and negative electrospray ionization (ESI) modes with an acquisition rate of 1.5 spectra/s in extended dynamic range (1700 m/z, 2 GHz). The spectra were internally mass calibrated in real time by continuous infusion of a reference mass solution using an isocratic pump connected to a dual sprayer feeding into an electrospray ionization source. Data were acquired with MassHunter Acquisition software B.05.01 and further processed with Agilent Profinder B.09 (Agilent, Santa Clara, CA).
For the data processing and chemometric analysis of the LC-MS untargeted data, the acquired raw data files (.d files) were first checked for quality in MassHunter Qualitative Analysis software (Agilent, version 7.0). Reproducibility of chromatograms was visually inspected by overlaying the total ion chromatograms (TICs) of all samples. Data files that exhibit outlier peaks, i.e., replicates with very dissimilar chromatograms (e.g., significant retention time shifts), were excluded for further processing. The raw data files were then converted to mzXML using ProteoWizard 3.0 (Kessner et al. 2008). Raw LC-MS data were analyzed by the MZmine 2 software (Pluskal et al. 2010) for chromatogram deconvolution, peak detection and alignment. The metabolites were called by batch-targeted feature extraction. The putative identification was achieved by online searching for the accurate m/z values of the peaks against HMDB and KEGG databases (Kanehisa and Goto 2000; Wishart et al. 2018). Those peaks were manually inspected for the quality of the extracted ion chromatograms (plausible adduct formation, max mass deviation 5 ppm, isotope ratios and peak shape) and for the remaining duplicate compound names. Data were calculated and presented as average log 2 -fold change in each independent experiment ± SD of at least two independent experiments performed in four to six replicates. Differences between treated and non-treated groups were assessed by two-way ANOVA (GraphPad Prism 8.4.3), followed by Bonferroni's comparison post hoc test including correction for multiple testing. Post hoc test was only performed vs. controls. Statistical significance is indicated as follows: *p < 0.05.

Results
To characterize the metabolic perturbation of FR, 3D rat primary neural cell cultures were obtained from 16-dayold fetal rat brains and were treated with FR from DIV7 until DIV14 or DIV 21 to cover critical periods of differentiation and maturation. The potential DNT effects of FR were assessed through the expression of specific genes selected from a previous work (Hogberg et al. 2009 that serve as markers for the structural and functional development during subsequent stages of neuronal and glial differentiation. Additional genes were selected based on IPP microarray data. Furthermore, untargeted metabolomics was performed using a quadrupole time-offlight liquid chromatography mass spectrometry (Q-TOF LC-MS). Metabolomics involves the analysis of metabolic profiles in living cells in response to physiological alterations triggered by endogenous or exogenous elements, such as chemicals and pathological and developmental factors (Nicholson et al. 2012;van Vliet et al. 2013van Vliet et al. , 2008).

Assessment of cell viability of 3D rat primary neural cell cultures exposed to FR
To determine non-cytotoxic concentrations of FR in rat brainspheres, aggregates were exposed to FR for 7 or 14 days starting from DIV7. Cytotoxicity was assessed using the resazurin cell viability assay at two different time points-DIV14 and DIV21. Initially, a wide range of concentrations for each FR was tested (Fig. 2, Table 2). Time and dose-dependent decrease in cell viability was induced by all FR (DIV14 vs. DIV21). At DIV14, 7 day exposure to 10 µM of BDE-47 (Fig. 2a), IDDP (Fig. 2d) and TMPP (Fig. 2e) induced significant reduction of the cell viability, while TPHP (Fig. 2B) only showed significant decrease at 20 µM. IPP (Fig. 2c) was only tested up to 10 µM where no significant effect was observed at DIV14. At DIV21, 14 days exposure to 5 µM of all FR showed a significant decrease in cell viability ( Fig. 2 and Table 2). The level of cytotoxicity caused by the FR was used as reference for the selection of concentrations for gene expression and metabolomics analyses. In conclusion, 0.1 µM, 1 µM (non-cytotoxic at 14 and 21 DIV) and 5 µM (non-cytotoxic at 14 DIV and lowest observed cytotoxicity effect at 21 DIV, ~ 70-80% viability vs. control) were selected for further experiments.

Exposure to FR significantly altered marker genes involved in neuronal morphology and function
To assess FR effects on neurons, neurofilament 200 (nf-200) as an intermediate filament highly expressed in neurons during the later stages of differentiation was chosen. It is an important cytoskeleton marker, whose expression can be used to detect neuronal morphology changes (Gupta et al. 1999;Tonnaer et al. 2010). The mRNA levels of nf-200 were significantly downregulated after exposure to all FR (5 µM), at DIV14 (no cytotoxicity detected) ( Fig. 3a and suppl. material 1). Further decrease was observed at DIV21 after exposure to 5 µM IPP and TMPP and already 1 µM TPHP exposure significantly decreased the mRNA expression of nf-200 at DIV21 (no cytotoxicity observed). In conclusion, all FR studies were toxic to neurons in the low µM range.
The mRNA expression of two receptors that play crucial roles in neuronal function were addressed to further assess neuronal impairment: subunits (grin1, grin2a and grin2c) of the ionotropic N-methyl D-aspartate receptor (NMDA-R) of the main excitatory neurotransmitter glutamate, (Blanke and VanDongen 2009;Busse et al. 2014) and subunit alpha 1 (gabra1) of the main inhibitory neurotransmitter gammaaminobutyric acid A receptor (GABA A -R) (Sigel and Steinmann 2012).
Exposure to the highest concentration (5 µM) of IPP, IDDP and TMPP significantly decreased mRNA levels of the NMDA-R subunit grin1 at DIV14 and was further decreased at DIV21 (Fig. 3a and suppl. material 1), although at DIV21, already the lower concentration (1 µM) of these FR (IPP, IDDP and TMPP) downregulated the expression of grin1. Exposure to the highest concentration (5 µM) also decreased the mRNA levels of grin1 at 21 DIV for BDE-47, and TPHP. The gene expression of the NMDA-R subunit grin2a was significantly downregulated after exposure to all OPFR both at 14 and 21 DIV (Fig. 3a and suppl. material 1). Exposure to BDE-47 did not alter the expression of grin2a. In contrast, the expression of    subunit grin2c was significantly upregulated after exposure to all FR (5 µM) both at DIV14 and DIV21 (Fig. 3a and suppl. material 1). Already the lower concentration of 1 µM increased the mRNA levels after exposure to TPHP, IPP and TMPP at 14 and/or 21 DIV) ( Fig. 3a and suppl. material 1). In conclusion, with slight differences in active concentrations, all FR decreased grin1and grin2a, but increased grin2c expression that could affect the affinity of the receptor. The mRNA expression of GABA A -R subunit alpha 1 (gabra1) was significantly downregulated after exposure to 5 µM of all OPFR at 14 DIV ( Fig. 3a and suppl. material 1). Already exposure to 1 µM of IDDP and TMPP significantly decreased the mRNA level of gabra1 at DIV14. After longer treatment (DIV21) only exposure to IPP and TMPP (5 µM) continued to decrease in the mRNA level of gabra1, while TPHP and IDDP treated cultures were reversed to control levels. Exposure to BDE-47 did not alter the expression of gabra1. In conclusion, the replacement FRs downregulated gabra1, while the same concentration of BDE-47 had no effect.
Furthermore, the expression of the neuronal marker synapsin 1 (syn1) was evaluated. SYN1 is a protein present in the membrane of synaptic vesicles (Lu et al. 1992). It modulates neurotransmitter release and is believed to exert an important role in the functional maturation of synapses (Harrill et al. 2011). Already the lower concentration of 1 µM decreased the expression of syn1 for most of the OPFR at DIV14 and DIV21 (Fig. 3a and suppl. material 1). Exposure to 5 µM of all the OPFR downregulated the expression of syn1 at both time points. BDE-47 did not have any effect on the expression of syn1 at those concentrations. In conclusion, the expression of syn1 was significantly downregulated after exposure to OPFRs, but not to BDE-47.
To characterize the perturbation of brainspheres on metabolic level, the neural-specific metabolite n-acetyl aspartate (NAA) was measured. NAA is considered a diagnostic molecule for patients with brain damage and neurodegenerative disorders (Alakkas et al. 2019;Baslow et al. 2003;Chitturi et al. 2018;Nordengen et al. 2015). Mass spectrometry revealed that NAA was significantly lower after exposure to IPP (5 µM), IDDP (1 and 5 µM) and TMPP (1 and 5 µM) at 14 and 21 DIV (Fig. 3b and suppl. material 1). l-Aspartic acid, the precursor of NAA, was also significantly lower after exposure to all OPFRs ( Fig. 3b and suppl. material 1). This indicates that OPFRs induce toxicity in rat brainspheres.
In summary, the exposure to OPFR induced stronger effects on NAA, l-aspartic acid and selected genes involved in neuronal morphology and function than BDE-47. In addition, expression of subunits of the NMDA-R were more affected by exposure to OPFR than the subunit of the GABA A -R.

Exposure to OPFR significantly decreased the level of neurotransmitters and downregulated the expression of genes involved in their production and transportation
To further understand the effect of FR on neuronal cells, we evaluated expression of the enzymes catalyzing neurotransmitter production and neurotransmitter transporters. Neurotransmitters GABA and glutamate (l-glutamic acid) play important regulatory roles in neuronal activities in the brain (Busse et al. 2014;Sigel and Steinmann 2012). Therefore, the expression of genes encoding two enzymes, glutamate decarboxylase gad1 and gad2 that catalyze the synthesis of GABA from glutamate and the key enzyme responsible for catalyzing GABA degradation, 4-aminobutyrate transaminase (abat), was investigated. Treatment with 5 µM of all OPFRs, but not BDE-47, significantly decreased the expression of gad1 and gad2 at both DIV14 and DIV21 ( Fig. 4A and suppl. material 1). The expression of abat was only affected by exposure to 5 µM IPP at DIV21.
Metabolomics analysis showed a significant decrease in the neurotransmitter glutamate at DIV14 after exposure to all FR ( Fig. 4b and suppl. material 1). The strongest effect was observed after exposure to IDDP and TMPP as already the lowest concentration of 0.1 µM was effective. After DIV21, the amount of glutamate was restored except for the highest concentration of IPP, IDDP and TMPP (5 µM). Similar effects were observed on GABA: exposure to IDDP and TMPP decreased levels of GABA already at the lowest concentration (0.1 µM) and was restored at DIV21 except for 5 µM of IPP, IDDP and TMPP ( Fig. 4b and suppl. material 1). Exposure to BDE-47 and TPHP did not alter the levels of GABA. Moreover, α-ketoglutaric acid a derivate of glutamate were significantly decreased at DIV21 after exposure to all FR except IPP, with TPHP having the strongest effect ( Fig. 4b and suppl. material 1). In conclusion, all FR decreased the neurotransmitters glutamate, GABA (except BDE-47 and TPHP) and α-ketoglutaric acid (except IPP).
Next, we addressed the dopaminergic neurotransmitter system. The plasma membrane dopamine active transporter (dat-slc6a3) and the vesicular monoamine transporter (vmat1-slc18a1 and vmat2-slc18a2) are critical components for dopamine neurotransmission: The plasma membrane dopamine active transporter terminates the signal of the neurotransmitter by providing dopamine reuptake from the synaptic cleft, while vesicular monoamine transporter packages cytoplasmic dopamine into vesicles for storage and future use (Miller et al. 1999;Yamamoto et al. 2007). Exposure to all OPFR (5 µM) at DIV14 decreased the mRNA expression of dat ( Fig. 4A and suppl. material 1). At DIV21, the expression of dat after exposure to IDDP was returned to control levels, while the decrease was similar to DIV14 levels for TPHP, IPP and TMPP. BDE-47 treatment did not modify the gene expression of dat. In conclusion, at similar concentrations, the FR currently in use decreased plasma membrane dopamine active transporter expression, while BDE-47 did not.
Treatment with OPFRs did not alter either vmat1 or vmat2 mRNA levels at DIV14 (Fig. 4a and suppl. material 1). Interestingly, treatment with BDE-47, at DIV14 significantly increased the vmat1 gene expression. At 21 DIV significant downregulation of vmat1 gene expression was observed only after exposure to 5 µM IPP, while vmat2 mRNA expression was significantly downregulated after exposure to 5 µM IPP, TMPP and TPHP. BDE-47 did not affect the expression of either transporter at DIV21. Metabolomics analysis showed a significant decrease in the neurotransmitter dopamine after exposure to all OPFRs, but not BDE-47. Again, exposure to IDDP and TMPP had the strongest effect, as the decrease was observed already at DIV14 (Fig. 4b and suppl. material 1). In conclusion, BDE-47 showed very different effects to the other FR increasing not decreasing vmat1and no effect on vmat2 expression.
In summary, genes involved in the enzymatic transformation of glutamate to GABA were affected by all FRs at non-cytotoxic concentrations, again with the strongest effects observed by exposure to OPFRs. This indicates that the ratio of these neurotransmitters might be altered. In fact, the neurotransmitters glutamate and GABA were decreased as well as α-ketoglutaric acid (glutamate derivate). In addition, the expression of the dopaminergic transporter dat and dopamine levels were decreased for all OPFRs. Less prominent effects were observed on genes encoding the vesicular Fig. 4 Heatmap (double gradient, green-minus; red-plus) illustrating a genes (measured with RT-qPCR) and b neurotransmitters (metabolites) (measured with Q-TOF LC-MS), involved in neurotransmitter transportation and production after exposure to FR. All data were normalized to mean of untreated control cells (0) and are displayed as means log 2 -fold change in each independent experiment, from at least three independent experiments performed in 2-3 replicates (genes) or at least two independent experiments performed in 4-6 replicates (metabolites) transporters and only at the highest concentrations. Minimal effects were observed after exposure to BDE-47 on selected genes.

Exposure to FR affects the gene expression of glial markers
One advantage of the 3D rat primary neural cell cultures is the presence of different cell types of the brain including glial cells, (astrocytes, oligodendrocytes and microglia). Aiming to assess glial toxicity due to exposure to FR, the gene expression of two specific markers for mature astrocytes were evaluated: glial fibrillary acidic protein (gfap) and the calcium-zinc-binding protein s100 beta (s100β).
The expression of gfap was upregulated after the exposure to all OPFR (5 µM) at DIV14 (Fig. 5a and suppl. material 1). However, this effect was lower or even reversed at DIV21. OPFR exposure altered the mRNA Fig. 5 Heatmap (double gradient, green-minus; red-plus) illustrating a genes (measured with RT-qPCR) and b metabolites (measured with Q-TOF LC-MS), identified in glial cells (astrocytes, oligodendrocytes and microglia) after exposure to FR. All data were normalized to mean of untreated control cells (0) and are displayed as means log 2 -fold change in each independent experiment, from at least three independent experiments performed in 2-3 replicates (genes) or at least two independent experiments performed in 4-6 replicates (metabolites) levels of the astrocytic marker s100β at a later stage with significant increase at DIV21. On the contrary, BDE-47 decreased gfap expression at 21 DIV. This might indicate an astrogliosis induced by the OPFR not seen for BDE-47.
The effect of FR on the initial stages of the CNS developmental process, including proliferation of progenitor cells, was assessed by the expression of nestin. This protein is a cytoskeletal intermediate filament that is gradually replaced by cell-specific intermediate filaments such as nf-68 and nf-200 in the neural cells and gfap in astrocytes during the course of the CNS differentiation and maturation (Jin et al. 2009;Wiese et al. 2004). Furthermore, nestin can be re-expressed in activated astrocytes in the event of brain or neuronal injury and for this reason it is also recognized as a sensitive marker for reactive astrocytes (Brook et al. 1999;Chen et al. 2002;Sergent-Tanguy et al. 2006).
Nestin expression at DIV14 was significantly upregulated at the highest concentration (5 µM) and increased even further at DIV21 after exposure to IPP, IDDP and TMPP and was not affected by the exposure to BDE-47 and TPHP at DIV21 (Fig. 5a and suppl. material 1). This is concordant with astrogliosis induced by the OPFR, but not BDE-47.
Lactic acid is mainly produced in astrocytes in the CNS and is an essential element of neuron-glia metabolic interactions (Pellerin 2003). We observed a decrease in lactic acid after exposure to IPP, IDDP and TMPP ( Fig. 5b and suppl. material 1). Creatine is a metabolite that has shown to be a marker of gliosis when increased (Konaka et al. 2003) as well as a marker of neurotoxicity when decreased (Diserens et al. 2018;van Vliet et al. 2008). Exposure to IPP, IDDP and TMPP significantly decreased the levels of creatine ( Fig. 5b and suppl. material 1). Moreover, a significant decrease in dl-serine was observed after exposure to IPP, IDDP and TMPP with stronger effects to IDDP and TMPP ( Fig. 5b and suppl. material 1). D-serine is considered an astroglia-derived neurotransmitter (Van Horn et al. 2013), though it is also produced marginally in neurons. However, with the method chosen here for mass spectrometry, it is not possible to distinguish between d-and l-isoforms, therefore, they are reported together.
In summary, expression of glial markers indicates that exposure to OPFR may activate astrocytes, either in a primary or secondary fashion, due to neuronal damage while less effect was observed after exposure to BDE-47 on selected genes and metabolites.
Myelin basic protein (mbp) is a protein expressed in mature oligodendrocytes and is a well-recognized oligodendrocytic marker. The highest concentration (5 µM) of TPHP, IPP and TMPP significantly downregulated mbp expression at both DIV14 and DIV21 (Fig. 5a and suppl. material 1). Exposure to BDE-47 and IDDP did not interfere with the expression of mbp mRNA. In conclusion, at the higher concentration three of the OPFR but not BDE-47 might interfere with myelination.
To further analyze the possible effects on myelination, myelin-associated glycoprotein (mag) was assessed. The mag is a transmembrane glycoprotein essential for the formation and maintenance of the myelin sheath by promoting glia-axon interactions (Paivalainen and Heape 2007;Quarles 2007). At DIV14, the expression of mag was upregulated by exposure to all FR (5 µM) except IPP (Fig. 5a and suppl.  material 1). At DIV21, the expression was instead significantly downregulated by IPP and TMPP (5 µM) and back to control levels in BDE-47, TPHP-and IDDP-treated samples. This indicates some possible restored effects on myelination by these FR.

Exposure to FR may trigger an inflammatory response
The inflammatory cytokines interleukin-6 (IL-6) and tumor necrosis factor-alpha (TNF-α) are synthesized by cells of the CNS, including neurons and glial cells, to exert neuroprotective roles (Carlson et al. 1999). A significant increase in the expression of il-6 was observed after exposure to the highest concentration (5 µM) of IPP, IDDP and TMPP at DIV21 (Fig. 5A and suppl. material 1). No change was observed after exposure to BDE-47 and TPHP. Expression of tnf-α was significantly increased only after exposure to 5 µM IPP at DIV21 (Fig. 5a and suppl. material 1).
To characterize further the inflammatory component, the macrophage colony-stimulating factor 1 receptor (csf1r), involved in microglial proliferation and activation (Mitrasinovic et al. 2005), and the allograft inflammatory factor (aif), upregulated in activated microglia due to inflammation (Deininger et al. 2002), were evaluated. The expression of csf1r and aif was only increased after exposure to TMPP ( Fig. 5a and suppl. material 1). In contrast, levels of the metabolite 2,3-pyridinedicarboxylic acid (quinolinic acid) produced by activated microglia was significantly decreased after exposure to all OPFRs ( Fig. 5b and suppl. material 1).
The increase of cytokine gene expression indicates that exposure to OPFRs may induce an inflammatory response in the rat primary neural model, but the microglia are only activated after exposure to TMPP. The latter was also supported by the metabolite observation. However, to fully confirm the lack of activated microglia, additional time course studies need to be performed.

Exposure to IPP modified the transcriptome
To get a more complete representation and to identify additional genes and pathways of interest, the whole transcriptome was analyzed after exposure to 1 µM IPP at 21DIV versus correspondent solvent control. We observed 1971 significantly different expressed genes (|FC|> 1.4 and p < 0.05 false discovery rate Benjamini-Hochberg procedure), which were further analyzed by pathway/category overrepresentation analysis using KEGG, WikiPathways and Reactome (Fabregat et al. 2018;Kanehisa and Goto 2000). The top perturbed pathways were associated with four main mechanisms by overrepresentation gene analyses (Table 3): (1) downregulation of neurotransmitter receptors including VMAT1 and associated intracellular signal transduction (G-protein, Ca 2+ , MAPK); (2) upregulation of immune response, inflammation and defense response mainly driven by MHC-I, FasL, Fas, complement system, FC IgG receptors (macrophages) and AIF (specific microglial activation biomarker); (3) upregulation of cell cycle; (4) changed fatty acid metabolism and transportation.
The transcriptome and RT-PCR analysis indicates similar effects.
In conclusion, after gene-wise careful evaluation, most pathways/categories could be classified into four categories: (1) transmission of action potential, cell-cell signaling, synaptic transmission, receptor signaling, (2) immune response, inflammation, defense response, (3) cell cycle and (4) lipids metabolism and transportation.

Discussion
The U.S. Consumer Product Safety Commission (CPSC) accepted a petition to ban furniture, children's products, electronic enclosures, and mattresses containing any member of the class of organohalogen flame retardants in 2017. In response to the petition, the National Academies of Sciences (NAS) released a report in May 2019, titled "Scoping Report for Conducting a Hazard Assessment of Organohalogen Flame Retardants as a Class" (National Academies of Sciences 2019). The authors concluded that there is a need to categorize flame retardants into broad classes for risk assessment rather than regulate on individual compounds.
Following the voluntary phase-out of the PBDEs, there has already been an increase in the use of the OPFR as evidenced by their ubiquitous detection in air and in dust in a variety of indoor environments (Bergh et al. 2011;Hoffman et al. 2015;Stapleton et al. 2009). A recent study showed the presence of OPFRs in children from residential exposures (Phillips et al. 2018). OPFR metabolites have also been detected in urine in the general population (Hammel et al. 2016;Hoffman et al. 2017;Ospina et al. 2018). Considering the recent petition on the ban on organohalogens, the use of OPFR is further projected to be on the rise. However, their toxicological hazard has not yet been well characterized. Several recent streams of data are converging on the conclusion that OPFRs show equivalent (or in some cases greater) toxicity compared to some of the phased-out PBDEs in several human-derived cell-based models . These compounds have also been shown to produce neurobehavioral deficits in complementary animal models (Alzualde et al. 2018;Bailey and Levin 2015;Glazer et al. 2018;Oliveri et al. 2015;Yan et al. 2017;Zhang et al. 2019) and might be a regretful substitution (Zimmerman and Anastas 2015).
The current study sheds light on possible mechanisms by which OPFR may exert their developmental neurotoxicity using a 3D primary rat neural culture. The aggregates were exposed to non-or low-cytotoxic concentrations of BDE-47, TPHP, IPP, IDDP and TMPP (0.1, 1 and 5 µM) to understand the toxicity mechanism and reveal potential perturbations to crucial neurodevelopmental processes, including neuronal and glial differentiation, and functional maturation. Metabolomics and gene expression were evaluated in control and FR-treated cultures at DIV14 and DIV21, after 7 and 14 days of exposure, respectively.
Exposure to FRs altered several markers involved in neuronal morphology and function, both at the gene and metabolome level. A decrease was observed in the mRNA levels of nf-200, an essential component of the axonal cytoskeleton. The synthesis of nf proteins and the timely suppression of nf gene expression are of functional importance for neuronal activity during the differentiation and maturation process of the CNS. Hence, nf-200, a late developing cytoskeleton protein, appears as a crucial element for neuronal function, since it is involved in the regulation of the axonal growth and is believed to provide nf stability and resistance to protein breakdown by retarding the slower axonal transport component (Goldstein et al. 1987;Liu et al. 1994). The decrease in the expression of neurofilaments after exposure to FRs indicates a structural disruption that can disturb the axonal organization and may cause degeneration of the axons (Perrone Capano et al. 2001) and the onset of neuropathology. Decreased expression of genes involved in cytoskeleton organization in neurite formation alongside with altered locomotor behavior has also been observed in zebrafish larvae after exposure to OPFR (Sun et al. 2016). Moreover, exposure to OPFR has previously been shown to affect neurite outgrowth at similar concentrations as in this study Hausherr et al. 2014) supporting our results. Neurite outgrowth is one of the most established in vitro developmental neurotoxicity tests  covering one of the key processes of CNS development that if perturbed likely leads to adverse outcomes.
Moreover, genes involved in neuronal function such as synaptogenesis and receptor expression were altered after exposure to all OPFR. Synapsin1 (syn1) has a functional Table 3 Whole transcriptome analyzes using KEGG, WikiPathways and Reactome after exposure to 1 µM IPP at 21DIV revealed alternations in genes belonging to four main mechanism, (1) neurotransmit-ters, (2) immune response, cell cycle and (4)  role during neuronal development and playing an important part in the formation of synapses and in the regulation of neurotransmitter release by control of the amount of synaptic vesicles ready for exocytosis at the axon terminal (Evergren et al. 2007;Ferreira and Rapoport 2002). The decrease in the expression of syn1 after exposure to OPFRs could impact the neuronal function and be linked to a defective neuron signaling system or even death. Previously, the expression of another synaptogenesis marker syn2a has been shown to be downregulated after exposure to OPFRs in zebrafish larvae (Sun et al. 2016) supporting disruption of neuronal signaling. Indeed, exposure to all FR used in this study has previously shown to decrease the neural network activity in rat primary cortical cultures acutely . Exposure to IPP and BDE-47 during development also showed decrease in activity in the network formation assay (NFA) (Shafer et al. 2019). However, exposure to TMPP, IDDP and TPHP did not induce any changes. The decreased activity (~ 6-16 µM) was reported close to the highest concentrations tested in the NFA (20 µM) and it is possible that a change in activity would be observed in the NFA if higher concentration had been tested. The decrease in syn1 in the rat 3D model is clearly at lower concentrations (1 µM). The window of development, cell populations, exposure scenario and culture condition (2D vs. 3D) differ between these studies and could contribute to the vulnerability. Additional experiments need to be performed to understand if the effects observed in the rat 3D brainsphere model also leads to functional changes. Still, the decrease in levels of the metabolites NAA and L-aspartic acid observed in this study is a clear indication of neuronal damage and impaired neuronal function. NAA is a common biomarker used in patients to diagnose stroke and neurodegenerative disorders (Alakkas et al. 2019;Baslow et al. 2003) and is produced from l-aspartic acid and acetyl-coenzyme A (Ariyannur et al. 2010).
In the rat brainsphere model, we also observed that neurotransmitters and mRNA levels of their receptors were altered due to exposure to the selected FR. The strongest effect was noted in the expression of the glutamate NMDA receptor, where subunits grin1 and grin2a were significantly downregulated at non-cytotoxic concentrations, especially after exposure to IPP and TMPP. The switch in subunit 2b-2a during development is crucial for proper maturation of the brain, as it is involved in neuronal function and synaptogenesis (Liu et al. 2004;Luthi et al. 2001). Interestingly, the expression of subunit grin2c was significantly increased after FR exposure and could be a sign of altered receptor affinity to the neurotransmitters (Yi et al. 2020). It was previously observed that exposure to TMPP reduces the response to glutamate in mouse cortical neurons (Hausherr et al. 2014). In this study, we observed a clear decrease in glutamate (l-glutamic acid) and the glutamate derivate α-ketoglutaric acid after the exposure to FR, which supports such theory. This was already observed at the lowest concentration (0.1 µM) and at the early time point (DIV14), indicating that this could be one of the mechanisms of the tested OPFR. Moreover, it is well documented that the NMDA receptor is a potential target during brain development and its inhibition can lead to decreased synaptogenesis, decreased neuronal network formation and ultimately impairment of learning and memory abilities (AOP 12, https ://aopwi ki.org/aops/12 and AOP 13, https ://aopwi ki.org/aops/13) (Sachana et al. 2018;Spinu et al. 2019;Wang et al. 2018). The altered gene expression of subunits of the NMDA receptor and the reduction in glutamate are therefore of high concern and of importance to investigate further.
GABA is the principal inhibiting neurotransmitter in the mature brain, but also exerts excitatory actions during the formation of the CNS (Ben-Ari et al. 2007;Li and Xu 2008). It is implicated in the proliferation of neural progenitor cells (Haydar et al. 2000), migration (Luhmann et al. 2015), differentiation (Barbin et al. 1993;Ganguly et al. 2001), outgrowth of neurites (Maric et al. 2001), and synaptogenesis (Ben-Ari 2002). The neurotransmitter GABA was decreased after exposure to IPP, IDDP and TMPP alongside with decreased expression of GABA A receptor subunit alpha1 (gabra1), but to a lower degree than glutamate and its receptor and mainly after prolonged exposure (14 days). Moreover, two isoforms of glutamate decarboxylase (gad1 and gad2), which are the enzymes responsible for catalyzing GABA synthesis from glutamate, were significantly decreased. This likely contributes to the decrease seen in GABA levels, but can also be a secondary effect due to the decrease in glutamate levels. Inhibition of the GABA A receptor has previously been suggested to contribute to the toxicity of TMPP (Gant et al. 1987) and TPHP (Flaskos et al. 1998;Gant et al. 1987). Reduction of GABAergic neurons is one key event identified in AOPs (AOP 10, https ://aopwi ki.org/aops/10 and AOP 54 https ://aopwi ki.org/aops/54) for neurotoxicity and developmental neurotoxicity (Bal-Price et al. 2015b;Spinu et al. 2019;Westerholz et al. 2010) and could contribute to the toxicity of OPFR. Moreover, the ratio of glutamate/GABA has shown to be disturbed in autistic children (Gaetz et al. 2014;Gogolla et al. 2009;Rippon et al. Rubenstein and Merzenich 2003), raising a concern for OPFR as replacements, especially in children products, particularly as exposure to organophosphate pesticides and PBDEs is a suggested environmental risk factor for autism (Mostafalou and Abdollahi 2018;von Ehrenstein et al. 2019;Vuong et al. 2018;Ye et al. 2017). Furthermore, decreases in the neurotransmitter dopamine and genes related to dopamine transportation were observed. It is known that dopaminergic neurons are more vulnerable to oxidative stress, which has been identified as potential mechanisms in the toxicity of some FR (Hendriks et al. 2014;Pellacani et al. 2014;Tagliaferri et al. 2010;Wu et al. 2016). Genes involved in oxidative stress were not altered in our transcriptomics data (data not shown), but the time points of sample collection might not be optimal to observe these effects.
Astrocytes are considered key participants in the brain maturation for providing neuronal structural, trophic and metabolic support, synthesis of growth factors and defense mechanisms, and for influencing synapse formation (Barker and Ullian 2008). Thus, the adverse effect of chemicals on astrocytes may interfere with the morphological development and functional performance of neurons in the CNS. Aiming to assess glial toxicity due to FR exposure, the gene expression of two specific markers for mature astrocytes were evaluated: glial fibrillary acidic protein (gfap) and the calcium-zinc-binding protein s100 beta (s100β). Both genes were significantly upregulated after exposure to all OPFRs. Gfap is necessary for many important processes in the CNS correlated with neuronal survival (Liedtke et al. 1996;Tardy et al. 1990). Increased gfap is an indication of gliosis/activated astrocytes, a common response of glial cells to neuronal injury and implicated in several neurological disorders such as schizophrenia, bipolar disorder and depression (Johnston-Wilson et al. 2000). In the developing CNS, s100β is believed to give support to growth, survival and differentiation of neurons (Wang and Bordey 2008), but also to play an important role in the recovery of the CSN after injury (Yardan et al. 2011). It is used as a biomarker of brain damage and as a parameter of activation of astrocytes (Esposito et al. 2006). The elevated level of s100β may also by itself produce adverse effects, including overgrowth of dystrophic neurites (Griffin et al. 1995), and is related to the occurrence of various neuropathologies (Yardan et al. 2011). The data from a study on zebrafish exposed to OPFRs are in concordance with our results by showing an increase in gfap expression (Sun et al. 2016).
In addition, the gene expression of the neural precursor marker nestin was significantly upregulated after the exposure to OPFRs. In the normal course of neurodevelopment, astrocytes and neurons differentiate from neural precursor cells leading to a decrease in expression of nestin over time (Shaltouki et al. 2013). However, upon injury to the CNS, it is transiently re-expressed in activated astrocytes (Brook et al. 1999;Chen et al. 2002;Michalczyk and Ziman 2005) and is therefore recognized as a sensitive marker for activated astrocytes (Hogberg et al. , 2009. The upregulation of nestin in reactive astroglia has been observed in several diseases, e.g., cerebral ischemia, hippocampal excitotoxicity lesions, traumatic brain injury, and after MPTP (1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine) exposure Wei et al. 2002). It is likely that the increase in nestin expression observed in this study is due to reactivated astrocytes as also other glia markers were increased.
As a response to injury or chemical insult, glia cells (astrocytes and microglia) become activated and begin to produce different proinflammatory and neurotoxic substances such as cytokines and free radicals (Bal-Price and Brown 2001). TNF-α is a cytokine considered to be a primary proinflammatory mediator capable of playing a dual functional role by promoting tissue regeneration/growth and destruction (Wajant et al. 2003). IL-6 is another factor that exerts distinct functions on the CNS by participating in inflammatory responses and infections, and modulating neural processes (Abreu et al. 2018;Scheller et al. 2011). It has been shown that IL-6 influences the differentiation of neurons and astrocytes (Oh et al. 2010), but it can also be neurotoxic and cause neuronal death (Brown and Bal-Price 2003;Conroy et al. 2004). The overexpression of IL-6 has been linked to the onset of neurodevelopmental and neurodegenerative diseases and mental disorders such as schizophrenia and autism (Conroy et al. 2004;Smith et al. 2007;Wei et al. 2012). We observed a significant upregulation of the il-6 mRNA after exposure to most OPFR at DIV21. After exposure to IPP, the gene expression of tnf-α was also increased. This indicates the onset of an inflammatory response in the 3D model after exposure to OPFR, possible as a secondary effect due to neuronal injury. Only exposure to TMPP led to possible microglia proliferation and activation as shown by upregulation in the csf1r and aif1 gene expression. It should be noted that the experimental setup is not ideal to detect an inflammatory response, as cytokine release is rapid within hours after a trigger and our focus was on prolonged effects on development after long-term exposure. In addition, the microglia population, that is mainly responsible for cytokine release, is minor in the model; thus the cytokines levels can be close to the detection limits. However, transcriptomics data for IPP showed a clear increase in genes involved in inflammatory pathways and implies that the immune system plays a role in the the toxicity of OPFR. Moreover, AOP 13 (https ://aopwi ki.org/aops/13) links NMDA receptor inhibition during brain development to neuroinflammation and impairment of learning and memory (Villeneuve et al. 2018). Critical effects observed in this study were decreased gene expression of the NMDA receptor and glutamate levels that can be linked to the observed neuroinflammation. This is the first time that potential underlying pathway(s) associated with developmental neurotoxicity of these OPFR have been investigated by utilizing a threedimensional rat primary neural model. Although there are some limitations with the cell model including but not limited to kinetics, metabolism, being non-human, nonetheless it is valuable in showing how alternate streams may aid toward mechanistic understanding of classes of compounds. The advent of human brainspheres (Pamies et al. 2017) now enables similar studies in human cell models (Pamies et al. 2018;Zhong et al. 2020) work on FR is ongoing. This model also captures several end points that are not currently evaluated in traditional DNT guideline studies including, but not limited to the role of neurotransmitters and critical receptors that have previously been implicated in neurodevelopment (e.g. glutamate and GABA). Importantly, using in vitro to in vivo extrapolation, these findings suggest that activity is noted at relevant human exposures (within model constraints) (Blum et al. 2019). Studies are underway at the National Toxicology Program to evaluate a couple of these compounds using traditional guideline DNT studies to help provide in vivo anchors. Nonetheless, it is not feasible or practical to perform guideline studies on every member of a class of compounds. Hence, approaches such as these in combination with other assays using a DNT battery (Aschner et al. 2017;Bal-Price et al. 2018;Behl et al. 2015;Fritsche et al. 2017) complement traditional animal testing by providing guidance on prioritization, and shedding light on possible mechanistic understanding which may contribute to putative AOPs, which is currently not a part of standard guideline DNT testing. Acknowledgement This work is based on model optimization supported by FDA (U01FD004230) "DNTox-21c-Identification of pathways of developmental neurotoxicity for high throughput testing by metabolomics". Some of the transcriptomics work was supported by the NIH transformative research grant "Mapping the Human Toxome by Systems Toxicology" (RO1ES020750).
Rita de Cássia da Silveira e Sá was supported by CNPq -Conselho Nacional de Desenvolvimento Cientifico e Technologico, of the Ministry of Science, Technology and Innovation of Brazil, Department of Physiology and Pathology, Federal University of Paraíba, João Pessoa, Brazil. Ozge Cemiloglu Ulker was supported by The Scientific and Technological Research Council of Turkey and Ankara University Faculty of Pharmacy.
We are thankful to the Division of National Toxicology Program, U.S. National Institute of Environmental Health Sciences (NTP, NIEHS, Research Triangle Park, NC) for providing the FR. We acknowledge Dr. Kristen Ryan, NIH, NIEHS for NIH internal review of the manuscript.

Compliance with ethical standards
Conflict of interest TH and HH are named inventors on a patent by Johns Hopkins University on the production of human mini-brains, which is licensed to AxoSim, New Orleans, LA, USA. They consult AxoSim and TH is shareholder.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.