Impacts of Gestational FireMaster 550 Exposure on the Neonatal Cortex Are Sex Specific and Largely Attributable to the Organophosphate Esters

Chemical flame retardants (FRs) are applied to a myriad of consumer products to comply with flammability standards aimed at delaying the spread of fire. Due to continued widespread use, FRs are prevalent pollutants in the natural and human environment including house dust and food sources [1‒3]. This is concerning because multiple FRs are thought to be neurotoxic and endocrine disrupting [4‒6]. In response to these toxicological concerns, over the last 2 decades, the FR landscape has changed considerably with some phasing out and others phasing in. Most significantly, use of some polybrominated diphenyl ethers (PBDEs) was discontinued after being identified as thyroid disrupting and linked to deficits in memory and cognitive function [7]. Since the PBDE phaseout, different brominated FRs (BFRs), sometimes called novel BFRs and other alternative FRs, have been developed and implemented as replacements [8, 9]. Also, the use of organophosphate esters (OPEs) as FRs (OPFRs) has rapidly grown in popularity largely because of their relatively rapid metabolism and assumed lower human toxicity, although that presumption of greater safety has been questioned [10‒13]. Notably, some OPEs have also been used in a variety of industrial applications for decades prior to their adoption as FRs, so global production was already high and human exposure widespread [14]. The newer BFRs and OPFRs present in commercial mixtures such as FireMaster 550 (FM 550) were designed to be less neurotoxic than their predecessors, yet their potential to disrupt neurodevelopment has not been thoroughly investigated. Here, we utilized a combination of transcriptomic and lipidomic approaches to understand how exposure to these “next generation” FRs impacts neurodevelopment in Wistar rats of both sexes.

Commercially available FRs are often mixtures, sometimes containing both BFR and OPFR components. FM 550 is one such mixture that has been used since the early 2000s as an additive to treat foam-based furniture, baby products, and car seats, among other common items [15, 16]. FM 550 is composed of two brominated compounds, 2-ethylhexyl-2,3,4,5-tetrabromobenzoate (EH-TBB) and bis(2-ethylhexyl) 2,3,4,5-tetrabromophthalate (BEH-TEBP), the OPE triphenyl phosphate (TPHP), and numerous isopropylated triarylphosphate isomers (ITPs) [17, 18]. FM 550 components can readily escape from products and multiple studies have demonstrated the environmental mobility, persistence, and bioaccumulation [2, 19]. Indoors, this includes contamination of house dust [1, 18, 20], where inhalation or ingestion is the most significant route of human exposure, especially for children [14]. Higher levels of FM 550 components are regularly found in children compared to adults, likely due to increased hand-mouth behavior [21‒23]. In humans, OPFRs are metabolized more quickly than BFRs, leading some to conclude they are, consequently, less likely to result in fetal exposure or toxicity [24, 25]. We and others, however, have detected them in animal and human placentas, along with evidence of placental stress and other biomarkers of disrupted placental function [26‒30].

Evidence of OPFR neurotoxicity is also growing [10, 11, 31, 32]. Of additional concern is that human OPFR exposure is now on par with, or even higher than, historical PBDE peak exposure levels [10]. Although the OPEs used as FRs were designed to have limited capacity to inhibit acetylcholine esterase, the primary mode of action by which the structurally similar OP pesticides have long been recognized to be neurotoxic, emerging evidence suggests they can still disrupt cholinergic and other neurotransmitter signaling pathways including cholinergic neurodifferentiation [12, 32, 33]. Despite their widespread use, proper risk assessment of OPFRs and the newer BFRs has yet to be conducted.

To date, we and others have shown in a variety of experimental models that FM 550 is disruptive of sexually dimorphic socioemotional behaviors [26, 34, 35] and their related neural pathways [26, 28, 29, 31, 35‒41]. Less is known, however, about how each of the chemical classes may contribute to those adverse outcomes. In a related prior study using Wistar rats, we explored how perinatal exposure to FM 550 or its component classes alters a suite of behaviors in both sexes as adults [35]. The animals were gestationally exposed to 1,000 µg of the BFR mixture, 1,000 µg of the OPFR mixture, or 2,000 µg FM 550 [35] via oral intake by the dam. We demonstrated that both the BFRs and the OPFRs can induce adverse behavioral outcomes, albeit differently, suggesting each class impacts the brain via different modes of action.

The present studies were conducted to further understand the underlying mechanisms driving these differences and to identify unique outcomes of exposure to the full FM 550 mixture. Based on previous work in rodents demonstrating dysregulated neurodevelopment of the fetal forebrain [28, 29] and behavioral deficits [26, 35] following gestational exposure to FM 550 and OPFRs, we chose to focus on the prefrontal cortex. Using postnatal day (PND) 1 cortices from siblings of the animals used for our previously published behavioral study [35], we employed a combination of transcriptomics, lipidomics, and validating PCR approaches to determine the possible pathways each chemical class disrupts. While transcriptomics is commonly utilized in exposure studies, lipidomics studies and combined multi-omic evaluations are rare due to the complexity of lipidomics measurements (e.g., numerous lipid isomers and large dynamic range). However, since lipids constitute a major structural component of brain tissue and are involved in numerous signaling and biochemical processes, assessing how specific lipid species change in response to exposure is a potentially powerful method of assessing developmental neurotoxicity. Here we combined new analytical techniques to provide an in-depth lipid analysis of the species affected by environmental exposure in the developing brain [42]. Collectively, our analyses explored how FM 550 and its component classes altered cortical lipid composition in combination with transcriptomics, to better understand molecular and structural perturbations associated with behavioral and cognitive deficits observed in prior studies.

In addition to its disruptive effects on sexually dimorphic socioemotional behaviors, FM 550 and its components have been found to be endocrine disrupting [43], disruptive of embryonic neural pathways [28, 29], and adipogenic with compromised bone composition as a potentially consequential outcome [44, 45]. Thus, we anticipated significantly altered pathways and genes, including peroxisome proliferator-activated receptors (PPARs) and their transcriptional partners, associated with these phenotypes. Given the rapidly emerging body of work demonstrating that OPFRs may not be as low risk as previously thought [11], this component class was of key interest.

The ARRIVE (Animal Research: Reporting of in vivo Experiments) Guidelines “Essential 10” Checklist for Reporting Animal Research was used in the construction of this manuscript with all elements met [46]. The ARRIVE guidelines were developed in consultation with the scientific community as part of an NC3Rs (National Centre for the Replacement, Refinement, and Reduction of Animals in Research) initiative to improve the standard of reporting of research using animals.

All animals were siblings of animals used for a previously published study [35]. Animal care, maintenance, and experimental protocols met the standards of the Animal Welfare Act and the US Department of Health and Human Services “Guide for the Care and Use of Laboratory Animals” and were approved by the North Carolina State University (NCSU) Institutional Animal Care and Use Committee (IACUC; approval number 20–150). All procedures were approved and monitored by a supervising veterinarian throughout the duration of the project. Wistar rats were obtained from Charles River (Raleigh, NC) and/or bred in-house as indicated in humidity- and temperature-controlled rooms, each with 12-h:12-h light:dark cycles at 25°C and 45–60% average humidity in the AAALAC-Approved Biological Resource Facility at NCSU. As in our prior studies [47, 48], and in accordance with recommended practices for endocrine disrupting chemical (EDC) research, all animals were housed in conditions specifically designed to minimize unintended EDC exposure including use of glass water bottles with metal sippers, soy-free diet, woodchip bedding, and thoroughly washed polysulfone caging [49‒51]. Animals were generated from breeding adult naive Wistar rats obtained from Charles River (56 females at PND 72 and 16 males at PND 90) and maintained on Teklad 2020 (phytoestrogen-free) diet. Dams were weighed and randomly assigned to the four treatment groups (control, BFR, OPFR, and FM 550) with an average weight per group equaling approximately 300 g (mean and standard deviation for each group: control 295 g ± 16; BFR 301 g ± 13; OPFR 293 g ± 14; FM 550 299 g ± 12).

All animals were orally dosed using concentrated solutions of FM 550, BFR or OPFR as previously described [35]. Exposure concentrations were based on previous work demonstrating transplacental/lactational transfer, bioaccumulation, and neurotoxicity at levels lower than the purported NOAEL of 50 mg/kg for the BFR mixture [26, 27, 29, 43]. A concentration of 1,000 µg/day was chosen for the BFR and OPFR groups to represent the approximate concentration of BFRs and OPFRs in the 2,000 µg/day FM 550 dose [18, 52]. All mixtures were prepared in Dr. Heather Stapleton’s laboratory at Duke University. Dosing was not done blinded because the FM 550 mixture was a larger volume, and the control animals were always dosed first to minimize the risk of cross-contamination. The individual who did the dosing was not the same as the individuals conducting the experiments to ensure experimental blinding. The FM 550 commercial mixture was obtained from Great Lakes Chemical (West Lafayette, IN) and a 50 mg/mL dosing solution was prepared by diluting the appropriate amount of FM 550 into sesame oil. The Stapleton laboratory also prepared the BFR and OPFR mixtures into working 50 mg/mL solutions. The OPFR mixture is the same as previously reported [29, 52] and includes organic components of TPHP and isopropylated triphenyl phosphates (ITPs; also abbreviated IPPs). The FM 550 brominated component mixture (BFR) contained 2-ethylhexyl-2,3,4,5-tetrabromobenzoate (EH-TBB) and bis (2-ethylhexyl)-2,3,4,5-tetrabromophthalate (BEH-TEBP). Each dosing solution (sesame oil vehicle, 1,000 µg/day BFR, 1,000 µg/day OPFR, and 2,000 µg/day FM 550) was stored in foil-wrapped scintillation vials at 4°C.

Dam exposure was once per day, orally via a food treat as previously described [35], beginning 72 h after pairing and continuing through PND 1. All animals were presumed pregnant by 72 h. The dams gave birth within 3 days of each other, an observation that supports this presumption. The oral dosing method was chosen to reduce handling and, consequently, possibly confounding prenatal stress [53]. As previously reported [28], dosing was based on average group body weights (bw) of 300 g before they were impregnated, producing exposures of approximately 3.3 mg/kg bw/day BFR or OPFR and 6.6 mg/kg bw per day FM 550. Relative exposure likely decreased slightly across gestation as body weight increased [53]. We elected not to adjust for weight gain across gestation to minimize animal handling. We concluded that the benefits of making small adjustments to dosing to account for increasing dam weight would not outweigh the risk of prenatal stress daily weighing could induce.

Parturition checks were done in the morning at the time of dosing and the first day a litter was present was designated PND 0. On PND 1, pups were euthanized via rapid decapitation 4 h after dosing (13:00 h ± 60 min) to control for time post-exposure and time of collection. Whole neonatal heads were collected, flash frozen, and stored at −80°C. To verify offspring sex, a single paw was collected, and PCR was performed for the sry gene as previously described [26, 54]. Once sex was verified, 6 male and 6 female PND 1 whole heads per exposure group were selected. For all compounds, only one pup per sex (littermates in most cases) per litter was used (litter is the experimental unit). The OPFR group only had 5 females because of limited litter availability.

Similar to our previous work [55, 56], the whole heads were first cryosectioned at 20 µm (Leica CM 1900) from the most anterior portion of the cortex and the forebrain portion of the sequential sections were collected with a brush (Fig. 1), transferred to an Eppendorf tube, and stored at −80°C for RNA extraction and sequencing. Anatomical landmarks were identified using the Developing Rat Nervous System rat brain atlas [57]. For transcriptomics, forebrain tissue was collected from the regions depicted on P0 coronal plates 212–216. Anterior landmarks included clear separation of cortex and olfactory bulb and collection stopped on plate 216 as identified by the position of the corpus collosum, fornix, and lateral ventricle. A 1.25 mm wide by 1–1.15 mm deep micropunch was then used to isolate the medial portion of the cortex for lipidomics. These isolated tissues represent plates 216–221. Prior, pilot work established the smallest amount of tissue we could use for each approach and still obtain robust data. The transcriptomics and lipidomics tissues were taken sequentially to obtain sufficient tissue amounts while maximizing anatomical specificity to the greatest degree possible. Moreover, while transcription is subregion and even cell-specific (here we focused on the prefrontal region of the cortex), lipid content is more regionally homogeneous. Thus, the lipid composition observed likely also applies to other cortical regions, including the prefrontal cortex.

RNA extraction was performed with the Qiagen RNEasy Miniprep kit according to the manufacturer’s protocol (Cat. 74134; Qiagen). RNA quality was determined using an Agilent 2100 Bioanalyzer and all samples were found to have RNA integrity numbers (RIN) ≥9.2. Sequencing libraries were prepared as described previously [58] by the NC State Genome Science Laboratory (GSL), using NEBNext Ultra Directional RNA Library Prep kit and NEBNext Poly(A) mRNA Magnetic Isolation Module (catalogs E7420 and E7490; New England Biolabs, Ipswich, MA, USA) for Illumina sequencing. Isolation, heat fragmentation, and priming were performed according to manufacturer instructions. cDNA synthesis was followed by purification and size selection. Finally, library clean-up was performed using AMPure XP beads (Cat. A63881; Beckman Coulter Genomics, Brea, CA, USA) and quality was assessed using the Agilent 2100 Bioanalyzer. For each sample, library sequencing was performed using a 50-bp paired-end protocol on a single lane of an NovaSeq6000 sequencer. Approximately, 50–90 million uniquely mapped reads were generated per fetal forebrain library. The raw data were deposited in NCBI’s Gene Expression Omnibus (GEO Accession: GSE211430).

Quality control of read data was evaluated with FastQC. We used Cutadapt to trim sequences associated with library indexes, and then alignment was performed using the STAR short read aligner [59] to the Rattus norvegicus (rn6) reference genome. The number of reads mapped to GENCODE was determined using featureCounts software. PCA plots and heatmaps were generated using the 500 most variable genes to assess data quality. It was observed that a small number of disproportionally over-expressed genes were potentially driving 11 samples to be outliers. To further examine this phenomenon, a heatmap of the top 50 expressed genes was created (online suppl. Fig. S1a; for all online suppl. material, see www.karger.com/doi/10.1159/000526959). Additional details of each gene ID listed in the heatmap were generated using Uniport (www.uniport.org) including gene symbol, full gene name, and function. The NCBI database (www.ncbi.nlm.nih.gov) was then used to confirm the presence and/or absence of those genes in neuronal tissue using RPKM data (online suppl. Table S1a). The primary tissue sources of some of the 50 genes in question were blood, hair, muscle, the eye, skin, and bone, but not the brain. Because the forebrain tissue was collected from whole head cryosections using a brush, the risk of contamination by minute amounts of ectopic tissue, particularly bone and eye, from sticking to the brush or the brain tissue as it was isolated is plausible. Based on this analysis and in consultation with a bioinformatician, rather than excluding the 11 samples, 13 genes expressed in bone and 19 genes expressed in the eye were removed from the differentially expressed gene (DEG) set prior to data analysis (online suppl. Table S1b) to minimize analytical bias from ectopic expression. Additionally, because the PCA plots for the DEGs revealed clear separation by sex (online suppl. Fig. S1b), the data were analyzed within sex to account for sex-specific effects.

The goal of the approach was to probe for potential sex-specific vulnerabilities and directional (up- or downregulated) effects across exposure groups. Dispersion was estimated using the DESeq2 Bioconductor package in the R statistical computing environment to normalize count data, estimate dispersion, and fit a negative binomial model for each gene. The overall mean of the normalized counts, the log2 (fold-change), the p value, and the adjusted p value (padj) were used to identify DEGs. Pathway analysis was then performed on each exposure group’s up- and downregulated DEGs with padj <0.05. Pathway analysis was performed with the online tool Enrichr using the modules for GO biological processes, Predicted Protein-Protein Interaction (PPI) Hub Proteins, and KEGG.

To validate the RNAseq analysis and test predicted hub PPIs, a subfraction of each RNA extraction was also analyzed using the NanoString nCounter Analysis System (NanoString Technologies, Seattle, WA), which detects multiplexed gene targets without amplification or reverse transcription [60, 61]. Genes were selected based on the results of the pathway analyses including the PPI (see online suppl. Table S1 for selection details). In total, the code set included 29 up- or downregulated DEGs and 16 PPIs, as well as 6 positive and 8 negative controls and 7 housekeeping genes for quality assurance (online suppl. Table S2). Samples were normalized to 20 ng/µL prior to hybridization according to the manufacturer’s instructions. Briefly, a master mix was created using 70 µL of hybridization buffer into the reporter code-set tube. For each individual sample, 8 µL of master mix, 5 µL of RNA sample and 2 µL of capture probe-set was added to NanoString strip tubes and allowed to hybridize at 65°C for 20 h in a thermal cycler. Following hybridization, all samples were run on the NanoString nCounter Max Analysis System. The nCounter Prep Station incubation time was set to high sensitivity and the nCounter Digital Analyzer was set to max count.

Lipids were isolated from cortical micropunches using a modified Folch extraction [62, 63]. Cortical tissue was combined with 750 µL of −20°C methanol and homogenized in 2.0 mL, 2.4 mm tungsten-carbine bead tubes for 5 min with a Fisherbrand 24 bead mill. Samples were then transferred to Fisherbrand glass culture tubes containing another 750 µL of −20°C methanol, where 3 mL of chloroform and 200 μL of water were added. The samples were then vortexed for 30 s, sonicated for 30 min, vortexed again for 30 s, and incubated for 1 h at 4°C. Following incubation, 1.2 mL of water was added and samples were centrifuged for 10 min at 1,000×g. A 300 µL aliquot of the bottom lipid layer was finally transferred to a Sorenson microcentrifuge tube and dried in vacuo. The dried lipids were reconstituted in 190 µL of −20°C methanol and 10 µL of chloroform and then stored at −20°C until analysis with liquid chromatography, ion mobility spectrometry, collision induced dissociation, and mass spectrometry (LC-IMS-CID-MS).

The cortical lipid extracts were evaluated using an Agilent 1290 UPLC coupled to an Agilent 6560 IM-QTOF platform (Santa Clara, CA) with a commercial gas kit and MKS Instruments precision flow controller (Andover, MA). Lipids were first separated by reversed-phase LC by injecting 10 µL onto a Waters (Milford, MA) CSH column (3.0 mm × 150 mm × 1.7 µm particle size). The 34-min gradient having a flow rate of 250 µL/min and mobile phase A of 10 mM ammonium acetate in 40:60 LC-MS grade acetonitrile/water and mobile phase B of 10 mM ammonium acetate in 90:10 LC-MS grade isopropanol/acetonitrile is detailed in (online suppl. Table S3). Both positive and negative mode ESI analyses were performed on all samples. Following electrospray ionization (ESI), the ions were analyzed using the Agilent 6560 IM-QTOF MS platform [64, 65]. Collision cross sections were collected for all features detected and collision induced dissociation (CID) was performed with high purity nitrogen by ramping collision energies based on the ion arrival times analogous to previous IMS experiments [42, 66, 67]. Alternating scans of no fragmentation and all-ions data independent acquisition (DIA) were used to collect precursor and fragmentation information at 1 s/spectra for a mass range of 50–1,700 m/z. All lipidomic raw data for this analysis is deposited and available in MassIVE (MSV000090203).

Spectra were annotated in Skyline by matching features to an in-house lipid library of 778 lipids with experimentally determined LC, IMS, MS, and MS/MS information [68‒71]. Lipid annotations were made based on LC retention time, IMS collisional cross section, and m/z tolerances for precursor and fragment ions. The mass errors for all precursor annotations were <2 ppm and fragment mass errors were <10 ppm. LC elution times were within 2 s of predicted elution times derived from Skyline’s iRT feature, and experimental collisional cross sections were within 1% of the library values. Using our LC-IMS-CID-MS platform enables the annotation of lipid fatty acyl and head group moieties, but fatty acyl backbone attachment to the sn-1 or sn-2 positions, and double bond orientations and placement are not differentiable with this platform and method [72]. To elaborate, lipids were annotated with “_” to denote ambiguous fatty acyl positions and “/” when stereochemistry is known (e.g., PC[0:0_18:0] versus PC[0:0/18:0] or PC[18:0/0:0]). Lipids were also annotated with “a” or “b” to denote potential isomers. Additionally, features with multiple potential lipid matches are noted with “;” (e.g., PC[18:1_16:0]; PC[18:0_16:1]). Lipids were also assigned annotations for their summed carbon and double bond numbers when individual fatty acyl information could not be obtained or annotated (e.g., PE[34:1]). Finally, peak areas of all annotated lipids were exported from Skyline for further processing and statistical analysis.

The peak areas were transformed to a log₂ scale and normalized against their total ion current (TIC) since tissue variability was possible. Samples were screened for potential outliers through analysis of the PAV-RMD algorithm [73], Pearson’s correlation and principle component analysis using the pmartR [74] package in R [75] (Version 4.0.3 Vienna, Austria). Analysis of covariance (ANCOVA) was conducted to assess whether exposure influenced lipid abundances by comparing mean lipid peak areas of the control group to the exposed groups while accounting for differences in sex. Planned comparisons were implemented to compare control group lipid abundances to each of the exposed groups (BFR, OPFR, or FM 550) separately and determine whether each lipid had a statistically significant fold change. Holm’s correction was used to adjust p values for multiple comparisons, and lipids having an adjusted p ≤ 0.05 were considered statistically significant [76]. Lipids were visualized using the SCOPE cheminformatics toolbox [77]. In SCOPE, lipids were matched to the Lipid Maps Database to derive SMILES, from which the ECFP6 fingerprints were calculated using the rcdk R package [78‒81]. Hierarchical clustering of molecular fingerprints utilizing Euclidean distances and average linkages were calculated using the phangorn R package, and heatmaps were plotted with the ggtree R package [82, 83]. The resulting SCOPE dendrograms enabled visualization of the comparisons between the control and different exposure groups (BFR, OPFR, or FM 550).

Details of dam and offspring body weights, maternal care, and litter composition are reported in Witchey et al. [35] 2020. Briefly, dam body weights and maternal care did not differ across exposure groups at any time examined throughout the study. Although overall litter size did not differ across exposure groups, BFR dams had more female pups compared to controls. No significant differences were observed in offspring PND 1 body weights.

For all compounds, the number of DEGs was greatest in males (Fig. 2) with differential expression of 2,633 FM 550, 2,421 BFR, and 2,858 OPFR genes at an adjusted p value of <0.05. The affected numbers of DEGs were drastically fewer in female offspring: 7 FM 550, 2 BFR, and 535 OPFR DEGs. To perform the pathway analyses, data were examined within sex and direction of change (up- or downregulated). For the males, there were 1,406, 1,101, and 1,412 upregulated and 1,168, 1,276, and 1,375 downregulated genes in FM 550, BFR, and OPFR-exposed groups, respectively. Because so few genes were affected in the female FM 550 and BFR groups, pathway analysis was only performed on the OPFR female group, which consisted of 345 upregulated and 138 downregulated DEGs (Fig. 2).

Using the Ensembl IDs to compare DEGs across groups, separate Venn diagrams were created for upregulated and downregulated genes within sex (https://bioinfogp.cnb.csic.es/tools/venny/) to identify overlapping DEGs. In males, 35% of the upregulated (703) and 30.1% of downregulated (629) genes were identified as DEGs in all three exposure groups (Fig. 3). The FM 550 and BFR groups had 8.4% upregulated and 8.5% downregulated genes in common, while the FM 550 and OPFR groups had 13.7% upregulated and 6.8% downregulated genes in common. For the subsequent analysis, within each DEG subset (sex and direction), gene names were identified from the Ensembl IDs using the online tool Biodnet. Genes without names (unannotated regions) and genes for ribosomal proteins (rp) were removed from and excluded from the pathways analysis. The modules in En­richr used to identify enriched targets/pathways were transcription factor PPI, KEGG 2021 Human, and GO Biological Process 2021.

For the GO biological process analyses, in males, there were 157 FM 550, 163 BFR, and 180 OPFR upregulated and 24 FM 550, 42 BFR, and 55 OPFR downregulated processes. In males, approximately 36% [84] of the upregulated and 26% [18] of the downregulated GO biological processes were significantly altered across all three exposure groups (online suppl. Table S4). When comparing the top 10 upregulated and downregulated GO biological pathways within each male exposure group, there were 5 and 4 consistently identified pathways, respectively (Fig. 4). The 5 upregulated were mRNA processing, RNA splicing via transesterification reactions, mRNA via spliceosome, mitochondrial ATP synthesis coupled electron transport, and respiratory electron transport chain. The 4 downregulated were nervous system development, positive regulation of synaptic transmission, regulation of cationic channel activity, and axonogenesis. When looking at unique GO biological processes within exposure groups, OPFR exposure affected more GO processes related to mitochondrial function (13 upregulated) than BFR exposure (5 upregulated). BFR exposure uniquely produced effects on GO processes related to RNA replication including RNA polymerase 1 promotor (3 upregulated) and helicase activity (2 upregulated). The GO biological processes specific to OPFR males were related to axon guidance/neuron growth (6 downregulated) and synaptic signaling (3 downregulated) (online suppl. Table S5).

In the OPFR females, 16 GO biological pathways were significantly upregulated. Of the 16, 11 were also identified in the OPFR-exposed males, 7 of which were related to mitochondrial regulation and in the top ten significantly altered pathways in males (Fig. 4). No GO biological pathways were significantly downregulated in the OPFR females.

In males, 24 FM 550, 31 BFR, and 42 OPFR potential hub PPIs were identified for the upregulated DEGs, with 18 commonly predicted across all three groups. These were (in no specific order) ESR1, HTT, POU5F1, ILF3, ILF2, POLR2A, TARDBP, BRCA1, MYC, NANOG, RAD21, CTCF, TBP, TRIM28, TP53, UPF1, HSF1, PML. In the male downregulated DEG sets, only 1 FM 550, 3 BFR and no OPFR PPIs were identified. Downregulation of RXRA was predicted in both FM 550 and BFR-exposed males (Table 1). No statistically significant PPIs were predicted for either the up- or downregulated DEGs in the OPFR females.

In males, 71% of upregulated KEGG pathways were shared across FM 550, BFR, and OPFR-exposed animals. Many of the 22 pathways common to all three were related to neuronal degenerative diseases including Huntington disease, Parkinson’s disease, Alzheimer’s disease, and oxidative phosphorylation (Fig. 5A). Pathways unique to component class included ribosome biogenesis in eukaryotes, circadian rhythm for OPFR males, and mitophagy for BFR males. There were also four KEGG pathways only observed to be upregulated in FM 550 males: nucleotide excision repair, pyrimidine metabolism, salmonella infection, and shigellosis. By contrast, 2 FM 550, 4 BFR, and 10 OPFR KEGG pathways were downregulated in males (online suppl. Table S6). The only pathway shared between all male exposure groups was endocytosis. In the OPFR males, 8 KEGG pathways were exclusively downregulated including cholinergic and glutamatergic synapse and long-term depression.

In the OPFR females, there were 13 upregulated and 6 downregulated KEGG pathways. Twelve of the 13 upregulated KEGG pathways were also upregulated in OPFR-exposed males (Fig. 5B). These were oxidative phosphorylation, Parkinson disease, thermogenesis, Prion disease, Huntington disease, diabetic cardiomyopathy, amyotrophic lateral sclerosis, nonalcoholic fatty liver disease, pathways of neurodegeneration, Alzheimer disease, retrograde endocannabinoid signaling, and cardiac muscle contraction. Two downregulated KEGG pathways were also common to both sexes: endocytosis and axon guidance.

Five upregulated KEGG pathways related to neurodegeneration were shared between all male exposures and OPFR females: Parkinson disease, Huntington disease, oxidative phosphorylation, Alzheimer’s disease, and pathways of neurodegeneration. Within those KEGG pathways, 18 DEGs were commonly to all (Fig. 5C). The majority of those shared genes are associated with cellular respiration including subunits of mitochondrial ATP synthase (ATP5), the terminal enzyme of the mitochondrial respiratory chain, cytochrome c oxidase enzyme (COX), NADH:ubiquinone oxidoreductase subunits (NDUFA) important for oxidative phosphorylation system in mitochondria, and components of the Ubiquinol-Cytochrome C Reductase (UCQR) complex, which is part of the part of the mitochondrial electron transport chain.

NanoString was used to validate aspects of the transcriptional data including 19 upregulated genes, 10 downregulated genes, and 16 predicted PPIs (online suppl. Table S1). The average gene expression count ranges were 5–92 for negative controls, 138–89,160 for positive controls, and 836–4,148 for housekeeping genes. The NanoString results were largely concordant with the RNAseq results and pathway analyses with some exceptions. Three genes of interest, POU5F1, OXTR, and FOXP3, were found to have gene expression counts below those of the negative controls and thus considered not detected. This is not surprising since none are highly expressed in PND 1 cortex. Similarly, 4 upregulated DEGs (GNG5, Ndufa2, CLOCK, HTT), 1 downregulated DEG (CACNG) and 3 PPIs (ESR1, BRCA1, and RXRA) had expression levels below the lowest expressing housekeeping gene, Gusb, and were thus considered low expressors or undetected. These genes are also not highly expressed in the developing cortex. Thus, confirmation of these PPIs would require examination of regions that more strongly express them. Full results are presented in online supplementary Table S7.

A total of 443 lipids (208 lipids in negative ion mode and 235 in positive) were identified from a target list of 778 lipids with LC-IMS-CID-MS. As illustrated in Figure 6a, the BFR-exposed pups had 29 significantly disrupted lipids, while those exposed to OPFRs had 38. As expected, the FM 550 exposed pups had the most lipidomic changes with 76 statistically significant lipids. Of those, 16 were shared by all exposure groups, and the OPFR-exposed group had more lipid dysregulation in common with the FM 550 exposed group. Further breakdown of the statistically significant lipids is shown in online supplementary Table S8. The disrupted lipids were then clustered by structural similarity and saturation or unsaturation of the fatty acyl chains (Fig. 6b). Saturated free fatty acids and saturated triacylglycerols (TGs) were observed to be downregulated for all the FM 550-exposed pups and the OPFR-exposed pups, albeit to a lesser extent. Additionally, several sphingolipids, a lipid class including ceramides (Cers) and sphingomyelins (SMs), were also affected by exposure. Interestingly, Cers were upregulated by all three exposures, whereas SMs were predominantly downregulated in only the FM 550 exposed group.

Analysis by sex revealed strong sex differences in cortical lipid composition and impact of exposure. In the unexposed controls, females had 101 statistically significant lipids with lower abundances than the males and only 2 with higher abundances. A high degree of dichotomy was also observed between the exposure groups within sex (Fig. 7; online suppl. Table S9). Specifically, while the BFRs did not show strong lipid class-based trends for either sex, OPFR exposure caused SMs to decrease more in males than females. Additionally, while the male pups did not show any Cer dysregulation, many Cer species increased in the females in all three exposure groups. Both female and male rat pups also showed considerable Cer upregulation when exposed to FM 550, but only the males had substantial dysregulation of SMs. Both males and females in the FM 550 group had upregulation of unsaturated TGs and downregulation of saturated TGs, illustrating an exposure-based effect on double bond presence and absence for the fatty acyl groups in the lipids. When comparing commonalities between sexes for the BFR and OPFR-exposed pups, there was very little overlap among the dysregulated lipids. Female and male pups therefore showed more dysregulation in common with respect to the whole FM 550 mixture as opposed to its individual components.

Both the transcriptomics and lipidomics analyses indicated that the OPFRs were the most biologically active component of the FM 550 mixture. Strong baseline sex differences were observed and, consequently, also highly sexually dimorphic responses to prenatal FR exposure. Both component classes and the full FM 550 mixture appear capable of impairing aspects of mitochondrial function including oxidative phosphorylation, electron transport, and respiratory complex assembly and function. Mitochondrial dysfunction likely underlies the exposure-related enrichment of other pathways and systems, especially the numerous neurodegenerative pathways identified. Particularly in males, genes essential to mitochondrial function (ex. ATP5, COX, NDUF, UCQR) involving all but one of the electron chain complexes (Fig. 8a) primarily drove neurodegenerative pathway enrichment. In males, OPFRs also downregulated KEGG pathways related to glutamatergic and cholinergic synapse function and the synaptic vesicle cycle. Similarly, axon guidance and endocytosis were downregulated by OPFRs in both sexes, demonstrating that developmental neurotoxicity is multimodal. The lipid analysis was largely concordant with the transcriptomics and indicative of sex-specific disruption of neurodevelopmental processes, especially myelination. Mitochondria are essential to support the energy-demanding nature of myelination, preserve membrane potential, and fuel the enormous amount of neural circuit reorganization and synaptic refinement that occurs postnatally [85, 86]. Accordingly, mitochondrial dysfunction has been linked to many different pathological conditions of the developing central nervous system more common in boys including schizophrenia and autism spectrum disorders (ASD) [87, 88]. Thus, the data generated herein support stated concerns that FR, particularly OPFR, exposure contributes to a greater risk of male-biased neurodevelopmental disorders [10, 89, 90].

Profound sex differences in cortical lipid composition and transcriptome were identified at birth and, unsurprisingly, strongly sex-specific effects of exposure were also observed. The fact that the perinatal period is a time of intense brain sexual dimorphism is well established [84, 91] and highly dependent on steroid hormones [92, 93]. In rodents, testicular androgens, generated by a perinatal testicular hormone surge, are largely aromatized in the developing male brain and masculinize it via estrogen receptors (ERs), particularly the alpha form of ER (ERα) [94]. In humans, although a similar androgen surge occurs in boys, brain masculinization requires androgen receptors, with ERs playing a vastly less significant role [95, 96]. Thus, at birth, both rodent and human males have higher levels of circulating steroid hormones than females that drive brain sexual differentiation.

Accordingly, many transcripts and lipids with sexually dimorphic expression at birth are steroid sensitive, especially to estrogens [97]. Significantly, at least one prior study using a combination of wildtype and ERα knockout mice found that OPFR exposure disrupts aspects of metabolic homeostasis with males displaying reduced energy intake and body weight but not ovariectomized females [98]. Also reported, was that perinatal exposure to the same OPFR mixture (tris[1,3-dichloro-2-propyl]phosphate [TDCPP], TPHP, and tricresyl phosphate [TCP], each at 1 mg/kg) via oral exposure to the mouse dam (GD 7 – PND 14) resulted in upregulation of hypothalamic ERα expression on PND 14 in females but not males [99]. This upregulation was not observed at birth. Hub PPI performed herein identified ESR1 as a possible upstream target for both the BFRs and the OPFRs, as well as the full mixture, an observation consistent with the published mouse data. This prediction is also concordant with evidence reported by our group and others using a variety of in vitro and in vivo models demonstrating that FM 550 and its components are estrogen disrupting and can bind ERs including ERα [28, 100, 101]. However, our attempt via NanoString analysis to test if ERα expression was disrupted was unsuccessful because it is not highly expressed in the PND 1 rat cortex, thus it would have to be assessed elsewhere in the newborn brain.

Both brain RNA expression patterns and lipid composition, including mitochondrial lipid composition, are strongly sexually dimorphic at birth, with region-specific differences [102, 103]. For example, cortical cardiolipin content is not dimorphic across early neonatal development but saturation is, with females having a higher saturation ratio at birth then dropping to male typical levels by PND 4 [104]. This is strongly concordant with neonatal patterns of brain ER expression, which can be highly sexually dimorphic at birth but rapidly change [105, 106]. In the cortex, the cardiolipin biosynthetic pathway is sensitive to estrogens [107], as are many other nervous system lipids including endocannabinoids [108] and their metabolism [97]. Critically, mitochondria are also highly sexually dimorphic with fatty acid utilization and functional capacities, including electron transport capacity, typically higher in adult female brain mitochondria than in male brain mitochondria [109]. Estrogens also profoundly alter mitochondrial biogenics via mechanisms involving both canonical ERs and GPER1 depending on age and tissue type [110].

Collectively, these brain structural and functional sexual dimorphisms likely underlie sex and region-specific vulnerabilities to neuropsychiatric disorders including ASD [111‒113], and adverse outcomes to environmental insults including chemical exposures [114, 115]. Here we found that the OPFRs in particular disrupted mitochondrial pathways in both sexes related to the respiratory electron transport chain, mitochondrial ATP synthesis, and NADH dehydrogenase complex assembly with the male groups also showing strong evidence of downregulated pathways associated with nervous system development, axonogenesis, synaptic transmission, and both cholinergic and gluatamatergic synapse function, suggesting that male neurodevelopment is more sensitive to exposure and, potentially, the consequences of mitochondrial disruption. This is consistent with evidence that prenatal chemical exposures increase ASD risk by perturbing mitochondrial physiology [116].

The BFRs and OPFRs both appear to have profound, albeit somewhat different, effects on mitochondrial function and, consequently, neurodevelopment, with the BFRs having a smaller sequala of effects than the OPFRs. Numerous mitochondrial GO processes were affected by both OPFR and FM 550 exposure (online suppl. Table S10), including mitochondrial respiratory chain complex assembly (especially Complex I), indicating that the OPFRs are likely the primary drivers of FM 550 mitochondrial stress. Additionally, OPFR exposure exclusively affected 2 mitochondrial GO process: establishment of protein localization to mitochondrial membrane (GO:0090151) and mitochondrial RNA metabolic process (GO:0000959). Our OPFR results are consistent with prior work in the nematode (Caenorhabditis elegans) led by the National Toxicology Program (NTP), which reported evidence of significant OPE-related mitochondrial toxicity, particularly for TPHP, BDHP, and TMPP, with various degrees of cytotoxicity at higher doses [32]. The authors concluded that OPFRs are comparable to PBDEs on mitochondrial toxicity. This is significant because PBDE exposure has been linked to heightened ASD risk [89], purportedly via mitochondrial toxicity [117].

By contrast, the BFRs primarily only affected GO processes related to mitochondrial ATP synthesis and electron transport, suggesting they may have a somewhat different mode of action on mitochondria than the OPFRs. BFRs were also linked to upregulated responses to hypoxia and catabolic processes including ubiquitin-dependent protein catabolic processes. At least to some extent, thyroid hormones regulate the maturation of brain mitochondria and the electron transport chain [118], and prenatal hypothyroidism has long been known to induce severe neurodevelopmental abnormalities [119]. Many legacy BFRs, most notoriously the phased out PBDEs, impair thyroid signaling. Prior studies on the newer BFRs, however, have produced inconsistent results and collectively suggest that the thyroid-disrupting potential of the BFRs tested herein is minimal at best, particularly compared to the PBDEs [43, 120, 121].

Upregulated KEGG pathways in both sexes and all exposure groups included numerous pathways related to neurodegeneration including Parkinson’s, Alzheimer’s, and Huntington diseases. While this may seem counterintuitive for a newborn brain, the DEGs shared across these pathways are strongly involved in the electron transport chain and oxidative phosphorylation (Fig. 8). Thus, critically, the enrichment of neurodegenerative pathways in the FR exposure groups likely does not indicate that developmental exposure contributes to a higher risk of later-in-life disorders such as Alzheimer’s Disease per se, but rather emphasizes the fundamental significance of mitochondrial health for nervous system health and cognitive development (for example, see Fig. 8b). Within the Alzheimer’s disease KEGG pathway, for example, exposure upregulated genes integral to four of the five electron transport chain complexes plus TUBB, which is associated with axonal transport defects, demonstrating the importance of these genes for not only neurodegeneration but also neurodevelopment. Finally, all FR groups produced evidence that fundamental processes such as those related to protein processing in the endoplasmic reticulum and RNA transport were upregulated, which also signifies neurodevelopmental disruption.

In the OPFR-exposed animals, aspects of endocytosis were disrupted as well as pathways for long-term depression (LTD) and vesicle, vacuolar and cytoplasmic vesicle membrane formation, suggesting disruption of neurotransmitter systems and function. Accordingly, downregulation of genes required for axon guidance was found in both sexes (Fig. 9a, b), as well as genes related to glutamatergic and cholinergic synapses specifically in males (Fig. 9c, d), including ACHE and CAM2KA. Downregulation of the axon guidance KEGG pathway by OPFRs (Fig. 9b) is consistent with our prior observation that prenatal exposure to OPFRs alters serotonergic innervation of the rat fetal forebrain [29]. Critically, glutamate, the most prevalent neurotransmitter in the frontal cortex, is packaged and released from presynaptic vesicles. Glutamatergic neurons are capable of many different forms of LTD and thus many forms of neuroplasticity [122] one of which, endocannabinoid-dependent LTD, is common in the forebrain to both glutamatergic and GABAergic synapses, particularly in the prefrontal cortex [123]. This form of synaptic plasticity can, among many things, modulate dopamine signaling and affect motivated behaviors including addiction. Forebrain cholinergic systems regulate behaviors including attention, habit formation, memory, motor activity, and reward, and arise from neuronal precursors that could either become GABAergic or cholinergic, with most cholinergic cortical innervation occurring just after birth in the rat. Acetylcholine is critical for cortical refinement including dendrite stabilization, neuritogenesis, and synaptogenesis [124]. Mitochondria have an important role in supporting synaptic activity and the synthesis and packaging of some key neurotransmitters, including glutamate and choline [125, 126]. Thus, generally, OPFRs appear capable via multiple mechanisms of altering vesicle-related synaptic activity including amino acid neurotransmitter packaging and release, and cortical organization and plasticity. Significantly, our results support prior work illustrating that OPEs, including TPHP, while perhaps not as capable of impacting acetylcholinesterase activity as their predecessors, impact other aspects of choline metabolism and signaling [127‒129].

A notable endpoint where the full mixture was more impactful than either component class was the exposure- and sex-dependent impact of Cer levels. While both the OPFRs and the BFRs upregulated certain Cer species in females, Cer upregulation in males was only observed in the FM 550 group. This suggests that males may be more resilient to the exposure of the individual FM 550 components, and that the full mixture is needed to overcome this resilience and produce disruption. Similarly, SMs, which are derived from Cers, essential for myelin formation, and ubiquitous in cell membranes, were profoundly downregulated in FM 550 males with females seemingly more resistant.

Cers and SMs are intensely studied in the context of neurodegenerative disease because they are associated with Alzheimer’s disease and Multiple Sclerosis, but their roles in development are less well delineated. Interestingly, a prior rat study identified a transient dip in brain Cer levels at birth but did not account for sex [130]. Similarly, sphingosine levels were shown to be markedly higher than dihydrosphingosine levels. Further work using a wider range of OPFR and BFR doses will be needed to more comprehensively explore the sex difference in Cer and SM sensitivity. However, the observation that dysregulation was so sex specific is compelling, especially since Cers are beneficial for the early growth and development of neuronal cells, and sphingolipids are critical components of cell membranes and potent signaling molecules in many pathophysiological processes including neuroinflammatory processes [131]. Accumulation of brain Cers has also been implicated in energy balance disorders and, because de novo synthesis occurs in the smooth endoplasmic reticulum and the mitochondria, their disruption lends further support to the conclusion that FM 550, and the OPFRs in particular, disrupt mitochondrial function [132]. Cer concentrations are precisely controlled with cellular growth or death resulting from very slight shifts [133]. Cers are highly sensitive to endogenous estrogens and essential to key aspects of neural and synapse development including myelination [134]. More work is needed to understand how brain Cer composition may be vulnerable to FR exposure and what the resulting consequences may be.

Another finding of interest from the lipidomic studies is the upregulation of TGs with unsaturated fatty acids and downregulation of TGs with saturated fatty acids for both the males and females exposed to FM 550. Saturated and unsaturated fatty acids are known to have different health effects and this split dysregulation indicates disruption of specific enzymatic activity occurring for the different TGs upon exposure to the FM 550 mixture. To our knowledge, no prior study has reported this split, likely because many fatty acid analyses are performed with gas chromatography coupled to mass spectrometry, which cleaves the fatty acyl groups from the lipid headgroups. This causes a pool of all fatty acids in the sample and ultimately the inability to distinguish specific lipid species trends. Therefore, we find this result quite novel and a unique benefit of our approach that we will follow up in future studies.

A few of the GO processes associated with FM 550 or BFR exposure involved modulation of viral and related processes. While this may at first appear odd, genes in this process include IGF2R and APOE, which have significant roles in embryonic development and neural function. APOE was downregulated in both the BFR and FM 550 males. APOE codes for apolipoprotein E (ApoE), which is expressed at the surface of TG-rich lipoproteins, chylomicrons, and very low density lipoproteins (VLDL), acting as a receptor binding ligand and is critical for lipid transport in the brain [135, 136]. In humans there are multiple APOE alleles, with apoE4 conferring the strongest genetic risk factor for Alzheimer’s disease [137].

The hub PPI analysis predicted downregulation of RXRA for both FM 550 and the BFRs but not the OPFRs. RXRs are common binding partners to many other nuclear receptors, especially PPARs, liver X receptors (LXRs), and vitamin D receptors (VDRs). The PPI result is concordant with prior work in rat placenta identifying disruption of LXR/RXR pathways as a primary outcome of FM 550 exposure [28]. Notably, the RXRA/PPARα heterodimer is required for PPARα transcriptional activity on fatty acid oxidation genes such as ACOX1 and the P450 system genes, and plays a significant role in oxidative stress, energy homeostasis, mitochondrial fatty acids metabolism, and inflammation in the brain [138]. Genes in the PPAR/RXR family are also involved in aspects of neural proliferation and differentiation, myelination, and synaptic plasticity and have been associated with neurodevelopmental disorders including schizophrenia [139]. Prior animal and in vitro studies suggest that FM 550 and the OPFRs, particularly TPHP and ITP, may have adipogenic effects via PPARγ activation and PPARγ-mediated aP2 enhancer activity [45, 100, 140]. The combined effects of OPFR exposure on PPARγ-related targets and BFR exposure on LXR/LXR endpoints could possibly account for some unique and/or absence of effects by the FM 550 mixture.

Collectively, the data revealed multiple adverse modes of action for the BFRs and OPFRs on neurodevelopment. Most significantly, both induced mitochondrial disruption, which is associated with numerous neurodevelopmental disorders, including ASD, and neurodegenerative disorders. The OPFRs were most disruptive, and via multiple mechanisms including disruption of cholinergic and glutamatergic systems. The discovery that the lipid composition of the PND 1 cortex is so highly sexually dimorphic provides additional and novel insight as to how and why the male and female brain responds differently to environmental and other exogenous exposures and stresses. That pups of both sexes showed more lipid dysregulation in common with respect to the FM 550 mixture than its individual components indicate that exposure to the different FR classes affects the sexes differently on a molecular level. Lipid classes with the most notable changes included the Cers, SMs, and TGs. Since dysregulated Cers are strongly associated with Alzheimer’s disease risk, the robust upregulation of Cers observed in the OPFR-exposed females could suggest heightened risk of brain metabolic disease at a very early age. Downregulation of phosphatidylcholines has also been noted in Alzheimer’s disease due to choline shuttling, and downregulation was observed herein in pups exposed to OPFRs and the FM 550 mixture. Finally, this work emphasizes how essential it is that future studies regarding the impact of FRs and other exposures on neurodevelopment focus on sex differences. Too many resources and published studies remain either highly male-biased or sex-agnostic [141].

The authors are grateful to Heather Stapleton, with whom we have collaborated for years, for preparing the dosing mixtures, as well as David (Andy) Baltzegar, Director of the Genome Science Laboratory at NC State, for transcriptomics design consultation, library preparation, and RNA sequencing. We also thank the staff in the Biological Resources Facility for their assistance with animal care and husbandry. All lipidomic measurements were performed in the Molecular Education, Technology, and Research Innovation Center (METRIC) at NC State University.

All animal work was reviewed and approved by the NC State Institute of Animal Care and Use Committee (IACUC; approval number 20-150) and supervised by the attending veterinarian. All animal work conformed to all applicable federal and state US laws.

The authors have no known or perceived conflicts of interest.

This work was supported by RO1ES028110 to H.B.P., U01ES026717 to D.L.A., and P30ES0251278 to Rob Smart. M.G.D., M.T.O., and E.S.B. were supported by P42 ES027704 and P42 ES031009 and a cooperative agreement with the United States Environmental Protection Agency (STAR RD 84003201).

Shannah K. Witchey led the design and execution of the studies; data collection; data analysis; supervision of assisting personnel; manuscript preparation including all figures, tables, and online supplementary materials; and gave final approval of the version to be published. Michael G. Doyle performed data analysis for the lipidomics portion of the project and contributed to manuscript preparation and final approval of the version to be published. Jacob D. Fredenburg performed the transcriptomics analysis in collaboration with Shannah K. Witchey and helped prepare the related figures and tables. Genevieve St. Armour assisted with animal dosing and care, necropsy, brain microdissection, and the nanostring-related work including the preparation of the related methods and results sections. Brian Horman organized and oversaw all of the animal dosing, housing, and necropsy as well as IACUC protocols, and also assisted with manuscript editing including approval of the final version. Melanie T. Odenkirk extracted the lipids, aided in analyzing the lipidomics portion of the project, and assisted with the manuscript including final version approval. David L. Aylor led, supervised, and partially funded the transcriptomics work including data analysis and interpretation, drafting of the transcriptomics portion of the manuscript, and revising it for intellectual content. Erin S. Baker led, supervised, and partially funded the lipidomics work including the development of appropriate methods, drafting of the lipidomics portion of the manuscript, and revising it for intellectual content. Heather B. Patisaul conceived of the project, funded the majority of the work, oversaw the entire study including the work in the coordinating labs, contributed to the drafting of the work, and gave final approval of the version to be published.

All data generated or analyzed during this study are included in this article, its online supplementary materials, and the GEO (GSE211430) and MassIVE (MSV000090203) repositories. Further inquiries can be directed to the corresponding author.