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₂/L, respectively. EE₂ (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₂/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₂ 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).

Water samples were collected from different tanks at nine separate occasions during the exposure and stored in darkness at −20°C until analysis. EE₂ 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₂ content. The quantification range was 0.5–100 ng/L, using EE₂-d4 as internal standard.

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 pre-heated 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.

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).

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₂, while for females, brains were from unexposed females and females exposed to 10 ng/L EE₂. 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™ 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™ 2 system using the Ion Proton™ Template OT2 200 v3 Kit (Life Technologies) chemistry. Enriching was conducted using the Ion OneTouch™ ES (Life Technologies). Samples were loaded on an Ion PI v2 Chip (2 samples per chip) and sequenced on the Ion Proton™ System using Ion Proton™ Sequencing 200 v3 Kit (200 bp read length, Life Technologies) chemistry.

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.

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).

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® Reverse Transcription Kit (Qiagen) from 1 μg RNA per sample and then diluted 1:10 for qPCR reactions. Bio-Rad C1000 Touch™ Thermal cycler, CFX96™ 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™ Universal SYBR® 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).

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 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.

The actual concentration of EE₂ was determined from water samples from 9 occasions during the 80 days exposure period. The determined actual concentrations of EE₂ (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₂.

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₂ 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).

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₂ had significantly increased latency to upper half (z = 2.50, 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.

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 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₂ 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).

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₂ 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 with 10 ng/L EE₂ compared to control females. Males exposed to 3 and 10 ng/L EE₂ 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.

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₂-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₂. 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₂ exposure in both male and female brains: sv2b, coding for a synaptic vesicle protein, was upregulated by EE₂ (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₂ (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 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).

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₂ 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₂-exposed males than in control brains (Chisq = 1.86, df = 1, p = 0.173).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.