The SR4G transgenic zebrafish line was kindly provided by Dr. Karl Clark (Mayo Clinic, Minnesota, USA) and Dr. Xiao-Yan Wen (University of Toronto, Ontario, Canada) (38). One-day post-fertilization embryos were transferred to Ottawa and placed in the quarantine room of the University of Ottawa aquatic facility per regulations of Animal Care and Veterinary Services (ACVS). The embryos were allowed to hatch and were raised in the quarantine facility until they reached six months of age. These zebrafish were used to establish our founder generation. The F₀ generation was derived from this founder line.

To study the life-long effects of developmental exposure to FLX on the central nervous system, we collected the larval heads, and brain tissues (Tel and Hyp) of male adult zebrafish. Two independent variables were studied: stress condition and FLX exposure ( Supplementary Figure 1 ). A total of 12 treatment and control groups were studied. For the larval study, all experiments were conducted with six replicates (experimental units) which each consisted of 22-28 larval heads (sampling unit) from individuals that were raised in separate containers to avoid pseudo-replication ( Supplementary Figure 1 ). Similarly, for adults, all experiments were conducted with six replicates (experimental units) prepared by pooling of 5 tissue samples (sampling unit: either Tel or Hyp) harvested from adult male zebrafish. The same adult was used to harvest the Tel and Hyp so that the pooled samples were from the same five adult zebrafish for Tel and Hyp in each replicate. The five individuals pooled together were raised in separate containers from other groups of 5 to avoid pseudo-replication ( Supplementary Figure 1 ). The details about FLX exposure and stress routine is provided in the following sections.

The stock of 1 mg/ml of fluoxetine (FLX; CAS# 56296-78-7, Millipore-Sigma, Burlington, MA, USA) was prepared in 99% ethanol as the vehicle. To confirm concentrations in the working solution (54 µg/L), fluoxetine was quantified by comparing its LC/MS peak area with a serial dilution of the fluoxetine standard compound. The serial dilution was prepared in 0.1% formic acid, 20% acetonitrile in water (LC: Agilent 1100 capillary HPLC, Agilent, Santa Clara, CA, USA; MS: Thermo LTQ Linear Ion Trap mass spectrometer, Thermofisher, CA, USA). The HPLC gradient was set to 15% B (0.1% formic acid in acetonitrile) to 45% B (0.1% Formic acid in acetonitrile) over 15 mins. The mass spectrometer was set on MS1 scan mode between 150-1500 m/z range. Chromatograms were plotted by monoisotopic ion at 310.2 m/z. Embryos from on the day of fertilization (0 dpf) were exposed to freshly prepared 54 µg/L FLX daily until they reached 6 dpf. The control groups were exposed to the vehicle only (final concentration of 0.005% ethanol) from 0-6 dpf. The total water volume of the rearing tanks was adjusted to 1 ml per larva. The treatment and control larvae were further divided into unstressed and stressed groups and were housed in separate shelves in the incubator (PR505750L, Precision™ plant growth chamber, 504L; Thermofisher) to avoid unwanted stressors. The unhatched or dead embryos/larvae were removed once daily. After Day 6, control and FLX-exposed lineages were raised in clean water until they reached six months. At that time, the sexes were segregated, and subsequent experiments performed on male adult zebrafish. Both the control and FLX-lineages were further divided into unstressed and stressed groups.

For both FLX- treated larvae and adults and their corresponding controls, the unstressed groups were not exposed to any handling or mechanical stress and were sacrificed immediately on the start of light cycle on the day of an experiment using ice-cold water. The stressed groups were exposed to a modified handling net stress described previously (14). Briefly, the larvae from each group were exposed to air for 1 minute trapped in a mesh net, afterwards returned to the embryo medium for 3 minutes for rest, followed by another 1-minute exposure to air trapped in a mesh net ( Supplementary Figure 2 ). The adults from each group were exposed to air for 3 minutes trapped in a mesh net, afterwards returned to the water for 3 minutes for rest, followed by another 3-minute exposure to air trapped in a mesh net. Sample collection was carried out 60 minutes after the last net stress based on our previous study (40). For 7 dpf larvae, a pooled sample of heads was acquired by removing the trunk and tail sections. This was done to reduce the number of abundant muscular and gastrointestinal-related genes in our data set. Each larval sample was a pool of 22-28 heads. The tissue samples were immediately frozen on dry ice and stored at -80 °C until further assessments. For adults, after careful dissection of the skull, the Tel and the Hyp were separately dissected and immediately placed in sample tubes and frozen on dry ice. Brain tissues from 5 adult fish were pooled and stored at -80°C until further experimentations. All pooled tissue samples were prepared in 6 replicates across all treatments (FLX-L: fluoxetine exposed larvae; FLX-Tel: Tel from adults exposed to FLX as larvae; FLX-Hyp: Hyp from adults exposed to FLX as larvae; cL: control larvae, cTel: control adult Tel; cHyp: control adult Hyp) ( Supplementary Figure 1 ).

Total RNA extraction from all samples was carried out using Qiagen RNeasy mini kits (Qiagen, Hilden, Germany) based on the manufacturer recommendations. The quality of the total extracted RNA was assessed using an Agilent 4150 TapeStation system (Agilent, Santa Clara, CA, USA). The RNA integrity number (RIN) was acquired for all samples and only those with RIN >8.5 were considered of sufficiently high quality and used for RNA sequencing (41).

Next-generation high-throughput RNA sequencing was carried out using Illumina technology (Illumina, San Diego, CA, USA). A total of 72 samples were sequenced, considering three tissues (L, Tel and Hyp), two independent factors (stressed vs. unstressed) and (FLX treatment vs vehicle control), and six experimental replicates for each studies condition. Standard single-read mRNA libraries were built using the TruSeq Stranded mRNA Library Prep kit (Cat. # 20020595, Illumina) with 150 ng of total RNA extracts. The poly-A enrichment technique was used to select mRNA. Thereafter, all mRNAs were fragmented, converted to cDNA using random hexamers, and indexed using TruSeq RNA CD Idx (Cat. # 20019792, Illumina) for downstream multiplex analysis. The prepared libraries were assessed using Agilent 4150 TapeStation system (Cat. # G2992AA, Agilent, Santa Clara, CA, USA). The approximate band size was 260 bp for each prepared library. The indexed libraries were normalized to 10nM and pooled in equal volumes. Sequencing was done using NSQ 500/550 Hi Output KT v2.5 (Cat. # 20024906, Illumina) kits with 75 cycle cartridges. Each pooled library was run on two separate flow cells. Read alignment, as well as gene mapping to the zebrafish reference genome (GRCz11; GCA_000002035.4) was performed using STAR (42). An average of 18 million reads was mapped across different samples (range considering all samples: 12-24 million reads).

For quality assessment, QA/QC plots were generated from BAM files using Qorts (43). Furthermore, raw counts were transformed to logarithmic pseudo counts, and the boxplot functions were used to visualize between sample distribution of pseudo counts for outlier identification. Boxplots showed similar distribution of pseudo counts between samples did not reveal outlier samples

Differentially expressed genes (DEGs) analysis was performed using the DESeq2 package (44). Given our study design had FLX exposure and stress status as factors, the DESeq function was used for comparing stressed and unstressed groups considering the FLX exposure. This function also reported an adjusted p-value (padj) based on the Benjamini-Hochberg procedure, which controls the false discovery rate (FDR). The threshold to identify DEGs was an adjusted p value ≤ 0.05 and a fold change (FC) ≥1.2, which has been reported as the optimal FC for nervous system transcriptomic analysis (45). A hierarchal clustering ( Figure 2 ) and a principle component analysis (PCA) ( Supplementary Figure 3 ) was performed on the differentially expressed set of genes to further evaluate and compare the transcription patterns between treatment conditions.

The lists of DEGs were analyzed using Qiagen Ingenuity Pathway Analysis (IPA; Qiagen, Hilden, Germany) to identify enrichment of canonical pathways among the differentially expressed genes. Since IPA does not have zebrafish genome annotation currently, the gene Ensembl IDs were converted to human orthologs using ID converter on DAVID (The Database for Annotation, Visualization and Integrated Discovery; https://david.ncifcrf.gov/tools.jsp). The duplicated paralogue genes in the zebrafish genome were assigned similar Homo sapiens orthologue Ensemble ID. All annotations reported here for pathways, networks, transcription regulators, and individual genes are according to Homo sapiens orthologue annotations.

A total of 27,162 (85%) gene IDs were successfully mapped by IPA. To avoid losing sensitive information, no attempt to delete or adjust duplicated probes was made. Therefore, the IPA function to assign a Z-score for activation or inhibition of any canonical pathway was not used. All analysis reported here is solely based on gene enrichment in canonical pathways. The IPA network prediction tool was used to cluster affected genes in functional molecular networks. This functionality was also used to determine potential upstream transcription regulators driving the observed alterations in the DEGs.

Our transcriptomic analysis showed that a total of 4456, 3172, and 5778 genes were significantly affected by developmental exposure to FLX in larvae (L), adult telencephalon (Tel), and adult hypothalamus (Hyp), respectively, regardless of stress condition. The total number of up-and down-regulated genes following developmental exposure to FLX in unstressed condition was 941, 1927, and 5055 for L, Tel, and Hyp, respectively ( Figure 1 ). In the unstressed state, 68% of DEGs were downregulated and 32% were upregulated in the FLX exposed larvae (unstr-FLX-L). Interestingly, this ratio was reversed in the same unstressed condition in Hyp (unstr-FLX-Hyp), with 63% and 37% of genes upregulated and downregulated, respectively. For adult Tel in the unstressed condition (unstr-FLX-Tel), 57% and 43% of genes were upregulated and downregulated, respectively. The total number of up-and down-regulated genes following developmental exposure to FLX in the stressed condition was 3515, 1245, and 723 for L, Tel, and Hyp, respectively ( Figure 1 ). In the stressed state, the ratio of up-and down-regulated genes was similar across all samples. The stressed larvae group (str-FLX-L) developmental exposure to FLX resulted in downregulation of 60% of the affected genes; the remaining 40% were upregulated. In the Hyp of stressed adult zebrafish (str-FLX-Hyp), 44% and 56% of genes were upregulated and downregulated, respectively, following developmental exposure to FLX. Moreover, 43% and 57% of genes were upregulated and downregulated in the str-FLX-Tel group, respectively ( Figure 1 ). Global transcriptional responses for each condition relative to matched controls are shown Supplementary Figure 4 .

Hierarchical clustering revealed the major effects of treatment (FLX vs control) and stress status (unstressed vs stressed) on gene expression patterns ( Figure 2 ). For larvae and adult Hyp, treatment (FLX vs control) was a primary driver of clustering, with subclustering based on stress status. For Tel, the stress status (unstressed vs stress) divided the main clusters ( Figure 2 ). The replicates within stressed and unstressed groups were highly correlated (tightly clustered) in FLX exposed groups across all tissue samples. In contrast, three samples from the unstressed larval control groups had highly similar expression patterns that were quite distinct from the rest of the controls in this group and were on a branch of their own.

Enrichment analysis in IPA was performed to identify the pathways affected by developmental exposure. The canonical pathways were ranked based on their p-value (p ≤0.05) within each condition. We ranked the most significantly affected pathways within each tissue sample. The top 10 pathways for each of the six main groups are shown for unstressed ( Figure 3A ) and stressed conditions ( Figure 3B ). Our data showed that the cholesterol biosynthesis was significantly dysregulated in the unstressed larvae by developmental exposure to FLX; with eight of the ten mostly affected pathways were associated with cholesterol metabolism ( Figure 3A ). Also, several pathways associated with neural cell signaling, neurotransmitter signaling, and regulation of synaptic junctions populated the highly affected pathways in both Hyp and Tel of the unstressed adults ( Figure 3A ).

In the stressed groups, DNA metabolism (4 out of 10 pathways) and neurotransmitter signaling (3 out of 10 pathways) were most significantly affected by developmental exposure to FLX in the stressed larvae ( Figure 3B ). In the Hyp and Tel of stressed adults’ pathways associated with neural cell-signaling, circadian rhythm, and energy metabolism were populating the mostly affected pathways; for example, Circadian Rhythm Signaling, Sirtuin Signaling Pathway, and Aldosterone Signaling, ( Figure 3B ).

We studied the life-long effects of FLX by tracing pathways that were significantly affected in all tissue samples (L, Hyp, and Tel) or at least in one of the adult samples (Hyp or Tel). This was performed separately for unstressed ( Supplementary Tables 1 , 3 ) and stressed ( Supplementary Tables 2 , 4 ) conditions. For example, fifteen pathways were affected in all studied samples in the unstressed conditions: including glucocorticoid receptor signaling pathway, Aldosterone signaling in epithelial cells, IL-22 signaling, ILK, signaling, and Myc mediated apoptosis signaling. ( Supplementary Table 1 ). Moreover, 47 pathways showed dysregulation in both Tel and Hyp of the unstressed adults which were exposed to FLX during their embryonic development including GNRH signaling, CREB signaling in neurons, Axonal guidance signaling, NF_κB signaling, STAT3 pathway, and synaptogenesis signaling pathway ( Supplementary Table 1 ). Also, 14 more pathways were only affected in the Hyp tissue of the unstressed adults with another 160 canonical pathways dysregulated in the Tel tissue of the unstressed adults after developmental exposure to FLX ( Supplementary Table 3 ). A similar result was observed in the stressed groups however, to a lesser extent. We observed 8 pathways dysregulated in the adult stress groups in Hyp and Tel; including Circadian Rhythm signaling, Sirturin signaling pathway, BAG2 signaling pathway, and Unfolded protein response ( Supplementary Table 2 ). Also, 32 more pathways were only affected in the Hyp tissue of the stressed adults with another 26 canonical pathways dysregulated in the Tel tissue of the stressed adults after developmental exposure to FLX ( Supplementary Table 4 ). Careful analysis showed that the top ten pathways within each group ( Figure 3 ) were amongst pathways showing life-long dysregulation upon developmental exposure to FLX ( Supplementary Tables 1–4 ).

Our findings showed that there were more significant affected pathways overall in the unstressed groups relative to the stressed groups in both FLX-exposed larvae and FLX-exposed adults ( Figure 3C and Supplementary Tables 1–4 ). Moreover, amongst the 40 pathways that were affected in both unstressed and stressed conditions more pathways were significantly affected (based on p-value) in the unstressed treatment groups compared to the stressed groups ( Supplementary Table 8 ).

An upstream regulator analysis was conducted in IPA to identify the potential key transcriptional regulators that drive the gene dysregulation observed ( Table 1 ). We found enrichment of a higher number of transcriptional regulators by FLX in zebrafish larvae in unstressed condition than the other groups. It was noted that CLOCK, a key transcription factor in regulating circadian rhythm, was affected in both larvae and adults and was more affected in stressed condition (unstr-FLX-L, p=0.02; str-FLX-L, p=0.04; str-FLX-Hyp, p=5^E-06; str-FLX-Tel, p=9^E-06). Such significant effects were not observed in the Hyp and the Tel in the unstressed condition (unstr-FLX-Hyp, p>0.05; unstr-FLX-Tel, p>0.05). Other transcriptional regulators that were observed in the adult Hyp and Tel in stressed fish include: NR3C1 (str-FLX-L, p=0.03; str-FLX-Tel, p=0.02), MAT1A (str-FLX-L, p=0.002; str-FLX-Tel, p=0.02), NPC1 (str-FLX-Hyp, p=0.02; str-FLX-Tel, p=0.04), SIRT1 (str-FLX-Hyp, p=0.02; str-FLX-Tel, p=0.02), and NCOA1 (str-FLX-Hyp, p=6^E-04; str-FLX-Tel, p=0.0001). Moreover, in the unstressed condition we observed several transcriptions regulators were affected in adult Hyp and Tel; for example, CPT1C (unstr-FLX-Hyp, p=0.03; unstr-FLX-Tel, p=0.002), SREBF1 (unstr-FLX-L, p=0.02; unstr-FLX-Tel, p=8^E-04), and SREBF2 (unstr-FLX-L, p=0.001; unstr-FLX-Tel, p=2^E-04), SCAP (unstr-FLX-L, p=0.01; unstr-FLX-Tel, p=0.002), and PPARGC1A (unstr-FLX-L, p=0.009; unstr-FLX-Tel, p=0.004).

The top 5 networks in each tissue sample and their IPA functional categories are shown ( Table 2 ). The scoring is based on the number of genes involved and the internal scoring algorithm of IPA. Our findings suggest that several predicted networks were affected in adults upon developmental exposure to FLX. For example, genes involved in “lipid metabolism” and “molecular transport” were significantly enriched in larvae and adult in both stressed and unstressed condition. Moreover, networks involved in “nervous system development” and “endocrine system development” were affected in the adult Hyp and Tel.

Many genes were dysregulated in both larvae and adults. We report that 46 and 117 genes were dysregulated (up- or down-regulated) upon developmental exposure to FLX in all tissue samples (L, Hyp, Tel) in the unstressed and stressed conditions, respectively ( Figure 4 ). The associated functional enrichment of these genes based on Gene Ontology annotations of Homo sapiens orthologues are shown in Supplementary Tables 5 , 6 for the unstressed and stress conditions, respectively. Downregulation predominated in the unstressed larvae (unstr-FLX-L, 38 out of 46 genes) and upregulation predominated in the unstressed adults (unstr-FLX-Hyp, 34 out of 46; unstr-FLX-Tel, 34 out of 46). In contrast, upregulation was more prevalent in the stressed larvae (str-FLX-L, 95 out of 117 genes) and downregulation was the main feature in the stressed adults (str-FLX-Hyp, 91 out of 117; str-FLX-Tel, 91out of 117). Interestingly, for most of the genes in the unstressed and stressed responses, the direction of the gene expression change in larvae was opposite to that in adults; 93% of genes and 82% of genes for unstressed and stressed conditions, respectively.

To evaluate the potential life-long impact of fluoxetine exposure on several pathways we clustered the pathways which were affected in both larval stage and adult stage (from either Hyp or Tel or both). Any given pathway that was significantly dysregulated compared to the control samples, given the description above, in either stressed or unstressed conditions was considered to be a potential within-generation affected pathway. A list of 512 genes associated with the major depressive disorder was extracted from the DisGeNET database (46). Our transcriptome analysis shows that more than 80% of the orthologues of these genes were affected by developmental exposure to FLX at least in one of the studied tissue samples ( Supplementary Table 6 ). The main dysregulation was in adult Hyp in the unstressed condition, with 271 depression-associated genes significantly affected. Among these, 73% were upregulated upon developmental exposure to FLX. Examples of genes downregulated in adult brain following developmental exposure to FLX in the unstressed condition include: nrgn (unstr-FLX-Hyp: FC= -1.4, p=0.04; unstr-FLX-Tel: FC= -1.3, p=0.008), fabp2 (unstr-FLX-Hyp: FC=-1.6, p=0.02; unstr-FLX-Tel: FC=-2, p=0.04), klf11a (unstr-FLX-Hyp: FC=-1.4, p=0.007; unstr-FLX-Tel: FC= -1.45, p= 0.005), and pvalb6 (unstr-FLX-Hyp: FC=-1.2, p=0.01; unstr-FLX-Tel: FC= -1.6, p= 0.03). Similar downregulation was also observed in the stressed adult zebrafish; for example, abcc1 (str-FLX-Hyp: FC= -1.7, p= 1.5^E0-8; str-FLX-Tel: FC= -1.5, p= 1.4^E0-12), fkbp4 (str-FLX-Hyp: FC= -1.4, p=1.2^E0-7; str-FLX-Tel: FC=-1.3, p=9^E-016), pclob (str-FLX-Hyp: FC= -1.4, p= 0.04; str-FLX-Tel: FC=-1.2, p= 0.02), srrt (str-FLX-Hyp: FC= -1.4, p= 0.00002; str-FLX-Tel: FC= -1.3, p= 0.002), and tspo (str-FLX-Hyp: FC= -2, p= 7^E0-7; str-FLX-Tel: FC=-1.7, p= 2^E0-9).

We also found that several genes associated with development of depressive disorders were upregulated in adult brain following developmental exposure to FLX. For example, elovl5 (unstr-FLX-Hyp: FC= 1.7, p=0.001; unstr-FLX-Tel: FC= 1.9, p= 0.04; str-FLX-Hyp: FC=1.9, p= 9^E0-6; str-FLX-Tel: FC=1.9, p=0.008), tefb (unstr-FLX-Hyp: FC=2, p= 4^E0-6; unstr-FLX-Tel: FC= 1.5, p= 0.002; str-FLX-Hyp: FC= 1.6, p= 0.002; str-FLX-Tel: FC= 1.8, p= 7^E0-10) and klf15 (unstr-FLX-Hyp: FC= 2.5, p= 2^E0-7; unstr-FLX-Tel: FC= 1.7, p= 0.002; str-FLX-Hyp: FC= 1.6, p= 0.006; str-FLX-Tel: FC= 1.4, p=0.03) were upregulated in both adult Tel and Hyp in both unstressed and stressed conditions. Several genes associated with development of depressive disorders were upregulated in adult brain following developmental exposure to FLX only in the unstressed condition; e.g., nr3c2 (unstr-FLX-Hyp: FC= 3.12, p= 2^E0-5; unstr-FLX-Tel: FC= 1.6, p= 0.01), clocka (unstr-FLX-Hyp: FC= 2, p= 2^E0-6; unstr-FLX-Tel: FC= 1.4, p= 0.03), grm3 (unstr-FLX-Hyp: FC= 2.3, p= 5^E0-6; unstr-FLX-Tel: FC= 1.8, p= 0.02), and ntrk2a (unstr-FLX-Hyp: FC= 1.7, p= 3^E0-12; unstr-FLX-Tel: FC= 1.4, p= 0.006). Additionally, genes associated with development of depressive disorders were upregulated in adult brain following developmental exposure to FLX only in stressed condition, e.g., grm2a (str-FLX-Hyp: FC=1.2, p= 0.01; str-FLX-Tel: FC= 2, p= 1.2^E0-7), per1b (str-FLX-Hyp: FC= 1.8, p= 5^E0-4; str-FLX-Tel: FC= 1.7, p= 2^E0-8) and per1a (str-FLX-Hyp: FC= 3, p= 5^E0-6; str-FLX-Tel: FC= 2.6, p= 2^E0-7) ( Supplementary Table 7 ).

Finally, a list of genes with implication is epigenetics was extracted from the Gene Ontology database using Homo sapiens GO: biological process annotations. Table 3 shows the genes that were significantly dysregulated following developmental exposure to FLX in larvae and adult zebrafish in both unstressed and stressed conditions. Upregulation of these genes in unstressed condition in =Hyp and Tel of adult zebrafish was the main effect.