Calcium signaling as a possible mechanism behind increased locomotor response in zebrafish larvae exposed to a human relevant persistent organic pollutant mixture or PFOS

Persistent organic pollutants (POPs) are widespread in the environment and their bioaccumulation can lead to adverse health effects in many organisms. Previously, using zebrafish as a model vertebrate, we found larvae exposed to a mixture of 29 POPs based on average blood levels from the Scandinavian population showed hyperactivity, and identified perfluorooctanesulfonic acid (PFOS) as the driving agent for the behavioral changes. In order to identify possible mechanisms, we exposed zebrafish larvae from 6 to 96 h post fertilization to the same mixture of POPs in two concentrations or a single PFOS exposure (0.55 and 3.83 μM) and performed behavioral tests and transcriptomics analysis. Behavioral alterations of exposed zebrafish larvae included hyperactivity and confirmed previously reported results. Transcriptomics analysis showed upregulation of transcripts related to muscle contraction that is highly regulated by the availability of calcium in the sarcoplasmic reticulum. Ingenuity pathway analysis showed that one of the affected pathways in larvae exposed to the POP mixture and PFOS was calcium signaling via the activation of the ryanodine receptors (RyR). Functional analyses with RyR inhibitors and behavioral outcomes substantiate these findings. Additional pathways affected were related to lipid metabolism in larvae exposed to the lower concentration of PFOS. By using omics technology, we observed that the altered behavioral pattern in exposed zebrafish larvae may be controlled directly by mechanisms affecting muscle function rather than via mechanisms connected to neurotoxicity.


Introduction
Persistent organic pollutants (POPs) are organic compounds that are resistant to environmental degradation. POPs include polychlorinated biphenyls (PCBs), dichlorodiphenylsichloroethanes (DDTs), brominated flame retardants (BFRs), dioxins and per-and poly-fluoroalkylated substances (PFASs) (UNEP, 2005) with perfluorooctanesulfonic acid (PFOS) being the one most frequently detected in the environment (Paul et al., 2009). Many of the POPs were used, or are still used today, as pesticides, industrial chemicals, solvents, and pharmaceuticals, and most of them are man-made. Their persistence in the environment and high lipid solubility means that POPs can readily bioaccumulate in fatty tissues (Ritter et al., 1998). As such, POPs are detected at high concentrations in wildlife, such as birds, fish, and marine mammals, but also humans (Giesy and Kannan, 2001;Boon et al., 2002;Chen and Hale, 2010;Porta et al., 2012;Kato et al., 2015;Olsen, 2015;Jepson et al., 2016). In humans, POP exposures mainly occur through food consumption, but also via drinking water, outdoor and indoor air, and from the working environment (EFSA, 2008;Guo et al., 2019;WHO, 2020).
Exposures to POPs have been found to lead to diverse effects on health. For example, POP exposures have been linked to the rise of type 2 diabetes and obesity (Ruzzin, 2012;Thayer et al., 2012;Guo et al., 2019), endocrine disruption (Li et al., 2008;Gregoraszczuk and Ptak, 2013), hormone-related cancers (Brody et al., 2007), the disruption of sexual development (Vallack et al., 1998;Sanderson, 2006), and cardiovascular disease (Lind et al., 2012;Ljunggren et al., 2014). Of more recent concern are early life exposure to POPs in utero or during breast feeding, since exposure during this sensitive period of development can have persistent effects on different life traits such as neurodevelopment. Early life exposures to POPs for instance have been associated with effects on cognitive functions such as poorer performance in intelligence tests, language, processing speed, short term memory, and association with traits relating to ADHD including problems in impulse https://doi.org/10.1016/j.envres.2020.109702 Received 19 February 2020; Received in revised form 30 April 2020; Accepted 18 May 2020 control, hyperactivity, and attention (Hoffman et al., 2010;Lam et al., 2017;Pessah et al., 2019).
Due to their persistent nature, POPs are rarely found as single elements within individuals, but rather occur in mixtures and can be passed on from mothers to developing embryos (Porta et al., 2012;Pumarega et al., 2016;Berntsen et al., 2017;Haug et al., 2018). This information makes the study of mixture effects of high importance since it is hard to predict toxic effects based on single chemical exposure as some of these chemicals may interact in a synergistic, additive, or antagonistic manner (Bopp et al., 2018).
Risk associations of exposure to POPs and related adverse health effects have not only been obtained from epidemiological studies, but also from in vitro and in vivo testing (Brody et al., 2007;Ruzzin, 2012;Guo et al., 2019;Pessah et al., 2019). Animal testing in particular has provided a plethora of evidence of the effects of chemicals on the physiological and molecular level (Parasuraman, 2011). Zebrafish (Danio rerio) are a freshwater fish that has been used extensively as a model organism in human and environmental toxicology due to its many advantages such as rapid development, small size, embryonic transparency, short generation time, and fully sequenced genome (Hill et al., 2005). POP exposures in zebrafish show similarity to effects in mammals (Hallgren et al., 2001;Keil et al., 2008;DeWitt et al., 2009;Chan and Chan, 2012;Martínez et al., 2019;Parsons et al., 2019), including effects on behavior (Johansson et al., 2008;Huang et al., 2010;Jantzen et al., 2016). As the fundamental processes of neurodevelopment are conserved among species (Lein et al., 2005), this makes zebrafish a powerful translational model for human health (Kalueff et al., 2016). Previous work on POPs has shown that larval zebrafish exposed to PCBs displayed a decreased avoidance behavior when presented with a visual stimulus (Lovato et al., 2016). In other studies, exposure to PFASs led to either hyperactive or hypoactive locomotor activity, burst activity and startle response depending on the concentration and age of testing in zebrafish larvae (Huang et al., 2010;Spulber et al., 2014;Jantzen et al., 2016;Khezri et al., 2017;Menger et al., 2020). Exposure to BDE-47 led to an increase in spontaneous movement, a decrease in touch response, and a decrease in swimming speed, in a dose-dependent manner (Chen et al., 2012). In a large-scale study by the National Toxicology Program, behavioral screening of BFRs revealed mainly hypoactivity of zebrafish larvae whereas exposure to DDT and other organochlorine pesticides causes abnormal behavior i.e. hyperactivity or constant movement (Dach et al., 2018).
Previously we observed behavioral effects on zebrafish larvae that were exposed from 6 to 96 h post fertilization (hpf) to a mixture of POPs . This POP mixture is based on the average levels of chemicals found in human blood of the Scandinavian population and consists of 29 chemicals including PCBs, BFRs, organochlorine pesticides and PFASs (Berntsen et al., 2017). The aim of this study was to employ zebrafish as a model organism to confirm the behavioral effect observed during developmental exposure to POPs previously observed in our group and employ another behavioral test, the thigmotaxis assay, to evaluate the effect of the POP mixture on anxiety levels of zebrafish larvae. We also included a single PFOS exposure since this chemical is suspected to be solely responsible for the behavioral effect attributed to the mixture exposure . Furthermore, we performed RNA-seq analysis on exposed larvae with the aim of elucidating possible modes of action for PFOS and the POP mixture. Finally, chemical analysis was performed on exposed larvae for the evaluation of the uptake and accumulation of each chemical within the POP mixture.

Fish maintenance and breeding
The experiments were performed at the Norwegian University of Life Sciences (NMBU), Oslo, Norway, which is licensed by the Norwegian Animal Research Authority (NARA) (www.mattilsynet.no).
The experiments were carried out under the regulations approved by the unit's animal ethics committee (Institutional Animal Care and Use Committee/IACUC) following Norwegian laws and regulations controlling experiments and procedures on live animals in Norway.
AB wild-type (AB) zebrafish were kept at 28 ± 1°C under a 14:10 light/dark photoperiod. Animal care was performed in accordance with local protocols and more information can be found in Supplementary data. For embryo production, adults were placed in breeding tanks in the afternoon with a divider separating males and females. The next morning the separator was removed as soon as the lights turned on (08:00) and embryos were collected 1 h later. Embryos were maintained in sterile embryo media (60 μg/mL Instant Ocean® sea salts) until the time of exposure.
Stock solutions of PFOS, dantrolene and the POP mixture were prepared in DMSO. A fresh stock solution of caffeine was made the day of the testing (diluted in embryo media) whereas stock solutions of PFOS, dantrolene, and the POP mixture were stored at −20°C.

Solutions preparation
For all experiments two concentrations of the POP mixture were used. The low concentration was equal to the levels of chemicals that are 10 times higher than what is found in average Scandinavian human blood levels and the high concentration corresponds to levels 70 times higher than the average (exposures will be referred to as POP10 and POP70 in the remainder of the manuscript). The solutions were prepared on the day of the experiments by diluting the stock solution (1,000,000x) in sterile embryo media and adjusting the DMSO concentration to 0.1%. The concentrations of the POP mixtures were chosen based on previous studies performed in our group where the 70x concentration was the lowest concentration with an observable behavioral effect whereas the 10x concentration was the highest concentration with no observed behavioral effect . Concentrations of PFOS were based on the nominal concentration of PFOS in the POP mixture. The low concentration of PFOS corresponded to 0.3 mg/L (0.55 μM, will be referred from now on as PFOS10) and the high concentration was equal to 2.06 mg/L (3.83 μM, will be referred to as PFOS70). PFOS solutions were prepared the day of the experiments by diluting the stock solution of PFOS (54.8 mM) in sterile embryo media and adjusting the DMSO concentration to 0.1%. The dantrolene concentration was chosen based on the current literature (Brennan et al., 2005;Yuen et al., 2013). Working solutions of dantrolene (50 μM) were prepared the day of the experiment by diluting the stock solution (7.4 mM) in sterile embryo media. Because of low solubility of dantrolene in DMSO, the final concentration of DMSO in the working solution was 0.68%. All exposures (Control, POP10, POP70, PFOS10, PFOS70) therefore for the dantrolene experiments were adjusted to a DMSO concentration of 0.68%. Caffeine concentration was selected initially at 250 μM based on Schnörr et al. (2012) but produced an overt toxic effect on larval zebrafish following chronic exposures (see Supplementary data Table S2), hence the concentration was adjusted to 50 μM. Working solutions of caffeine were prepared on the day of experiment in sterile embryo media and the concentration of DMSO was adjusted to 0.1%. Before and after all the experiments, the embryos were reared in an incubator at 28 ± 1°C. The light cycle in the incubator was set at 14:10 light/dark (lights on 07:30/lights off 21:30).

Exposures for behavioral tests
For the light-dark transition test, fertilized embryos were transferred into clear polystyrene 96 well plates (Nunc™ MicroWell™) with one embryo per well and exposed under static conditions from 6 hpf until the time of testing at 96-100 hpf (between 9:00-13:00) in 200 μL of media. For the thigmotaxis assay embryos were placed in 24 well plates (Corning® Primaria™) with one embryo per well in 1 mL exposure media. Embryos were exposed either to the POP mixtures in two concentrations, POP10 and POP70, or two concentrations of PFOS, PFOS10 and PFOS70. Each well plate also included a control with embryos immersed in sterile embryo media (DMSO 0.1%). For the thigmotaxis assay two controls were used, one for the POP mixture treatment and one for the PFOS treatment respectively, as the plate layout meant that each group could not be equally represented on each row and column without the addition of the extra controls. All well plates were preincubated with the respective exposure for 24 h. The media was removed before the embryos were placed in the well and replaced with a fresh working solution. All groups were spread equally on each row and column to avoid position-bias during behavioral testing. For the lightdark transition test, each well plate included 10 embryos per condition and was repeated in five independent experiments. For the thigmotaxis assay, each well plate contained three embryos per condition with three well plates per replicate and was repeated in three independent experiments.

Exposure for transcriptome analysis
Fertilized embryos were exposed to five different conditions (Control, POP10, POP70, PFOS10, PFOS70) in six well plates (Falcon® 6-well). Each well plate included only one condition with four pseudoreplicates. In each well 15 embryos were placed with 3 mL exposure media. All wells were pre-saturated for 24 h prior to the experiment and the media was replaced with fresh working solution at the start of the experiment. Embryos were exposed in static conditions from 6 hpf until 96 hpf. Replicates of 10 non-deformed embryos per experimental condition were collected, snap-frozen in liquid nitrogen and stored at −80°C until further analysis (RNA extraction for high-throughput sequencing analysis and RT-qPCR for gene validation purposes).

Exposures for behavioral tests with RyR agonist and antagonist
Based on the transcriptome analysis, we designed an experiment to uncover a potential mechanism of action of POPs on the ryanodine receptor (RyR) using the light-dark transition behavioral test. Caffeine, an established agonist, and dantrolene, an established antagonist of RyR were used (Brennan et al., 2005;Hernández-Fonseca and Massieu, 2005;Yuen et al., 2013). Each compound was used in two different exposure scenarios. In the chronic exposure the larvae were exposed from 6 hpf until the behavioral testing at 96-100 hpf (between 9:00-13:00). In the acute exposure the compound was added in the well 10 min prior to the behavioral testing. Each compound and exposure scenario (chronic and acute) were tested in a separate well plate. Six conditions were included in each well plate, Control (DMSO 0.1% for the caffeine study or DMSO 0.68% for the dantrolene study), Control with DAN or CAF, POP70, POP70 with DAN or CAF, PFOS70 and PFOS70 with DAN or CAF. We chose the highest concentration of the chemicals in question because only the high concentrations had an effect on the behavioral outcome of zebrafish. For each test, 96 well plates (Nunc™ MicroWell™) were used with one embryo per well in 200 μL exposure media. Three independent experiments with 16 individuals per condition were performed. All groups were spread equally over each row and column. All wells were pre-incubated prior to the start of each experiment.

Behavioral assays
Behavioral assays were performed in a ViewPoint® Zebrabox and its tracking software (ViewPoint Life Sciences, Lyon, France). Behavioral tests were conducted between 9:00-13:00 in 96-100 hpf zebrafish larvae. For the light-dark transition test and the thigmotaxis test, the cumulative distance travelled and the time spent active were measured simultaneously for all larvae in a well-plate. The light-dark cycle lasted 20 min and consisted of 10 min of light and 10 min of dark. Prior to the test the larvae were acclimatized in the behavioral chamber for 10 min with the lights on. For the light-dark transition test the mean swimming speed was calculated by dividing the cumulated distance travelled by the total time spent active. For the thigmotaxis test the percent of the total distance moved in the outer zone was calculated. For this purpose, each well in a 24 well plate was divided in two zones (total diameter of each well 16.2 mm). The width of the outer zone was set at 5 mm relative to the border of the well. The light level was set to 100% on the ViewPoint software (7.45 Klux, TES 1337 light meter). Infrared light (850 nm) tracks larval activity during the "dark" periods. The threshold for determining movement was set at 5 mm/s. Dead and deformed (coagulated, unhatched, notochord deformations, yolk sac or cardiac edema, swim bladder development, loss of equilibrium) larvae were excluded from the analysis (Supplementary data, Table S2).

RNA purification
RNA from samples for RNA-seq analysis and gene expression for validation were purified using NucleoSpin® RNA extraction kit (Macherey-Nagel, Germany). RA1 lysismix (Macherey-Nagel, Germany) was added to each sample (10 larvae per sample) and samples were passed through a 21-gauge needle (HSW HENKE-JECT®, Germany) until complete homogenization. Total RNA was extracted from samples following manufacturer's instructions. Each sample was eluted in 50 μL RNase-free water and stored at −80°C until further analysis. RNA purity for all samples was assessed using a Nanodrop-1000 Spectrophotometer (NanoDrop Technologies, DE, USA). For samples sent for RNA-seq analysis, RNA integrity number (RIN) was determined with Agilent 2100 Bioanalyzer (Agilent Technologies, Ca, USA) using RNA Nano LabChip Kit (Agilent Technologies, Ca, USA). All samples were found to be of acceptable quality for sequencing (RIN > 9.0).

RNA-seq and transcriptome analysis
Samples for sequencing were sent to Novogene (Hong Kong). A total of 2 μg was used for library preparation. A quality check (QC) of total RNA was performed with a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), Agilent 2100 (Agilent Technologies, CA, USA) and agarose gel electrophoresis prior to library construction. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations (more information about RNA-seq analysis can be found in Supplementary data).
For the transcriptome analysis raw fastq files were adapter trimmed using trim_galore (v0.4.5, Babraham institute, UK) under standard parameters. We used the STAR aligner (v2.5.4 b) (Dobin et al., 2012) to align and map sequences to the zebrafish genome (GRCz11, https:// www.ensembl.org) with a recent release of the zebrafish transcriptome GTF (v92, www.ensembl.org). After alignment, the generated BAM files were loaded to SeqMonk sequence analysis tool (v1.41, Babraham institute, UK) and mRNAs were quantified using the built-in mRNA seq pipeline. Transcriptome analysis was performed as described in Hurem et al. (2018) with minor modifications (see Supplementary data). Differentially expressed genes (DEG) were chosen based on a false discovery rate, FDR < 0.05 and an absolute fold change, FC > 1.5.
A principal component analysis was performed in all expressed genes that had a log 2 ≥ 0 expression in all groups (Control, POP10, POP70, PFOS10, PFOS70) using ClustVis, a web tool for visualizing clustering of multivariate data (Metsalu and Vilo, 2015). PCA scores were loaded to R version 3.6.1 (R Development Core Thayer et al., 2012, http://www.r-project.org) and biplots of principal components were designed with "ggplot2" library, while the stat_ellipse argument within the "ggplot2" library was used to compute 95% confidence ellipses. Venn diagrams of DEG were created using Venny (version 2.1, https://bioinfogp.cnb.csic.es/tools/venny/).

Pathway analysis
Ingenuity Pathway Analysis (IPA) (version 49,309,495, Qiagen) was used to uncover enriched pathways in each condition (PFOS10, PFOS70, POP10, POP70). Differentially expressed genes were imported and used with the user dataset as background. Around 60% of the genes were annotated as having a human orthologue, and these genes were used for pathways analysis. IPA uses human orthologues for pathway analysis; hence we used the IPA nomenclature inside the context of pathways (e.g. MYH4) and use the official zebrafish gene nomenclature when referred to outside IPA context (e.g. myhz1.1). IPA calculates over representation of genes and gene lists involved in known pathways, upstream regulators, and diseases, using Fisher's exact tests. Furthermore, it uses the direction of the differentially expressed genes to predict activation or inhibition of pathways and upstream regulators by means of Z-scores. P values < 0.01 were considered significant. Another tool was employed to explore affected pathways, differentially expressed genes were also imported in Webgestalt (www.webgestalt. org) for KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis and gene ontology (GO) analysis.

RT-qPCR for gene validation
For gene validation analysis equal amounts of RNA were reverse transcribed with the high-capacity cDNA RT kit (Applied Biosystems, NY, USA), followed by a 10 times dilution of the cDNA reaction with nuclease free water. QPCR on the diluted cDNA were performed in 10 μL, containing 5 μL FastStart Essential DNA Green Master (Roche, Norway), 250 nM of forward and reverse primers, 2 μL of diluted cDNA, and nuclease free water, in technical duplicates. PCR was performed on a Roche Lightcycler 96 (Roche, Norway), with 10-min denaturation at 95°C, followed by 45 cycles of 10 s at 95°C and 30 s at 60°C. After the run a melting curve was generated from 65 to 95°C. Primers for validation genes were developed using the Primer-BLAST software from NCBI (https://www.ncbi.nlm.nih.gov) (Supplementary data Table S3). Primer sequences for reference genes (beta-actin, ef1a) were taken from a previous publication . All primers were validated for specificity by melting curve analysis and gel electrophoresis. Efficiency of primers was determined against a dilution curve of pooled zebrafish cDNA. Cq values were calculated using linreg (Ramakers et al., 2003). These values were used for calculating normalized gene expression using the geometric average of the 2 reference genes (betaactin, ef1a). Correlation was calculated with non-parametric Spearman correlation in Graphpad (v8.3, CA, USA).

Chemical analysis
Embryos from 6 until 96 hpf were exposed statically to either POP70 or PFOS70 in a 96 well plate (1 plate per condition) with one embryo per well in 200 μL exposure media. All wells were pre-saturated with the respective exposure for 24 h. Larvae were then collected in 2 mL eppendorf tubes. In total, 10 tubes were collected with each tube having an ascending number of embryos (the 1st tube contained 1 embryo, the 2nd tube contained 2 embryos, and so on, up to the final tube that contained 10 embryos). The reason was to make certain that chemicals would be detected during chemical analysis. The exposure media was also collected and sent for chemical analysis. Pools of 6, 7, and 8 embryos were used as replicates to evaluate the chemical uptake in the tissue of larvae zebrafish.
The measurements were performed at the Norwegian University of Life Sciences (NMBU), Department of Food Safety and Infection Biology, Laboratory of Environmental Toxicology. For the lipophilic group, extraction with cyclohexane/acetone and water was followed by gel permeation column or sulphuric acid clean-up. Separation and detection of the OCPs and PCBs were performed on a GC coupled to Electron Capture Detector (ECD) and low-resolution mass spectrometry (LRMS). Detection of BDEs and HBCD was performed on a HRGC-LRMS (Polder et al., 2014). For perfluorinated compounds, the samples were extracted with methanol and clean up was accomplished using active carbon. Further, the samples were separated by high-performance liquid chromatography (HPLC) and detection achieved by tandem mass spectrometry (MS-MS) (Bytingsvik et al., 2012). Details from the extraction, clean-up, and instrument run for the samples and quality control parameters can be found in the supplementary data.

Statistical analysis for behavioral data
Behavioral data were transferred to R version 3.6.1 (R Development Core Thayer et al., 2012, http://www.r-project.org). Dead and deformed larvae were discounted for behavioral analyses. Here, no more than 8 larvae were removed (chronically exposed to dantrolene only, Supplementary data Table S2). For all test scenarios, only motility during the dark phase was analyzed as movement was minimal during the light periods as expected based on previous work (e.g. Fraser et al. (2017)). Linear mixed effect (LME) models within the "nlme" package were used to assess behavior. For the initial light dark transition tests on exposed larvae, the dependent variable was either the cumulative time spent active (seconds), the cumulative distance travelled (mm), or average swimming speed, with group (5 levels, control, POP10, POP70, PFOS10, PFOS70) as a categorical independent variable, and replicate as a random effect. For the thigmotaxis assay, the same models were used, with the addition of a fourth dependent variable, thigmotaxis (% of time spent in the outer zone) and the independent variable group had 6 levels (control-POP, POP10, POP70, control-PFOS, PFOS10, PFOS70). For the mechanistic work, the dependent variable was either cumulative distance moved, cumulative time active, or swimming speed, with dantrolene or caffeine, (2 levels, Y/N) and treatment (3 levels, control, POPs, PFOS) as categorical independent variables, and replicate as a random effect. Here, the variables dantrolene or caffeine, and treatment were initially allowed to interact. We then compared the interaction model to a model without the interaction using the Akaike Information Criterion with a correction for small sample size (AICc). The model with the lowest AICc score was then considered the true model (Aho et al., 2014). If the interaction was the true model, this provides evidence of a mechanistic action, whereas no interaction provides evidence against a mechanistic action. For model validation, all models were also compared to a null model to verify results were not based on the random effect, also using AICc. More information about the statistical analysis parameters can be found in Supplementary data. Significance was assigned at p = < 0.05. All graphs were plotted in Graphpad (v8.3, CA, USA).

Results
3.1. Behavioral outcome of larvae exposed to the POP mixture or PFOS Statistical analysis of the light-dark transition test revealed that the cumulative distance moved and swimming speed were significantly affected by chemical exposure, but not the time spent active (Supplementary data Fig. S1). The cumulative distance moved was higher compared to controls for PFOS70 and POP70 and swimming speed increased when exposed to PFOS70, POP10, and POP70. The thigmotaxis test revealed that the percent of the total distance moved in M. Christou, et al. Environmental Research 187 (2020) 109702 the outer zone (thigmotaxis) was affected by exposure to chemicals. Larvae exposed to PFOS10, PFOS70, and POP70 moved significantly further in the outer zone of the well compared to controls (Supplementary data Fig. S2).

Sequencing analysis
Analysis with FastQC revealed high quality sequences with phred scores generally above 35 over all reads except for the paired-end read towards the 3' end where reads were still above the acceptable phred score of 20 (data not shown). Average mapping efficiency was over 90% unique reads (Supplementary data Table S4).
Quality analysis, performed in SeqMonk on aligned reads, showed a high proportion of reads that fell into genes (90%). Additionally, a high proportion of those reads fell into exons (around 94%) which means that the library contained very few unspliced transcripts. A low percent of reads was present in ribosomal and mitochondrial RNA and 87% of the annotated zebrafish genes were mapped. The reads mapped equally on the sense and anti-sense strands which confirms that our library was non-strand specific (Supplementary data, Fig. S3). A cumulative distribution analysis of expressed genes over all samples revealed the same profile over all levels of expression, which indicated highly similar sequencing libraries, and that further normalization based on reads per million (RPM) was not biased (Supplementary data, Fig. S4).   separated the PFOS10 group from all the other groups whereas we did not observe any separation along the first and the second principal component, indicating marginal effects in the other exposure groups. Furthermore, the confidence ellipses containing each group revealed that there is higher variation in PFOS10 and POP70 than in the Control, POP10 and PFOS70 groups. Venn diagrams showed that PFOS10 and PFOS70 share the highest number of overlapping DEGs with 54 common genes (Fig. 1F). Interestingly these two groups show a clear separation in the PCA analysis, since all the common genes are downregulated in PFOS10 and upregulated in PFOS70. Between POP10 and POP70 the number of mutually DEGs was much lower, only 5 genes (Fig. 1F). Finally, myhz1.1 was mutually differentially expressed in all groups (Supplementary File 2).

Pathway analysis with IPA and Webgestalt
DEGs were imported to Webgestalt using a custom background of all measured genes for gene enrichment analysis. Results of enriched pathways are summarized in Tables 1 and 2. The most enriched pathways were evident in the PFOS10 and PFOS70 groups. This result was expected due to the fact that these two groups had the highest number of DEGs. KEGG analysis revealed that the most enriched pathways in the PFOS10 group were involved in metabolism, peroxisome proliferator-activated receptors (PPAR) signaling (Supplementary data Fig.  S5) and extracellular matrix-receptor (ECM) interaction. Gene ontology (GO) analysis for biological functions revealed multiple pathways involved in lipid metabolic and homeostatic processes. KEGG analysis for PFOS70 revealed that significantly enriched pathways (p < 0.05) were involved in ECM-receptor interaction, cardiac muscle contraction and calcium signaling. GO analysis revealed multiple pathways involved in muscle contraction. Interestingly, albeit not significant, pathways involved in muscle contraction (both cardiac and skeletal) and calcium signaling were also observed in the remaining 2 groups POP10 and POP70 (Tables 1 and 2 and Supplemental File 3). Using these results, we used IPA to further explore the calcium signaling pathway. Datasets of PFOS70, POP10 and POP70 were overlaid ignoring the DE analysis cutoff values, to include all genes involved in the pathway in question, and the molecule activity predictor (MAP) function revealed activation of muscle contraction via increased influx of calcium through ryanodine receptors (RyR) (Fig. 2).

Validation with RT-qPCR
We performed a validation of the sequencing results by measuring 10 genes that were differentially expressed in either or all conditions. Validation of these 10 genes showed a significant correlation (p < 0.05) between qPCR and sequencing data (Supplementary data Fig. S6, r = 0.902 and Table S5).
3.6. Behavioral outcome of larvae exposed to RyR agonist and antagonist Based on the results of the transcriptomic analysis, we performed an experiment using the light-dark transition behavioral assay to uncover a potential mechanism of action of POPs on the calcium signaling pathway via ryanodine receptors (RyR). Caffeine, an established agonist, and dantrolene, an established antagonist of RyR were used (Hernández-Fonseca and Massieu, 2005). All details of the statistical results are summarized in Supplementary data Table S6.
In the acute (10 min before test) and chronic (6-96 h post fertilization, hpf) exposure scenarios, three groups were tested with or without the addition of 50 μM caffeine; Control, PFOS70, and POP70. For the controls, both the chronic and acute exposure to caffeine resulted in a significant decrease in the time spent active and an increase in swimming speed, but caffeine had no effect on the distance moved. Although the raw data and model selection suggested a possible interaction between PFOS and POPs with chronic caffeine on swimming speed and distance moved, the interaction was not significant for any endpoint (Supplementary data Fig. S7). However, all groups responded similarly to acute caffeine exposure with no tendencies towards any interaction effect (Supplementary data Fig. S7).
We followed the same exposure scenario as the caffeine experiments for the dantrolene chronic and acute exposure. Three groups were tested with or without the addition of 50 μM dantrolene, Control, PFOS70, and POP70. Acute and chronic dantrolene exposure resulted in a decrease in distance moved and time spent active. There was no Table 2 GO analysis of canonical pathways involved in biological processes for the exposed groups, FDR false discovery rate.  interaction on any endpoint following acute exposure to dantrolene (Fig. 3). On the contrary there was an interaction effect between PFOS and POPs following chronic dantrolene exposure (Fig. 3C). Here, whereas swimming speed was reduced in the control following dantrolene exposure (p-value < 0.01), there was no effect of dantrolene in either the PFOS or POPs exposed fish.

Chemical analysis
We used a variety of methods to determine chemical accumulation in zebrafish larvae tissue after 90 h of exposure (6-96 hpf). Results of the chemical analysis are summarized in Table 3. In general, the POP mixture analysis revealed that the chemicals with the highest accumulation in the zebrafish larvae were PFUnDA, PFDA, PFOS, HBCD, αchlordane, α-HCH, β-HCH and γ-HCH. For the single PFOS exposure, results showed an accumulation of 31.5% of the nominal input concentration in the tissue of larvae zebrafish (Table 3). The accumulation of PFOS incorporated in the POP mixture, in the tissue of larvae, was lower than the one from the single PFOS exposure and reached a level of 12.7% of the nominal input concentration.

Discussion
In the present study, we observed behavioral parameters in zebrafish larvae after exposure to a POP mixture or PFOS for 90 h, followed by transcriptome analysis to try and identify potential toxic mechanisms. Subsequently, we identified the calcium signaling pathway via the activation of RyRs as a potential mechanism of POPs and PFOS exposure and found interactions between POPs and PFOS with the RYR antagonist dantrolene on our behavioral endpoint. As such, although POPs and PFOS resulted in some differences in transcriptome profiles, a common mechanism may be behind behavioral aberrations.
The mixture of POPs used in the present study was based on the average levels of chemicals found in human blood of the Scandinavian population and includes 29 chemicals (Berntsen et al., 2017). Following chemical analysis, and based on the recovered percentage in larvae and media we estimated that the low concentration of PFOS was equal to 106.5 μg/L and the high concentration was equal to 731.3 μg/L. In comparison, the concentration of PFOS incorporated in the POP mixture was equal to 55.3 μg/L and 380.1 μg/L for POP10 and POP70 exposures respectively. These values are very relevant for human toxicity. For example, in human blood, PFOS serum levels range from 1   3. Behavioral responses following the light-dark transition test in larval zebrafish exposed to PFOS or a POP mixture and co-exposed to RyR antagonist dantrolene chronically (90 h) or acutely (10 min prior to the behavioral test). (A) Distance moved. (B) Time spent active. (C) Swimming speed. Statistics are from linear mixed effect models. Data shown are least square means ± 95% CI. Different lowercase letters indicate general exposure effects (i.e. Control vs POPs vs PFOS). Asterisk indicates an effect of dantrolene within exposure group (for the significant interaction model only).
to 10 μg/L in multiple surveys of general populations, and from 100 to 1000 μg/L for highly exposed populations (Kato et al., 2015;Olsen, 2015). In addition, our PFOS measurements within the larvae (52.3 ng) are similar to the levels found by Huang et al. (2010) in zebrafish exposed to 2000 μg/L between 0 and 96 hpf (≈38 ng). Results from the chemical analysis showed that the highest recovery based on tissue accumulation in larva for the POP mixture was 21% for PFUnDA, followed by PFDA at 16.5%, PFOS at 12.7%, HBCD at 11.8%, β-HCH at 14.7% and γ-HCH at 16% of the nominal input concentrations. For the single PFOS exposure tissue, we recovered 31.5% of the nominal exposure. It is unclear where the remaining quantities went, since the chemical analysis revealed that most of the remaining quantities of chemicals were only recovered to some extent (up to 35%) from the media. Xenobiotic metabolism by CYP450 enzymes is active by 72 hpf in zebrafish larvae but it is highly unlikely that zebrafish at this stage could metabolize the missing quantities of chemicals (Chen et al., 2012). Additionally, even though test vessels were pre-conditioned with exposure media some adsorption might have also occurred, for instance in the case for PFASs high affinity to polystyrene vessels has been already documented (Llorca et al., 2018). Finally, we cannot rule out interactions between chemicals in our mixture since the recovered level in larval tissue for the single PFOS exposure was more than double of what was recovered from the POP mixture. Nevertheless, the levels of POPs recovered from the chemical analysis are closer to human exposure levels than would be expected from the nominal input.
We tested the effect of exposure to the POP mixture or PFOS alone on the behavioral parameters of zebrafish larvae at 96 hpf using a light-dark transition test and a thigmotaxis assay. The increase in "wallhugging" identified by the thigmotaxis assay, and the increase in locomotor activity identified in the light-dark transition test, are both perceived to indicate increased anxiety (MacPhail et al., 2009;Ali et al., 2011;Schnörr et al., 2012). Our results matched those from a previous study in our lab , whereby the POP mixture and PFOS led to similar increases in swimming speed. Here, it is noted that Khezri et al. (2017) failed to identify a behavioral effect following exposures to separate sub-mixtures of the POP mixture that didn't contain PFOS. It was notable that although there was a twofold reduction in tissue concentrations, there was no difference in the behavioral response between the highest concentration of the POP mixture and PFOS. This could suggest that either the chemicals within the POP mixture had an additive or synergistic effect with PFOS on behavior, or that the behavioral assay is not sensitive enough to detect such relatively small differences in chemical burden. Nevertheless, our results are in agreement with previous studies where developmental PFOS exposure has been linked with hyperactivity in zebrafish (Huang et al., 2010;Spulber et al., 2014;Jantzen et al., 2016) and in mammalian models like mice (Johansson et al., 2008) where similar concentrations of PFOS were used.
Transcriptomics analysis revealed that groups exposed to PFOS10 (0.55 μM) and PFOS70 (3.83 μM) were most affected having 879 and 164 DEGs, respectively. The low number of DEGs in the POP10 and POP70 exposed groups (41 and 18 genes respectively) suggests that larvae exposed to the mixture do not show as strong of a response compared to our singular PFOS exposure. Whether this is through unknown additive, synergistic, or antagonizing interactions between the chemicals present in our mixture or the fact that larvae exposed to PFOS had a higher accumulation of PFOS in their tissue than the level found after exposure to the mix cannot be determined at present. Nevertheless, we observed mainly a downregulation of genes in PFOS10 and an upregulation of genes in PFOS70, but also high numbers of unique DEGs in PFOS10 (814 genes) and PFOS70 (101 genes). This observation agrees with the notion of non-monotonic response, which states that a chemical can produce a decrease in response relative to controls at a low dose and a higher response at higher doses or vice versa, also known as a U-shape or inverted U-shape response (Conolly and Lutz, 2004). In a previous study with adult zebrafish exposed to bisphenol A an inverted U-shape relationship was observed between exposure and expression of cytochrome P450 aromatase (cyp19a) (Molina et al., 2018). Additionally, low-dose effect studies with chemicals representing environmental concentrations may produce a greater response as seen in this study where more DEGs were detected in the low PFOS exposure and can affect different systems as shown by the number of unique number of DEGs and enriched pathways in the two conditions (Vandenberg et al., 2012). This needs to be taken into consideration and further studies could focus on the testing of multiple concentrations of PFOS and the effects on the transcriptome level. Focusing on the PFOS70, POP10 and POP70 groups, IPA analysis identified the calcium signaling pathway as one of the significantly affected pathways. Although our initial hypothesis was that this would be mostly linked to brain development and behavior, we surprisingly found limited pathways linked to behavior, but strong enrichments in genes related to muscle development. Genes involved in this pathway were significantly upregulated and included actin alpha cardiac muscle 1 (ACTC1), myosin heavy chain 4 (MYH4, zebrafish orthologue myhz1.1), myosin heavy chain 7 (MYH7, zebrafish orthologue myh7l), myosin light chain 4 (MYL4, zebrafish orthologue myl4), myosin light chain 7 (MYL7, zebrafish orthologue myl7), troponin C1 (TNNC1, zebrafish orthologue tnnc1a), troponin I2 (TNNI2, zebrafish orthologue tnni2a.1) and troponin T2 (TNNT2, zebrafish orthologue tnnt2a) with many of these genes, such as MYH4, MYH7, TNNC1, TNNI2, being commonly expressed in two or all three groups. In particular, MYH4 was expressed in all four groups but showed a downregulation in PFOS10. KEGG analysis confirmed the results of IPA that the calcium signaling pathway was one of the commonly affected pathways. Calcium is loaded in the endoplasmic reticulum (ER) through the activity of the sarcoendoplasmic reticulum Ca +2 ATPase (SERCA). Postsynaptic calcium signals in vertebrate skeletal and cardiac muscles are generated via mechanical coupling of specialized calcium channels, the ryanodine receptors (RyRs) (Brennan et al., 2005;Hernández-Fonseca and Massieu, 2005). Troponin and myosin genes provide instructions for making myosin and troponin proteins. The troponin and myosin complex are part of a structure called sarcomere, which is the basic unit of muscle contraction. The troponin complex, together with calcium, helps regulate contraction of the cardiac and skeletal muscle (Schiaffino and Reggiani, 1996). Activation of RyR from cerebral cortex microsomes and increased influx of calcium has been observed after exposure to chemicals, including PCBs, polybrominated diphenyl ethers, and PFOS (Dusza et al., 2018). Although our canonical pathway results point to the direction of calcium release by muscle specific RyR isoforms, the results of Dusza et al. (2018) seem to agree with our study. An upregulation of genes related to muscle function such as actins, myosins and tropomyosins has also been observed in another study where zebrafish were exposed to 1 mg/L of PFOS from 2 to 5 dpf (Martínez et al., 2019). The upregulation of these genes might also suggest a higher translation of muscle related proteins (actins, myosins and troponins) which in turn might relate to the longer distance moved and higher swimming speed of zebrafish larvae exposed to PFOS or POP mixture.
Based on the transcriptomic findings, we looked for treatment interactions with the RyR agonist caffeine and the RyR antagonist dantrolene. Results for the caffeine experiments failed to reveal an interaction effect between exposure to chemicals and caffeine, however, dantrolene exposure revealed an interaction with chemical exposure for swimming speed in chronically exposed larvae. Dantrolene can suppress intracellular calcium release from the sarcoplasmic reticulum by inactivating the RyR (Paul-Pletzer et al., 2002) and it has been shown to reduce the number of spontaneous muscular contractions in zebrafish embryos exposed to 10-1000 μM (Brennan et al., 2005). In controls, chronic dantrolene exposure led to a reduction in swimming speed that was not observed in the POPs or PFOS treated larvae. This could indicate that the chemicals used in this study have a higher affinity to RyR binding sites than dantrolene. Further work into endpoints within this pathway, such as mRNA transcription levels and protein, or Crispr knockouts that could pinpoint the exact mechanisms and genes that are responsible for these effects. would therefore be of interest.
Lipid metabolism, transport, and homeostasis were the main biological processes affected by the exposure to the lowest concentration of PFOS (0.55 μM), through the involvement of PPAR signaling pathway (Supplementary data Fig. S5). Transcriptomic changes in the genes involved in the PPAR signaling pathway may lead to inhibition of all lipid biological processes since all genes in our dataset were downregulated (Supplementary data Table S7). Previous studies have also proposed the mechanism of action for PFOS toxicity via an interaction with nuclear receptors associated with metabolic regulation like PPAR alpha (Lau et al., 2007;Shi et al., 2009;White et al., 2011) with dysregulation in lipid metabolism and homeostasis often observed after exposure to PFOS (Cui et al., 2017;Das et al., 2017;Martínez et al., 2019).
In this study we aimed to assess the effect of exposure to a complex POP mixture or single PFOS exposure to zebrafish larvae behavior and transcriptome response. Exposure to chemicals led to hyperactivity in zebrafish larvae and this might be correlated with the activation of muscle contraction observed during pathway analysis. Whilst further research is necessary, our results point to a possible involvement in muscle related toxicity via calcium signaling, rather than via neural related pathways. Lipid biological processes were also affected by PFOS exposure. Further studies should focus on the effect of PFOS on the lipid profile of zebrafish larvae to understand how dysregulation of lipid processes correlates with changes of lipid composition.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data
Supplementary data to this article can be found online at https:// doi.org/10.1016/j.envres.2020.109702.

Funding
This project received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Innovative Training Network (ITN) program [Grant agreement No. 722634].