Persistent Effects of Developmental Exposure to 17α-Ethinylestradiol on the Zebrafish (Danio rerio) Brain Transcriptome and Behavior

The synthetic estrogen 17α-ethinylestradiol (EE2) is an endocrine disrupting compound of concern due to its persistence and widespread presence in the aquatic environment. Effects of developmental exposure to low concentrations of EE2 in fish on reproduction and behavior not only persisted to adulthood, but have also been observed to be transmitted to several generations of unexposed progeny. To investigate the possible biological mechanisms of the persistent anxiogenic phenotype, we exposed zebrafish embryos for 80 days post fertilization to 0, 3, and 10 ng/L EE2 (measured concentrations 2.14 and 7.34 ng/L). After discontinued exposure, the animals were allowed to recover for 120 days in clean water. Adult males and females were later tested for changes in stress response and shoal cohesion, and whole-brain gene expression was analyzed with RNA sequencing. The results show increased anxiety in the novel tank and scototaxis tests, and increased shoal cohesion in fish exposed during development to EE2. RNA sequencing revealed 34 coding genes differentially expressed in male brains and 62 in female brains as a result of EE2 exposure. Several differences were observed between males and females in differential gene expression, with only one gene, sv2b, coding for a synaptic vesicle protein, that was affected by EE2 in both sexes. Functional analyses showed that in female brains, EE2 had significant effects on pathways connected to the circadian rhythm, cytoskeleton and motor proteins and synaptic proteins. A large number of non-coding sequences including 19 novel miRNAs were also differentially expressed in the female brain. The largest treatment effect in male brains was observed in pathways related to cholesterol biosynthesis and synaptic proteins. Circadian rhythm and cholesterol biosynthesis, previously implicated in anxiety behavior, might represent possible candidate pathways connecting the transcriptome changes to the alterations to behavior. Further the observed alteration in expression of genes involved in synaptogenesis and synaptic function may be important for the developmental modulations resulting in an anxiety phenotype. This study represents an initial survey of the fish brain transcriptome by RNA sequencing after long-term recovery from developmental exposure to an estrogenic compound.

The synthetic estrogen 17α-ethinylestradiol (EE 2 ) is an endocrine disrupting compound of concern due to its persistence and widespread presence in the aquatic environment. Effects of developmental exposure to low concentrations of EE 2 in fish on reproduction and behavior not only persisted to adulthood, but have also been observed to be transmitted to several generations of unexposed progeny. To investigate the possible biological mechanisms of the persistent anxiogenic phenotype, we exposed zebrafish embryos for 80 days post fertilization to 0, 3, and 10 ng/L EE 2 (measured concentrations 2.14 and 7.34 ng/L). After discontinued exposure, the animals were allowed to recover for 120 days in clean water. Adult males and females were later tested for changes in stress response and shoal cohesion, and whole-brain gene expression was analyzed with RNA sequencing. The results show increased anxiety in the novel tank and scototaxis tests, and increased shoal cohesion in fish exposed during development to EE 2 . RNA sequencing revealed 34 coding genes differentially expressed in male brains and 62 in female brains as a result of EE 2 exposure. Several differences were observed between males and females in differential gene expression, with only one gene, sv2b, coding for a synaptic vesicle protein, that was affected by EE 2 in both sexes. Functional analyses showed that in female brains, EE 2 had significant effects on pathways connected to the circadian rhythm, cytoskeleton and motor proteins and synaptic proteins. A large number of non-coding sequences including 19 novel miRNAs were also differentially expressed in the female brain. The largest treatment effect in male brains was observed in pathways related to cholesterol biosynthesis and synaptic proteins. Circadian rhythm and cholesterol biosynthesis, previously implicated in anxiety behavior, might represent possible candidate pathways connecting the transcriptome changes to the alterations to behavior. Further the observed alteration in expression of genes involved in synaptogenesis and synaptic INTRODUCTION 17α-ethinylestradiol (EE 2 ) is the most potent endocrine disrupting compound (EDC) pollutant contaminating the aquatic environment (Aris et al., 2014;Laurenson et al., 2014). Being 10-fold more active than estradiol itself (Thorpe et al., 2003;Denny et al., 2005), and having a predicted no adverse effect level of 0.1 ng/L for water-living organisms (Caldwell et al., 2012), its ubiquitous worldwide presence in sewage treatment plant effluents and reclaimed water is a matter of great concern (Mompelat et al., 2009;Aris et al., 2014). It is detected in effluents in levels from <1 ng/L (the present detection limit of such analyses) to as much as 2-300 ng/L (Ternes, 1998;Kolpin et al., 2002;Laurenson et al., 2014). EE 2 is persistent in the environment and signs of bio-magnification have been observed (Aris et al., 2014).
In addition to the effects of EDCs on the reproductive axis in fish we and others have reported a significant influence also on non-reproductive behaviors in adult fish of importance for fitness, such as risky behavior, aggression, anxiety, and shoaling (Espmark Wibe et al., 2002;Majewski et al., 2002;Bell, 2004;Colman et al., 2009;Xia et al., 2010;Hallgren et al., 2011;Reyhanian et al., 2011;Dzieweczynski et al., 2014;Heintz et al., 2015). In recent studies, we found that low-dose EE 2 exposure of both zebrafish and guppies during development resulted in increased anxiety as adults (Volkova et al., 2012(Volkova et al., , 2015b, suggesting irreversible changes of the neuroendocrine system. We have also found that the effects on anxiety behavior from developmental exposure of EE 2 was transmitted to the first and the second generation of unexposed progeny in both guppy and zebrafish (Volkova et al., 2015a,b). Effects of EDCs at the level of tissue organization in the developing brain has in mammals been observed to affect fertility and behavior (McLachlan, 2001;Crews and McLachlan, 2006;Dickerson and Gore, 2007). Altered non-reproductive behavior in rodents developmentally exposed to EDC's have been shown to be transmitted over generations of unexposed progeny accompanied by alterations in the brain transcriptome (Skinner et al., 2008;Wolstenholme et al., 2012). Alterations in the epigenome, shown to be target of hormone modulations (Simerly et al., 1997;McCarthy et al., 2009;Casati et al., 2015;McCarthy and Nugent, 2015), are indicated as mediators of the changes in brain imprinting. Epigenetic mechanisms are dynamic and influenced by the environment, especially Abbreviations: EDC, Endocrine disrupting compound; EE 2 , 17α-ethinyl estradiol; NT, novel tank; RNAseq, RNA sequencing; qPCR, quantitative Real Time polymerase chain reaction; ncRNA, non-coding RNA; miRNA, micro RNA. during the early life stages of an organism (Casati et al., 2015).
In this study, we analyze if persistent effects on gene expression in the zebrafish brain could be discerned in the adult fish after developmental EE 2 exposure, accompanying an anxiogenic phenotype. We postulated that effects of anxiety should affect the regulation of genes related to the stress axis in the brain, but presumably also several other functions. To achieve this, we repeated the study mentioned above (Volkova et al., 2015b) with slightly different exposure protocol and, as it turned out, higher actual exposure levels. We analyzed three non-reproductive behaviors after 120 days remediation, and performed RNA sequencing of the complete brain transcriptome in male and female zebrafish. The purpose was to extend the understanding of the persistent effects of estrogenic substances that linger after discontinued exposure.

Animals and Treatments
Fertilized zebrafish eggs from eight separate parental pairs of the wild type strain AB were obtained from the Karolinska Institute Zebrafish Core Facility, Solna, Sweden. Eggs from each parental pair were divided into three lots and assigned to treatment groups of 0, 3, and 10 ng EE 2 /L, respectively. EE 2 (17αethinylestradiol, Sigma-Aldrich) was dissolved in acetone and stock solutions were mixed with pre-heated fish maintenance system water to final nominal concentrations of 0, 3, and 10 ng EE 2 /L. The nominal concentrations were considered relevant as they are within the range of concentrations measured in effluents and surface waters (Ternes, 1998;Kolpin et al., 2002;Laurenson et al., 2014). The final concentration of acetone was 1 ppm in control and EE 2 solutions. The embryos were exposed for 80 days (1-80 days post fertilization) in 2 L tanks in a flow-through system with an exchange rate of 5-6 ml/h. All solutions in the flow-through system were pre-mixed and exchanged daily. Concentrated paramecia culture was fed to each tank once a day for 10 days. Artemia was added to the diet once a day from day 10. Sera Dry Flakes (Vipan, Germany) were added to the diet twice daily from week 6. Exposure was terminated at 80 dpf and the fish transferred to clean water were they were kept under normal maintenance conditions until adulthood, resulting in a 120 day remediation period. Sexes were separated based on secondary sexual characteristics at 4 months of age and thereafter re-checked weekly for changes. Animals were kept in a 12/12 h light/dark cycle at 25-27 • C, pH 7.0. All experiments and fish handling were performed according to the Swedish Animal Care legislation and approved by the Southern Stockholm Animal Research Ethics Committee (Dnr. S130-09).

EE 2 Concentrations
Water samples were collected from different tanks at nine separate occasions during the exposure and stored in darkness at −20 • C until analysis. EE 2 concentrations were analyzed in single or duplicate samples as previously described in Volkova et al. (2015b). Briefly, water samples (100 mL) were extracted on 100 mg Strata X-33µ Polymeric Reversed Phase cartridges, reconditioned with MeOH. Dionex ultimate 3000 LC system (Thermo Scientific, San Jose, CA, USA), coupled to a triple quadruple mass spectrometer (TSQ Vantage, Thermo Scientific, San Jose, CA, USA) was used for quantitation of EE 2 content. The quantification range was 0.5-100 ng/L, using EE 2 -d4 as internal standard.

Behavior Studies
The behaviors were analyzed in the novel tank test (NT; Egan et al., 2009; Figure S1), shoaling test (Moretz et al., 2007; Figure S1), and scototaxis test (Maximino et al., 2010b; Figure S1), as previously described (Volkova et al., 2015b). Briefly, The NT and shoaling tests were performed one after the other in the same test episode. At the right end of the test tank (20 × 20 × 40 cm) filled with 15 L pre-heated tap water, a transparent Plexiglas screen trapped five untreated male or female littermates. A black plastic sheet prevented visual contact with the test compartment. One horizontal and one vertical middle line divided the tank into upper/bottom and right/left halves. The NT test was initiated by introducing a fish to the test tank, and latency time before the first crossing of the horizontal line to the upper half, number of transitions to the upper half and total time in the upper half was recorded. Swimming activity was quantified as number of lines crossed a grid, both horizontally and vertically, during 60 s during the last minute of the NT test. Behavior was videomonitored for 5 min, after which the black screen was removed, revealing the hidden shoal. The shoaling test started when the test fish made contact with the shoal, and latency to leave the shoal (cross the vertical line), number of transitions and time spent in the opposite half of the tank was recorded for 5 min. Fish that did not make contact with the shoal within 5 min were excluded from the analyses.
In the scototaxis test, the test tank (20 × 20 × 40 cm) was divided into one black and one white half and filled with preheated tap water up to 10 cm. The tank had two transparent central sliding doors, creating a compartment of 5 × 20 cm. The test fish was introduced into the central compartment, and after a 5 min habituation period the sliding doors were raised and the behavior of the fish was video recorded from above for 5 min. Latency to first entrance into the white half, number of transitions to white half and total time spent in white half was recorded.
All tests were video recorded and manually analyzed. Behavior analyses were not blinded for treatment due to logistical reasons. Three tanks were operating in parallel, and the behavior tests were performed between 9:00 a.m. and 1:00 p.m.

Dissection and Sex Verification
Fish previously macroscopically determined as males or females were re-examined based on secondary sexual characteristics, and separated before behavior testing. After behavior testing, two fish of each sex from each family group and treatment group were weighed and then sacrificed by anaesthetization in 0.5‰ 2-phenoxyethanol (Sigma-Aldrich) followed by immediate decapitation. During dissection, gonads and livers were removed and weighed and Gonadosomatic index (GSI) and hepatosomatic index (HSI) were calculated by the formula: organ weight/fish weight × 100. Brains were removed and stored at −80 • C in RNA later (Sigma-Aldrich).

Ion Proton RNA Sequencing
Based on the results in the novel tank test, RNA sequencing analysis was selected to be performed on control male brains and brains from males exposed to 3 ng/L EE 2 , while for females, brains were from unexposed females and females exposed to 10 ng/L EE 2 . Three biological replicates were used in the RNA sequencing analysis for each treatment group and sex. In the female analysis, fish from the same family groups were represented in both the control and treated samples; meaning that the control females and exposed females were siblings. Male samples were however taken from different family groups and the exposed males were therefore not the siblings of the control males. Brains were homogenized and total RNA extracted with TriReagent (Sigma-Aldrich) according to the manufacturer's protocol. Total RNA from selected samples was submitted to Bea Core Facility (Karolinska Institute, Huddinge, Sweden) for quality checking with R6K Screen Tape. IonProton RNA sequencing was performed by SciLifeLab (Uppsala University, Sweden) collecting ∼40 million reads per sample.
The RNA was treated with Ribo-Zero TM rRNA Removal Kit (Epicentre/Illumina) to remove ribosomal RNA, and purified using Agencourt RNAClean XP Kit (Beckman Coulter). The RNA was then treated with RNaseIII according to the Ion Total RNA-Seq protocol (Life Technologies) and purified with Magnetic Bead Cleanup Module (Life Technologies). The size and quantity of RNA fragments were assessed on the Agilent 2100 Bioanalyzer system (RNA 6000 Pico kit, Agilent) before proceeding to library preparation, using the Ion Total RNA-Seq kit (Life Technologies). The libraries were amplified according to the protocol and purified with Magnetic Bead Cleanup Module (Life Technologies). Samples were then quantified using the Agilent 2100 Bioanalyzer system (Agilent) and Fragment Analyzer (Advanced Analytical) and pooled followed by emulsion PCR on the Ion OneTouch TM 2 system using the Ion Proton TM Template OT2 200 v3 Kit (Life Technologies) chemistry. Enriching was conducted using the Ion OneTouch TM ES (Life Technologies). Samples were loaded on an Ion PI v2 Chip (2 samples per chip) and sequenced on the Ion Proton TM System using Ion Proton TM Sequencing 200 v3 Kit (200 bp read length, Life Technologies) chemistry.

Bioinformatics and Biostatistics
The software FastQC (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/) was used to assess read quality, and after that data was mapped to the version 9 of the Danio rerio genome using the software Star (Dobin et al., 2013). The mapped reads were converted to count data with the script htseq-count (Anders et al., 2015) using the ensembl annotation of the zebrafish genome sequence. Statistical modeling of gene expression was done with the library edgeR (Robinson et al., 2010) in R (McCarthy et al., 2012) following the workflow suggested by the authors. Before the modeling was done all genes that had expression levels lower than 1 count per million mapped reads in at least three different libraries were omitted from the data set. In the detection of differentially expressed genes family structure were taken into account for the female samples. False discovery rates (FDR) corrected p-values were estimated with a Benjamini-Hochberg procedure. After this correction 175 genes had p-values lower than 0.05 and were regarded as differentially expressed.

Functional Analysis
Classification of genes and predictions of biological gene function were performed manually, due to the observed low detection capacity of several automated classifications tested. Manual ortholog search was done in Ensembl. GO-terms were found in Zfin (for zebrafish), Entrez gene (for human and rodent orthologs) and NGNC (human orthologs).

Quantitative Real-Time PCR (qPCR)
For verification of the RNAseq results, the expression of selected differentially expressed genes in the male brain were analyzed with qPCR. RNA isolation was performed with individual brains from 10 control males and 10 males exposed to 3 ng/L, homogenized in TriReagent (0.8 ml/sample, Sigma-Aldrich, Germany) according to the manufacturer, quantified using NanoDrop ND-1000 spectrophotometer, and RNA quality verified by the 260/280 nm absorption ratio. Reverse transcription was performed from Quantitect R Reverse Transcription Kit (Qiagen) from 1 µg RNA per sample and then diluted 1:10 for qPCR reactions. Bio-Rad C1000 Touch TM Thermal cycler, CFX96 TM Real-Time System was used for quantitative Real-Time PCR (qPCR). qPCR reactions were run in triplicates and contained 5 µL cDNA template (25 ng RNA), 2.5 µL each of forward and reverse primer and 10 µL iTaq TM Universal SYBR R Green Supermix (BIO-RAD), using cycling parameters as recommended by the manufacturer. Oligonucleotide primers ( Table 1) were designed with Primer-Blast primer designing tool (Ye et al., 2012). Normalized expression was calculated for quantification according to the 2 − CT method (Schmittgen and Livak, 2008) with two internal control genes, 18S rRNA and Elongation factor 1 alpha (Elf1α), genes found suitable for zebrafish tissue analysis (Tang et al., 2007).

Statistical Analysis
Data from the behavior tests, GSI and HSI, and qPCR was analyzed in the statistical package R (R Core Team, 2013) with generalized linear mixed-effects models using package lme4 (Bates et al., 2015) and followed by post-hoc tests where necessary using the multcomp package (Hothorn et al., 2008). To separate maternal family effects from effects of treatment and sex, family
*Designed by Primer-BLAST primer designing tool.
was used as random factor as several individuals from the same family where represented at different treatment levels. Sex, treatment and treatment × sex where included as fixed factors. All models were analyzed using LR-tests. In cases were females and males were analyzed separately we used treatment as a fixed factor and family as a random factor. For models with heteroscedastic residuals the response variable was log transformed. Data in the form of counts were analyzed assuming Poisson distribution of residuals.

EE 2 Concentrations
The actual concentration of EE 2 was determined from water samples from 9 occasions during the 80 days exposure period. The determined actual concentrations of EE 2 (mean ± SEM) in samples taken from the treatment tanks were 2.14 ± 0.33 and 7.34 ± 1.42 ng/L for 3 and 10 ng/L, respectively. This represents 71 and 73% of the nominal concentrations. Control water samples contained no detectable levels (<0.5 ng/L) of EE 2 .

Body Weight, GSI, and HSI
Of the dissected animals all fish phenotypically classified as females contained ovaries, and all fish phenotypically classified as males had testes (N males = 47, N females = 46). After 120 days of remediation in clean water there were no differences in body weight between control fish and treated fish (Chisq = 5.16, p = 0.08). Zebrafish males treated with 3 ng/L EE 2 had significantly lower GSI compared to control males (z = −2.86, p = 0.012), while no significant difference was observed in GSI between controls and males treated with 10 ng/L (z = −0.214, p = 0.98). However, males treated with 10 ng/L had significantly higher HSI compared to control males (z = 3.015, p = 0.007) whereas males treated with 3 ng/L did not (z = 1.66, p = 0.22; Table 2). No significant differences in GSI or HSI were found between treated and control females ( Table 2).

NT Test
There were significant interactions between treatment and sex for latency to upper half (Chisq = 7.48, p = 0.024), and number of transitions to upper half (Chisq = 14.17, p < 0.001). There were no significant interactions found for the time spent in the upper half but it was affected by both treatment (Chisq = 13.77, p = 0.001) and sex (Chisq = 12.77, p < 0.001; Figure 1, Table 3).
To further evaluate the effects, we analyzed males (N Control = 40, N 3 ng = 39, N 10 ng = 39) and females (N control = 38, N 3 ng = 40, N 10 ng = 39) separately (Figure 1). Females exposed to 3 ng/L spent significantly less time in the upper half (z = −2.96, p = 0.009). Females exposed to 10 ng/L had significantly increased latency to upper half (z = 2.39, p = 0.045), fewer transitions to upper half (z = −3.0, p = 0.008) and spent significantly less time in upper half (z = −3.51, p = 0.001). Males exposed to 3 ng/L EE 2 had significantly increased latency to upper half (z = 2.50, FIGURE 1 | Novel tank behavior in male and female zebrafish developmentally exposed to 0, 3, or 10 ng/L EE 2 and remediated in clean water for 120 days. (A) Latency time before crossing the horizontal midline to the upper half. (B) Number of transitions to upper half and (C) Total time spent in upper half. Data represent mean ± 95% CI, # p < 0.05, *p < 0.05, **p < 0.01. *Significantly different from control fishes of the same sex, # control females significantly different from control males. p = 0.034). No significant treatment effects were observed in the novel tank test for males developmentally treated with 10 ng/L compared to the control males (Figure 1).
Swimming activity measured during the last minute of the novel tank test revealed average ± SD for females: 70 ± 21, 66 ± 19, 69 ± 15 for 0, 3, and 10 ng/L, respectively, and average ±SD for males: 74 ± 25, 70 ± 16 and 69 ± 14 for 0, 3, and 10 ng/L, respectively. The swimming activity was not significantly different between the treatments or between males and females.

Shoaling Test
Fish that failed to make contact with the group within 5 min were excluded from the analyses which slightly lowered the number of observations (N Control = 57, N 3 ng = 58, N 10 ng = 51). Twenty-one (8 female and 13 male fish), 21 (13 female Analyses include zebra fish males and females treated with EE 2 during development followed by 120 days remediation. Behaviors were analyzed in the novel tank test, shoaling test and scototaxis test. *p < 0.05, **p < 0.01, ***p < 0.001. and 9 male fish), and 27 (24 female and 3 male fish) fishes where excluded for this reason in the 0, 3, and 10 ng/L treatment groups, respectively. EE 2 exposure significantly decreased the number of transitions made away from shoal (Chisq = 24.8, p < 0.001). There were also significant differences by sex in latency to first transition away from shoal (Chisq = 16.3, p < 0.001), number of transitions away from shoal (Chisq = 202.5, p < 0.001) and total time away from shoal in the opposite half (Chisq = 45.4, p < 0.001; Table 3). Even if there were no significant interactions we also analyzed the two sexes separately; males (N Control = 27, N 3 ng = 30, N 10 ng = 36) and females (N Control = 30, N 3 ng = 28, N 10 ng = 15; Table 3). Effect of treatment was found only in females. Number of transitions away from the shoal was decreased in females treated with 3 ng/L (z = −4.38, p < 0.001) and 10 ng/L (z = −3.57, p < 0.001; Figure 2).
FIGURE 2 | Shoaling behavior in male and female zebrafish developmentally exposed to 0, 3, or 10 ng/L EE 2 and remediated in clean water for 120 days. (A) Latency time before leaving shoal (crossing the vertical midline), (B) Number of transitions to the opposite half and (C) Total time spent in the opposite half. Data represent mean ± 95% CI, # p < 0.05, ***p < 0.001. *Significantly different from control fishes of the same sex, # control females significantly different from control males.

Scototaxis Test
The models including both males and females (N control = 69, N 3 ng = 68, N 10 ng = 79) reveal significant interactions between sex and treatment in number of entries into white half (Chisq = 46.2, p < 0.001) and total time spent in white half (Chisq = 9.28, p = 0.0097; Table 3). Latency was not significantly affected by treatment but it was very close (Chisq = 5.75, p = 0.057) however there was a significant difference between males and females (Chisq = 21.3, p < 0.001; Table 3).
Analyzing males (N Control = 39, N 3 ng = 38, N 10 ng = 40) and females (N Control = 30, N 3 ng = 29, N 10 ng = 39) separately (Figure 3) the results showed that both sexes were affected of treatment. Females exposed to 3 ng/L EE 2 had increased latency to first transition into white half (z = 2.41, p = 0.042) and a significantly decreased total number of entries into the white half (z = −4.77, p < 0.001). No effect was observed in females treated FIGURE 3 | Scototaxis behavior in male and female zebrafish developmentally exposed to 0, 3, or 10 ng/L EE 2 and remediated in clean water for 120 days. (A) Latency time before entering the white half, (B) Number of transitions to the white half, and (C) Total time spent in the white half. Data represent mean ± 95% CI, # p < 0.05, *p < 0.05, **p < 0.01, ***p < 0.001. *Significantly different from control fishes of the same sex, # control females significantly different from control males.
with 10 ng/L EE 2 compared to control females. Males exposed to 3 and 10 ng/L EE 2 spent significantly less time in the white half compared to control males (z = −2.90, p = 0.011) and (z = −3.11, p = 0.005), respectively. While males exposed to 10 ng/L made significantly fewer entries into the white half (z = −4.59, p < 0.001) compared to control males.

RNA Sequencing
Of the 33,373 gene features that were annotated in genome sequence, 18,792 in the females, and 18,445 in males had gene expression levels high enough to be retained for analysis. One hundred and forty-six sequences, 95 protein-encoding, and 51 non-coding sequences, were observed to be differentially expressed in brains from developmentally EE 2 -exposed zebrafish after recovery, compared to brains from unexposed controls. In the brain samples from males, 34 coding genes showed significantly different expression levels between control males and males exposed to 3 ng/L EE 2 . For the female brains, 62 coding genes were differentially expressed due to the developmental exposure to 10 ng/L EE 2. Only one gene was affected by developmental EE 2 exposure in both male and female brains: sv2b, coding for a synaptic vesicle protein, was upregulated by EE 2 (log fold change: 0.86) in males, while in females the same gene was down-regulated (log fold change: −1.2). In males, another copy of this gene, svbb, located upstream of sv2b on chromosome 25 was down-regulated by EE 2 (log fold change: −0.68). No effect on svbb was detected in females. The lists of differentially expressed protein coding genes are shown in Table 4 (females) and Table 5 (males). In addition, 48 differently expressed sequences in females and three in males were annotated as non-coding RNA (ncRNA; Table 5).

Functional Analysis of Protein Coding Genes
Functional analyses were performed for males and females separately. As several automated classifications tested yielded low percentage of genes classified to function, a manual search was performed based on information from human and rodent ortholog data offered through Entrez gene and NGNC databases in addition to the zebrafish database Zfin. This search identified 73% of all differentially expressed genes (45 out of 62 genes for females and 25 out of 34 genes for males). Differentially expressed genes in Tables 4, 5 are presented under the categories that were manually assigned. For female brains, the search identified "Circadian Rhythm" and "Cytoskeleton and motor proteins" as top putative biological pathways, with seven affected genes each ( Table 4). "Cholesterol biosynthesis" was identified as the top putative pathway in male brains, with four differentially expressed genes ( Table 5).

Circadian Rhythm
A manual homolog search in Ensembl identified a total of seven differentially expressed genes encoding proteins presumably involved in the circadian rhythm in brains from developmentally EE 2 -exposed compared with unexposed female brains ( Table 4). The genes were dbpb (log fold change: 0.92), nr1d1 (log fold change: 2.21), per1a (log fold change: 1.83), ciarta (log fold change: 1.99), ciartb (log fold change: 0.81), cipca (log fold change: 1.65), and bhlhe41 (log fold change: 1.01). None of the differentially expressed genes in male brains were associated with the circadian rhythm.

Cholesterol Biosynthesis and Transportation
Cholesterol biosynthesis was identified as the top putative pathway affected by developmental EE 2 -exposure in zebrafish male brains after 120 days of remediation (Table 5). Four genes, sqlea (log fold change: 1.84), lss (log fold change: 1.34), msmo1 (log fold change: 0.98), and hmgcs1 (log fold change: 1.01) were significantly upregulated by EE 2 . None of these genes were significantly altered by EE 2 exposure in female brains. However, the manual search identified three genes involved in binding and transporting cholesterol ( Table 3) which were significantly down-regulated in brains of females developmentally exposed to EE 2 compared with brains from control females: apoeb (log fold change: −1.12), abca12 (log fold change; −3.08), and prom2 (log fold change: −1.96).

Synapses
Several genes involved in synaptogenesis and synapse function were affected by developmental EE 2 exposure in both male and female zebrafish brains. In male brains, in addition to the altered expression of the two synaptic vesicle genes mentioned above, sv2b (log fold change: 0.86) and svbb (log fold change: −0.68), also a decreased expression of arhgef9b (log fold change: −1.05) encoding collybistin was observed in response to developmental EE 2 exposure ( Table 5). Females developmentally exposed to EE 2 showed a down-regulation of brain sv2b expression (log fold change: −1.20) as mentioned above, and also of the expression si:ch211-120g10.1 (log fold change: −1.93) ( Table 4). In addition, females showed an upregulation of brain rtn4rl2a (log fold change: 1.08) and bai1b expression (log fold change: 0.82), coding for a negative regulator of angiogenesis in the brain, in response to EE 2 .

Cytoskeleton and Motor Proteins
Developmental EE 2 exposure affected seven genes related to the cytoskeleton and motor proteins in zebrafish female brains ( Table 4). The expression of four axonemal dynein genes were significantly down-regulated; dnah3 (log fold change: −1.51), dnah8 (log fold change −1.86), dnah12 (log fold change: −1.10) and AL935046.1, a suggested ortholog of the human DNAH2 (log fold change: −1.38). Down-regulated expression of two tubulin genes, tuba8l (log fold change: −0.90) and zgc:55461, a suggested ortholog to the human TUBB4A and TUBB4B genes (log fold change: −1.48), was also detected. Expression of the myosin complex-encoding gene myhb (log fold change: 1.67) was upregulated in the brains of females developmentally exposed to EE 2 . No genes related to cytoskeleton or motor proteins were found to be affected in male brains.

Heme Metabolism and Cardiovascular System
Several genes coding for heme-binding proteins were affected by developmental EE 2 treatment in both male and female zebrafish. In male brains, EE 2 exposure during early life resulted in an upregulation of hbaa1 (log fold change: 1.48) and si:ch211-5k11.6 (log fold change: 2.79) expression (Table 5); both gene products are predicted to be part of the hemoglobin complex. In female brains there was an upregulation of the mb gene expression (log fold change: 3.13) encoding myoglobin, by EE 2 (Table 4). Further, a down regulation of fech expression, encoding the terminal enzyme of the heme biosynthesis pathway (log fold change: −0.94) as well as of hmox2a (log fold change: −0.72) was observed. Also bai1b expression (log fold change: 0.82), a negative regulator of angiogenesis, was upregulated after developmental EE 2 exposure followed by 120 days of remediation (Table 4), while in males, the expression of tspan5b, coding for a protein suggested to be involved in angiogenesis, was down-regulated (log fold change: −2.02, Table 5).
Frontiers in Behavioral Neuroscience | www.frontiersin.org TABLE 4 | List of protein coding genes significantly affected in brains of female zebrafish developmentally exposed to 10ng/L EE 2 compared to non-exposed fish after 120 days of recovery in clean water.   5 | List of protein coding genes significantly affected in brains of male zebrafish developmentally exposed to 3ng/L EE 2 compared to non-exposed fish after 120 days of recovery in clean water.

Immune System
The brain expression of several genes in different ways involved in the immune response and inflammation were found to be affected by developmental exposure to EE 2 . The macrophage associated lectin lgals3bpb (log fold change: 1.17) gene expression was upregulated in female brains by the exposure, while the expression of socs3b (log fold change: −1.21), b2m (log fold change: −0.80), and ftr02f (log fold change: −3.20) were downregulated (Table 4). In males, downregulation of the expression of ptgr1 (log fold change: −2.32), encoding a protein involved in the inactivation of the pro-inflammatory factor leukotriene B4, and an upregulation of serpinb1l2 (log fold change: 5.71) expression was observed in zebrafish male brains by EE 2 ( Table 5).

DNA Repair and Chromatin Remodeling
The brain expression of two genes involved in DNA repair and recombination, brca2 (log fold change: 1.20) and smc1a (log fold change: 1.01), were upregulated in the brains of zebrafish males developmentally exposed to EE 2 ( Table 5). These genes were not affected in developmentally exposed female zebrafish. Females showed a marked down-regulation of brain cbx8a (log fold change: −4.59) and also of sirt7 (log fold charge: −1.34) expression by EE 2 (Table 4), encoding two proteins believed to be involved in chromatin remodeling and DNA repair.

Non-coding Sequences
Of the differentially expressed sequences 48 of the 110 in the females and 3 of 34 in males were annotated as non-coding RNA (ncRNA). A significant fraction, 19 of the 48 ncRNA in the female brain were novel microRNAs (miRNAs) whilst none of the 3 in males were miRNA ( Table 6).

Verification of Differential Expression by qPCR
Verification of the differential gene expression observed by RNA sequencing analysis was performed by qPCR for four selected target genes in the cholesterol biosynthesis pathway that was differentially upregulated by 3 ng/L EE 2 in male brains. Ten brain samples each in exposure and control groups were analyzed. As shown in Figure 4, significant increased expression was observed for msmo1 (Chisq = 9.00, df = 1, p = 0.003), sqlea (Chisq = 6.41, df = 1, p = 0.011), and lss (Chisq = 4.94, df = 1, p = 0.026). The expression of hmgsc1 was however not significantly higher in brains from EE 2 -exposed males than in control brains (Chisq = 1.86, df = 1, p = 0.173).

DISCUSSION
The present study verifies and extends our previous findings (Volkova et al., 2015b) that developmental exposure of zebrafish to EE 2 increase anxiety and shoal cohesion as adults. We also present an analysis of the whole brain transcriptome in exposed and control fish by RNAseq in an attempt to understand the mechanisms of the persistent behavior phenotype. The transcriptome analysis revealed marked differences in gene expression in males and females, with 94 protein coding genes significantly affected by the developmental exposure to EE 2 (Tables 4, 5). We used the same nominal concentrations as in Volkova et al. (2015b), but the actual measured concentrations were higher in the present study due to differences in exposure protocol. Although the experimental procedures were not identical, and the source of the AB strain was different among the studies, the similarities in the behavior results are obvious. Clear anxiogenic effects of EE 2 on behavior were observed in both the NT test and scototaxis test (Figures 1, 3). Both tests have been characterized by means of pharmacologic drugs in zebrafish (Maximino et al., 2015) however, only when exposed as adults. Scototaxis and NT have been suggested to be complementary rather than identical tests for anxiety based on differences in stimuli (Maximino et al., 2010a;Blaser and Rosemberg, 2012). Also, novelty in the NT has been suggested to be a conditional stimulus, and as such might be sensitive to minor external factor fluctuations, as discussed in Maximino et al. (2010a) and Blaser and Rosemberg (2012).
The current study verified that developmental EE 2 exposure of zebrafish increases anxiety as adults, after a long remediation 6 | List of non-coding genes significantly affected in brains of zebrafish developmentally exposed to 3 ng/L (males) and 10 ng/L (females) EE 2 compared to non-exposed fish after 120 days of recovery in clean water.  period in clean water. Increased anxiety is well-established as an EDC effect in mammals (Dugard et al., 2001;Ryan and Vandenbergh, 2006;Gioiosa et al., 2007;Skinner et al., 2008;Gonçalves et al., 2010), although the outcome differs by sex, age at treatment, exposure dose, and length (Gioiosa et al., 2013). Anxiety due to EE 2 has been observed by us and others in several fish species Reyhanian et al., 2011;Heintz et al., 2015). The behaviors affected are of great importance for fitness of fish populations in the wild, affecting survival, and predator avoidance, but also mating success and foraging.
In the shoaling test, developmentally exposed fish ventured away from the shoal less in agreement with previous results (Volkova et al., 2015b). The effects of EDCs on shoaling behavior are unclear at present, both increased, decreased, and unaffected shoal cohesion have been shown in EDC-exposed fish (Espmark Wibe et al., 2002;Ward et al., 2006;Xia et al., 2010;Reyhanian et al., 2011). A stress component in the behavior is expected, as tight shoaling is an anti-predator strategy, but no characterization of the test system by means of anxiogenic and anxiolytic drugs have been performed, however such studies are needed.
To learn more about what changes in the brain could be discerned to accompany the behavior phenotype, we performed RNA sequencing of the brain transcriptome. We acknowledge that this transcriptome analysis is hampered by the few biological replicates, and conclusions from the data have to be drawn with caution. Also, the use of whole brains might hide region-specific alterations. Transcriptome data are challenging to interpret, and we can only speculate on connections between identified pathways and imprinting of increased stress sensitivity due to developmental EE 2 exposure. However, RNA sequencing analyses of zebrafish brain are still few, and this study represents, to the best of our knowledge, the first study of estrogenic effects induced during development and persisting as a behavior phenotype in the adult fish. Further studies are needed to support, or contradict, the current findings.
No effects were observed in genes involved in the stress axis, neither on CRH, related genes in the hypothalamo-pituitaryinterrenal (HPI) axis nor the monoamine system shown to be involved in their regulation. We did, however, identify some possible candidate pathways that might indirectly affect anxiety behavior. It is thus possible that the alterations in anxiety behavior induced during development are not mediated via the most obvious target genes in the brain. The connection between estrogens, stress and the regulation of HPI axis and monoaminergic system in fish is well-established (Winberg and Nilsson, 1993;Pottinger et al., 1996;Winberg et al., 1997;Carpenter et al., 2007;Egan et al., 2009). Previous studies observed that acute EE 2 exposure result in altered expression of genes involved in the glucocorticoid receptor signaling pathway in the pituitary of sub-adult female Coho salmon (Harding et al., 2013) as well as in the hepatic tryptophan metabolism pathway, a precursor to serotonin, in adult zebrafish (Wit et al., 2010). The fish in this study were, however, not under EE 2 exposure, nor were they subjected to stressors prior to tissue collection, thus decreasing the likelihood of acute effects on the stress axis. The current study does, however, not support that direct effects on stress signaling underlie the adult anxiety phenotype observed in this study.
Only one gene, sv2b, was affected in both sexes of the treated animals. No other genes were significantly altered in both sexes after recovery from developmental EE 2 exposure, though some similarities between the affected pathways exist. Both males and females had alterations in pathways related to the synapses, immune response, lipid metabolism, and heme biosynthesis/degradation, even if the affected genes within these pathways were not identical. As the behavior phenotype is similar, this might suggest that different mechanisms mediate the anxiety phenotype in males and females.
The most significantly affected pathway in females was the circadian rhythm (Table 4) with differential expression of seven genes in response to EE 2 . Among them were per1a, bhlhe41, and nr1d1. Altered transcription of per1, bhlhe41, and nr1d1 has been observed in the pituitary of sub-adult female Coho salmon exposed to EE 2 , and circadian rhythm genes were among the top affected pathways also in that study (Harding et al., 2013). In zebrafish, regulation of circadian rhythm genes in whole brain differs between fish with high-and low-stress sensitivity (Rey et al., 2013). Although circadian rhythm genes are not FIGURE 4 | Brain mRNA expression of (A) msmo1, (B) sqlea, (C) lss, and (D) hmgcs1 in zebrafish males developmentally exposed to 0 or 3 ng/L EE 2 and remediated in clean water for 120 days. The mRNA expression was quantified with qPCR and normalized against the expression of 18S RNA and Elf1a, displaying the relative gene expression. Data represent mean ± 95% CI of 10 samples/group, *p < 0.05, **p < 0.01. * Significantly different from control group.
widely recognized as being estrogen responsive, several studies have found alterations in clock genes in response to estrogen exposure. In rodents, clock genes are altered in both the uterus and SCN in response to estrogen (Nakamura et al., 2001(Nakamura et al., , 2005(Nakamura et al., , 2008. Also, a connection has been found between the onset of locomotor activity and the different phases of the estrous cycle (see Krizo and Mintz, 2014 for review). Recent research in rodents shows that the circadian rhythm is highly responsive to gonadal steroids. Estrogens administered to ovariectomized rats have an effect on the period length and onset of locomotor activity (Krizo and Mintz, 2014). Photoperiod and melatonin, the major hormone in the circadian rhythm, is implied in the modulation of anxiety in rodents (Golombek et al., 1996;Nava and Carta, 2001;Trainor et al., 2007;Bilu and Kronfeld-Schor, 2013;Adamah-Biassi et al., 2014;Kumar et al., 2014;Nagy et al., 2014), possibly mediated by interaction with the GABAergic system (Golombek et al., 1996). Taken together, the observed effects on genes related to the circadian rhythm in female brain might represent a plausible connection to the stress sensitive behavior phenotype. No differential effects on clock genes were observed in male brains, but acute EE 2 exposure has been previously found to upregulate hepatic CRY2A expression in male zebrafish (Wit et al., 2010) and downregulate CLOCK in fathead minnow testis (Garcia-Reyero et al., 2009). This study is, to the best of our knowledge, the first presenting RNA sequencing data in EE 2 exposed male brains.
In males, the top putative pathway identified was cholesterol biosynthesis with four genes upregulated ( Table 5). Cholesterol biosynthesis and metabolism has been found to be affected by estrogen in several previous studies (Cypriani et al., 1988;Hoffmann et al., 2006;Sharpe et al., 2007;Flores-Valverde et al., 2010;Hogan et al., 2010). Zebrafish embryos exposed to low concentrations of EE 2 for 48 h (Schiller et al., 2013), as well as livers of adult zebrafish exposed to 30 ng/L EE 2 for 4 or 28 days (Wit et al., 2010) showed up-regulation of genes involved in steroid biosynthesis and metabolism Also, steroidogenic enzymes were altered in fish brains by exposure to estrogenic compounds (Arukwe, 2005;Lyssimachou and Arukwe, 2007). This study identified up-regulated expression of four genes directly involved in cholesterol biosynthesis, sqle, lss, msmo1, and hmgcs1. RNAseq of brains of fluoxetine-exposed zebrafish males has recently revealed the same genes to be down-regulated, concomitant with anxiolytic effects in NT (Wong et al., 2013). Also, in a follow-up to the fluoxetine study, Wong et al. compared brain transcriptomes of two zebrafish lines specifically bred for different stress-coping behaviors (Wong et al., 2015). The results showed that fish with high stress sensitivity also had upregulation of sqle, lss, msmo1, and hmgcs1 compared to fish with lower stress sensitivity (Wong et al., 2015). In rodents, dietary cholesterol affects anxiety behavior (Hu et al., 2014), and reduced brain cholesterol was observed following chronic mild stress (Sun et al., 2015). Locally synthesized cholesterol is the precursor for neuro-steroids implied as potential mediators of anxiety in mammals (Gunn et al., 2011;Zorumski et al., 2013), augmenting GABA A signaling (Gunn et al., 2011) but also interacting with NMDAR  and 5-HT1A receptor (Sun et al., 2015). If an association between cholesterol biosynthesis and anxiety in zebrafish males exist must await further confirmation, but it is intriguing to note a possible connection between GABA signaling and behavior in the top putative pathways in both males and females. Also, the association between cholesterol and neuro-steroids represent the only connection observed with the serotonergic system, via interaction with 5-HT1A.
In females, three genes with the ability to bind lipids were found to be down-regulated by developmental EE 2 exposure ( Table 4). One of these genes, apoeb, codes for the protein apolipoprotein Eb and is present in both brain and plasma. The human/mouse ortholog, apolipoprotein E, is produced largely by astrocytes and is believed to be one of the most important cholesterol transport proteins in the brain (Björkhem and Meaney, 2004). Apoeb mRNA levels has previously been found to be downregulated in adult zebrafish male, but upregulated in female, livers after exposure to 30 ng/L EE 2 for 4 days (Wit et al., 2010). Downregulation of abca12, coding for an ABCtransporter, was also observed in the brains of EE 2 -exposed female zebrafish. Members of the ABC-transporter superfamily may be responsible for excluding circulating cholesterol from the brain (Björkhem and Meaney, 2004). However, the function of abca12 remains largely unknown.
The current study identified six differentially expressed genes with recognized functions in the synapse. One of these genes, sv2b, encoding a synaptic vesicle protein (Bajjalieh et al., 1993) was upregulated in males developmentally treated with EE 2 while sv2b was downregulated. Interestingly, two sv2b genes were also shown to be affected by fluoxetine treatment in male brains (Wong et al., 2013). In resemblance to the current study (Table 5), one sv2b paralogue was upregulated, while the other was down-regulated in exposed males (Wong et al., 2013). Also, cholesterol is essential in the process of myelination of axons (Björkhem and Meaney, 2004) and genes involved in cholesterol biosynthesis appear to play an important role in brain development and plasticity (Mathews et al., 2014). Among them is hmgcs1, one of the genes mentioned above to be upregulated by EE 2 in male brains ( Table 5). Knock-down of the encoded enzyme caused oligodendrocytes to migrate past their target axons (Mathews et al., 2014). Furthermore, the down-regulation of arhgef9b expression in male brains a gene, encoding collybistin, which has a pivotal role in the formation of postsynaptic glycine and inhibitory GABA receptor clusters (Kalscheuer et al., 2009), represent an additional change that might affect inhibitory GABA signaling (Kalscheuer et al., 2009). In females, an upregulation of brain rtn4rl2a expression, a gene whose product in humans and rodents is associated with prevention of axonal outgrowth (McDonald et al., 2011) and a regulator of synaptic plasticity (Lee et al., 2008), and of bai1b expression, coding for a negative regulator of brain angiogenesis that has recently been associated with synaptogenesis (Duman et al., 2013) was observed. In addition, si:ch211-120g10.1, a gene where the human ortholog RNF39 encodes a protein suggested to play a part in synapse plasticity (Matsuo et al., 2001), showed decreased expression in exposed female brains. Taken together, EE 2 exposure during zebrafish development induces several persistent alterations in the brain that likely affect the outcome of nerve signaling in the adult fish.
In the female brain, a surprisingly high number of differentially expressed novel miRNA was observed ( Table 6). miRNAs cause modification of gene expression on posttranscriptional level, and are involved in sex differentiation in the brain (McCarthy and Nugent, 2015). Estrogens have been shown to regulate the expression and biogenesis of miRNA (Gupta et al., 2012;Klinge, 2012). It is possible that the differentially expressed miRNAs contributes to further effects at the protein level in the EE 2 -exposed female brain. As the function of the observed miRNAs are not known, no conclusions can be drawn and we can only report our findings awaiting future clarification of their function.

CONCLUSIONS
This study has verified our previous findings of persistent effects on stress behavior in response to developmental exposure to low doses of EE 2 . RNA sequencing of male and female brain transcriptome revealed differences in differential expression between the sexes. No effects on genes belonging to the stress axis was observed in any of the sexes, but the expression of several genes in the regulation of circadian rhythm in females and cholesterol biosynthesis in males were found to be affected. Both pathways have previously been implied in anxiety regulation. Also, altered expression of several genes associated with synaptic function was observed, which might turn out to be important for the developmental modulations resulting in an anxiety phenotype. Further studies are needed to evaluate the significance of these findings, representing an initial survey of the effects of developmental exposure to the ubiquitous environmental contaminant EE 2 on the brain transcriptome in the adult zebrafish.

AUTHOR CONTRIBUTIONS
TP planned, designed, and conducted experiment, performed qPCR, prepared figures and tables and completed the manuscript. KV planned, designed, and conducted experiment, performed manual functional classifications of coding genes and prepared tables, and wrote first draft of the paper and thereafter reviewed drafts of the paper. NR planned, designed, and conducted experiment and reviewed drafts of the paper. TK performed bioinformatics and biostatistics and reviewed drafts of the paper. PD performed statistical analyses and prepared figures and reviewed drafts of the paper. IP participated in and held main responsibility for the project from planning of experiment to the final manuscript.