Differential Molecular Responses of Zebraﬁsh Larvae to Fluoxetine and Norﬂuoxetine

: The occurrence of psychopharmaceuticals in aquatic ecosystems is a growing problem. Fluoxetine (FL) and its metabolite norﬂuoxetine (NF) are selective serotonin reuptake inhibitors. Although they may be potentially harmful to non-target species, available knowledge on the effects of NF is sparse, relative to FL. This study aimed at contributing to the body of knowledge about the modes-of-action (MoA) of these compounds and their underlying mechanisms eliciting hazardous effects during the early development of the teleost model zebraﬁsh ( Danio rerio ). One hour post-fertilisation (hpf), embryos were exposed up to 80 hpf to these compounds at levels found in surface waters and higher (FL, 0.0015 and 0.05 µ M; NF, 0.00006 and 0.0014 µ M). Developmental anomalies were observed at 8, 32 and 80 hpf. Larvae were collected at 80 hpf to assess the expression of 34 genes related to FL and NF MoA and metabolism, using qPCR (quantitative reverse transcription PCR). Results showed that both compounds elicited an increased frequency of embryos exhibiting abnormal pigmentation, relative to controls. Gene expression alterations were more pronounced in FL- than in NF-exposed larvae. Cluster Analysis revealed two groups of genes discriminating between the drugs. for their marked opposing responses. Globally, downregulation of gene expression was typical of FL, whilst upregulation or no alteration was found for NF. These clusters identiﬁed were linked to the adrenergic pathway and to the retinoid and peroxisome proliferator-activated nuclear receptors. Overall, our data contradict the prevailing notion that NF is more toxic than FL and unveiled the expression levels of genes drd2b , 5-ht2c and abcc2 as possible markers of exposure to FL.


Introduction
The detection of psychopharmaceuticals in aquatic ecosystems is a recognised problem of concern among the public and the scientific community. This is related to the increasing concentrations detected in water and sediments as a consequence of their growing consumption, linked with the massive increase in depression occurring worldwide [1]. Among these psychopharmaceuticals are selective serotonin reuptake inhibitor (SSRI) antidepressants. Despite the increasing number of studies published in the last decade, the effects of SSRIs on non-target organisms are still not fully understood [2]. Fluoxetine (FL) is an SSRI antidepressant [3] prescribed to treat obsessive-compulsive disorders, moderate to severe within the range found in environmental samples, with bioaccumulation factors higher than those found for fluoxetine [30,31].
This study investigated and compared the effects of exposure to FL or NF on zebrafish embryotoxicity and gene expression levels of newly hatched larvae. This model was chosen due to the transparency of its eggs, which allows easy monitoring of embryonic development, as well as the availability of the full genome in public databases. Thirtyfour target genes coding for proteins with different functions in detoxification and several neurohormonal receptors involved in the MoA of both FL and NF were assessed. The outcome of this work adds knowledge about the MoA and toxic effects of both compounds in aquatic organisms, which are still not fully known, and identifies possible molecular biomarkers of exposure to these compounds in fish.

Zebrafish Rearing and Reproduction
Zebrafish breeders (Wildtype AB) were maintained in in-house certified facilities at CIIMAR-Interdisciplinary Centre of Marine and Environmental Research (Matosinhos, Portugal)-in 70 L aquaria with water circulation and continuous aeration, at a stocking density of 1.5 fish L −1 . The temperature was kept at 27 ± 1 • C and the photoperiod was 14 h light:10 h dark. The animals were fed with Tetramin XL Flakes two times a day, a total of 2.5 g. For reproduction, breeders were placed in a breeding aquarium (1 female:2 male) bearing a bottom net covered with glass marbles.

Zebrafish Embryo Toxicity Test
The assays were carried out as described in Cunha et al. [19]. Briefly, the tests were done in 24-well plates containing 10 embryos in each well. A total of 40 embryos were exposed per treatment in each assay. The embryos were exposed from 1 hpf (hour postfertilisation) until 80 hpf to 0.0015 and 0.05 µM of FL (FL1 and FL2, respectively) and 0.00006 and 0.0014 µM of NF (NF1 and NF2, respectively). Tested concentrations were selected to include one concentration in the range found in environmental samples (FL1 and NF1) and a higher one for each chemical. The FL stock solution was in DMSO and that of NF was prepared in water. The experimental design included two control groups: water and DMSO (DMSO 0.004%). The assays were performed in duplicate each time and repeated three times, each with a different batch of embryos. Observations were made at 8, 32 and 80 hpf using a Nikon Eclipse TS100 inverted microscope to check for embryo lethality and abnormal development by assessing and recording a list of endpoints defined according to the existing literature [32]. Table 1 presents the list of endpoints recorded at each developmental stage.
The day before each assay, the microplate wells were filled with the respective test solution, to prevent decreasing bioavailability of the compound during the assay due to adsorption to the testing plates. The test solutions were renewed daily to avoid microbial contamination. At the end of the assays, the larvae were pooled in groups of 40 per experimental condition and preserved in RNAlater and kept at −80 • C until the quantification of gene expression levels.

RNA Extraction and cDNA Synthesis
Total RNA was extracted with the Illustra RNAspin Mini RNA Isolation kit (GE Healthcare, Chicago, Illinois, USA), following the manufacturer instructions. A BioTek microplate spectrophotometer, equipped with Take3 (2 µL volume), was used to quantify the extracted Pericardium oedema X X Unhatched embryos X
The identities of the amplicons obtained were confirmed through cloning and sequencing of the fragments of DNA, as described previously [33] (Appendix A). Efficiency curves were determined using eight dilutions ranging from 0.05 to 50 ng/µL of cDNA. qRT-PCR was done in an Eppendorf Mastercycler realplex 4 (Eppendorf, Hamburg, Germany) using PerfeCTa SYBR ® Green SuperMix from QuantaBio (Beverly, ME, USA) and 20 µL reaction volume (Appendix B). Duplicated were run for all reactions, according to the following cycles: one cycle at 95 • C for 3 min followed by 40 cycles at 95 • C for 10 s; 51 • C (vmat2), 55 • C (5-ht2c and drdb1) or 54 • C (remaining genes) for 30 s; and 72 • C for 30 s. Blank samples and melting curves were run for each gene. Gene expression was quantified by normalisation to reference genes using the algorithm Normfinder [34], which indicated actb1 and rpl8 as the most suitable reference gene combination. The method of Pfaffl, which accounts for the efficiency of the primers, was used to calculate the relative expression [35].

Statistical Analysis
Differences among experimental conditions in cumulative mortality and embryotoxicity endpoints were analysed with the χ2 test. For the embryotoxicity endpoints, the data analysed were the total number of embryos/larvae showing malformations at 32 and 80 hpf. Data were plotted as means and respective standard errors (SE). Considering that no significant differences were found between the water and the solvent controls, for any of the endpoints or expression levels analysed, all statistical comparisons were made against the solvent control, which for simplicity was indicated as control (ctr). Cluster analysis was used to investigate groups of genes with similar behaviour across substances and test concentrations. The elements used for the cluster analysis were the replicate expression values (in fold change), and the analysis was based on Pearson correlation values. Multivariate Analysis of Variance (MANOVA; based on the Wilks' Lambda criterion) was used to investigate differences among treatments for the most relevant gene clusters (fold change) identified. ANOVA followed by the Tukey HSD was used to identify homogeneous subsets. Statistical analysis was performed on TIBCO ® Statistica software (14.0.0.15). The significance level was set at 5% for all statistical tests. Raw data are available in [36].

Zebrafish Embryo Toxicity Test
No statistically significant differences were found between the two control groups for any of the parameters analysed throughout the work. Cumulative mortality in the water control was 14%, while the solvent control had cumulative mortality of 16%. For all tested concentrations, mortality occurred between 8 and 32 hpf ( Figure 1). After that period, no mortality was observed. Differences in mortality in relation to ctr (X 2 (2,720) = 8.22, p = 0.016) were found at 8 hpf for NF1 and at 32 hpf for FL1 and FL2 (X 2 (2,720) = 15.0, p = 0.0006). For NF, a significant reduction in the mortality rate was found; for FL, mortality was significantly increased, compared to the control group. For gross malformations, significant differences were found (p < 0.05, compared to the control) for both compounds ( Figure 1). A significant increase in gross malformations (~3-fold) was found at 32 hpf for FL2 (X 2 (2,240) = 23.3, p < 0.00001), while for NF, a significant increase in gross malformations in relation to control was found at 80 hpf for both tested concentrations (2.5-fold for NF1 and 3-fold for NF2, X 2 (2,656) = 26.6, p < 0.00001). Significant differences in gross malformations were mainly due to the appearance of anomalies in pigmentation in FL-(~8% for FL2 vs. 0% for ctr) and NF-exposed larvae (~7% for NF1 and~9% for NF2 vs. 0% for ctr at 80 hpf) ( Figure 2); skeletal abnormalities were also observed in exposed embryos, though at a much lower frequency (<5% vs. 2% for ctr).
hpf. Data were plotted as means and respective standard errors (SE). Considering that no significant differences were found between the water and the solvent controls, for any of the endpoints or expression levels analysed, all statistical comparisons were made against the solvent control, which for simplicity was indicated as control (ctr). Cluster analysis was used to investigate groups of genes with similar behaviour across substances and test concentrations. The elements used for the cluster analysis were the replicate expression values (in fold change), and the analysis was based on Pearson correlation values. Multivariate Analysis of Variance (MANOVA; based on the Wilks' Lambda criterion) was used to investigate differences among treatments for the most relevant gene clusters (fold change) identified. ANOVA followed by the Tukey HSD was used to identify homogeneous subsets. Statistical analysis was performed on TIBCO ® Statistica software (14.0.0.15). The significance level was set at 5% for all statistical tests. Raw data are available in [36].

Zebrafish Embryo Toxicity Test
No statistically significant differences were found between the two control groups for any of the parameters analysed throughout the work. Cumulative mortality in the water control was 14%, while the solvent control had cumulative mortality of 16%. For all tested concentrations, mortality occurred between 8 and 32 hpf (Figure 1). After that period, no mortality was observed. Differences in mortality in relation to ctr (X 2 (2,720) = 8.22, p = 0.016) were found at 8 hpf for NF1 and at 32 hpf for FL1 and FL2 (X 2 (2,720) = 15.0, p = 0.0006). For NF, a significant reduction in the mortality rate was found; for FL, mortality was significantly increased, compared to the control group. For gross malformations, significant differences were found (p < 0.05, compared to the control) for both compounds ( Figure 1). A significant increase in gross malformations (~3-fold) was found at 32 hpf for FL2 (X 2 (2,240) = 23.3, p < 0.00001), while for NF, a significant increase in gross malformations in relation to control was found at 80 hpf for both tested concentrations (2.5-fold for NF1 and ~3-fold for NF2, X 2 (2,656) = 26.6, p < 0.00001). Significant differences in gross malformations were mainly due to the appearance of anomalies in pigmentation in FL-(~8% for FL2 vs. 0% for ctr) and NF-exposed larvae (~7% for NF1 and ~9% for NF2 vs. 0% for ctr at 80 hpf) ( Figure 2); skeletal abnormalities were also observed in exposed embryos, though at a much lower frequency (<5% vs. 2% for ctr). . Cumulative mortality and rates of gross malformations recorded in embryos and larvae of zebrafish exposed to 0.0015 and 0.05 μM fluoxetine (FL1 and FL2, respectively) (left graph) or 0.00006 and 0.0014 μM norfluoxetine (NF1 and NF2, respectively). Results are expressed as mean ± SE. * indicates statistically significant differences in relation to the solvent control (ctr) (p < 0.05). Figure 1. Cumulative mortality and rates of gross malformations recorded in embryos and larvae of zebrafish exposed to 0.0015 and 0.05 µM fluoxetine (FL1 and FL2, respectively) (left graph) or 0.00006 and 0.0014 µM norfluoxetine (NF1 and NF2, respectively). Results are expressed as mean ± SE. * indicates statistically significant differences in relation to the solvent control (ctr) (p < 0.05).

Gene Expression.
The Cluster Analysis identified eight groups of related genes (clusters a to h, Figure 3). Analysis of the heatmap for gene expression showed that clusters e and h were the most relevant to distinguish between NF and FL ( Figure 3). In contrast, cluster g was useful to distinguish the two FL concentrations evaluated. MANOVA confirmed that the treatment significantly influenced the expression of these groups of genes (cluster e: F (15,23) = 6.36, p < 0.0001; cluster g: F (  The Cluster Analysis identified eight groups of related genes (clusters a to h, Figure  3). Analysis of the heatmap for gene expression showed that clusters e and h were the most relevant to distinguish between NF and FL ( Figure 3). In contrast, cluster g was useful to distinguish the two FL concentrations evaluated. MANOVA confirmed that the treatment significantly influenced the expression of these groups of genes (cluster e: F(15,23) = 6.36, p < 0.0001; cluster g: F(12,24) = 2.63, p = 0.02; cluster h: F(27,12) = 3.95, p = 0.007).
Cluster e was composed of the genes drd2b, 5-ht1a, dat, mao and cyp3a65. Of these, dat, mao and 5-ht1a were the most indicative genes, in terms of expression patterns and differences among group treatments, as indicated by the ANOVA results (p < 0.05, Table S2). Three separated homogeneous subsets were found for the expression of dat, distinguishing the effects of NF and FL, and the two concentrations tested for NF (Figures 3 and 4). NF exposure caused upregulation of dat expression, which was more intense in the highest test concentration. FL exposure elicited downregulation of dat expression. The expression of mao was clearly different under the effect of NF (upregulated) and that of FL (downregulated) (Figures 3 and 4). In addition, 5-ht1a expression was significantly different in FL2, compared to the remaining groups, for which clear downregulation was noted. Additionally, in larvae exposed to FL, significant differences in expression were detected for all genes (p < 0.05, compared to the solvent control, Figure 4). The most relevant changes were the decreased expression observed for 5-ht1a (~2-fold), dat (3-fold) and mao (~2-fold) in the FL1 group, relative to the solvent control.
Cluster g comprised the genes net, serta, pparb and ahr2 (Figures 3 and 5). Globally, the two FL concentrations had opposing effects on the expression of these genes; downregulation was found for FL1 and upregulation was observed in larvae exposed to FL2 (Figures 3 and 5). No significant differences relative to the solvent control were found for these genes, in larvae exposed to NF. For FL, significant downregulation of pparb and serta (~1-fold for both genes) was found for FL1 treatment, while significant upregulation was found for ahr2 in the FL2 treatment (1-fold) ( Table S2, Supplementary Material).
Cluster h was composed of genes raraa, rxrga, gst, 5-ht2c, rxrab, rxraa, rarga, rarab, drd1b and abcc2. Of these, rarga, rarab and rxraa were the most indicative genes, characterising the expression patterns and differences among group treatments of the whole group. Three homogeneous subsets were found for the expression of rarga, completely separating the effects of NF and FL, and the two concentrations tested for FL (Figures 3 and 6). NF exposure caused upregulation of rarga expression, while FL2 elicited its downregulation at 80 hpf. For rarab, two homogeneous subsets were identified, separating the larvae exposed to NF from those exposed to FL2. NF exposure also caused upregulation of rarab expression, compared to the downregulation found in larvae exposed to FL2. For rxraa, two homogeneous subsets were found, also pointing to upregulation of its expression, particularly in larvae exposed to NF2, compared to downregulation in larvae exposed to FL2 (Figure 4). These larvae also exhibited decreased expression compared to solvent control (Table S2, Supplementary material). Compared to solvent control larvae, both FL concentrations elicited inhibition of rarab, rarga, 5-ht2c and abcc2 expression ( Figure 6, Table  S2).  Cluster e was composed of the genes drd2b, 5-ht1a, dat, mao and cyp3a65. Of these, dat, mao and 5-ht1a were the most indicative genes, in terms of expression patterns and differences among group treatments, as indicated by the ANOVA results (p < 0.05, Table S2). Three separated homogeneous subsets were found for the expression of dat, distinguishing the effects of NF and FL, and the two concentrations tested for NF (Figures 3 and 4). NF exposure caused upregulation of dat expression, which was more intense in the highest test concentration. FL exposure elicited downregulation of dat expression. The expression of mao was clearly different under the effect of NF (upregulated) and that of FL (downregulated) (Figures 3 and 4). In addition, 5-ht1a expression was significantly different in FL2, compared to the remaining groups, for which clear downregulation was noted. Additionally, in larvae exposed to FL, significant differences in expression were detected for all genes (p < 0.05, compared to the solvent control, Figure 4). The most relevant changes were the decreased expression observed for 5-ht1a (~2-fold), dat (3-fold) and mao (~2-fold) in the FL1 group, relative to the solvent control.    Cluster g comprised the genes net, serta, pparb and ahr2 (Figures 3 and 5). Globally, the two FL concentrations had opposing effects on the expression of these genes; downregulation was found for FL1 and upregulation was observed in larvae exposed to FL2 (Figures 3 and 5). No significant differences relative to the solvent control were found for these genes, in larvae exposed to NF. For FL, significant downregulation of pparb and serta (~1-fold for both genes) was found for FL1 treatment, while significant upregulation was found for ahr2 in the FL2 treatment (1-fold) ( Table S2, Supplementary Material).
Cluster h was composed of genes raraa, rxrga, gst, 5-ht2c, rxrab, rxraa, rarga, rarab, drd1b and abcc2. Of these, rarga, rarab and rxraa were the most indicative genes, characterising the expression patterns and differences among group treatments of the whole group. Three homogeneous subsets were found for the expression of rarga, completely separating the effects of NF and FL, and the two concentrations tested for FL (Figures 3 and 6). NF exposure caused upregulation of rarga expression, while FL2 elicited its downregulation at 80 hpf. For rarab, two homogeneous subsets were identified, separating the larvae exposed to NF from those exposed to FL2. NF exposure also caused upregulation of rarab expression, compared to the downregulation found in larvae exposed to FL2. For rxraa, two homogeneous subsets were found, also pointing to upregulation of its expression, particularly in larvae exposed to NF2, compared to downregulation in larvae exposed to FL2 (Figure 4). These larvae also exhibited decreased expression compared to solvent control (Table S2, Supplementary material). Compared to solvent control larvae, both FL concentrations elicited inhibition of rarab, rarga, 5-ht2c and abcc2 expression ( Figure 6, Table S2). Relative expression of genes from cluster e. The small letters indicate the homogeneous subsets indicated by the Tukey HSD for the comparison of norfluoxetine based on MANOVA on for fold change (NF1 and NF2) and fluoxetine treatments (FL1 and FL2) (p < 0.05). The stars indicate statistically significant differences (p < 0.05) relative to the solvent control (from the results of the ANOVA on normalised expression, Table S2, Supplementary material). The horizontal line (x = 1) indicates the baseline expression for the solvent control.

Discussion
Embryo assays are very relevant to investigate potential hazardous effects of chemicals on aquatic organisms because animals are commonly more sensitive in their early developmental stages than juveniles or adults [15,16,19,20,[27][28][29]. Furthermore, they allow monitoring effects with impact on population maintenance that would otherwise be

Discussion
Embryo assays are very relevant to investigate potential hazardous effects of chemicals on aquatic organisms because animals are commonly more sensitive in their early developmental stages than juveniles or adults [15,16,19,20,[27][28][29]. Furthermore, they allow monitoring effects with impact on population maintenance that would otherwise be impossible to analyse later in life (i.e., teratogenic effects) [19,20]. This study evaluated the effects on zebrafish embryos and larvae of an antidepressant (FL) and its active metabolite (NF). These compounds were selected since more information about their MoA in non-target organisms is urgently needed. They are frequently detected in environmental samples [22,37,38] and FL consumption is expected to increase given the predicted trends of incidence of psychiatric disorders. Whilst they already represented a very high burden of disease, the incidence of mental health problems is also increasing significantly with the COVID-19 pandemic [1]. The study of pharmaceutical metabolites, such as NF, is also very pertinent since they are often active, exhibiting pharmacological properties and can also be more dangerous to aquatic fish than their parent compounds [39]. However, until very recently [27][28][29], their investigation has been neglected in scientific studies, hampering the assessment of their risk to aquatic animals due to a lack of toxicity, as also highlighted by other authors for NF and other pharmaceutical metabolites [40].

Zebrafish Embryo Toxicity Test
The zebrafish embryotoxicity tests were performed in tightly controlled conditions, including bathing of the assay microplates with the test solutions for 24 h before the embryos were exposed, to prevent losses of chemicals during the assay due to adsorption to the exposure material, and daily renewal of the media. Moreover, the pharmaceutical compounds under investigation are considered stable. Fluoxetine was shown to be photolytically and hydrolytically stable in several aqueous solutions, including natural waters, for at least 30 days, with calculated half-life times higher than 100 days [41]. Even though the information available for norfluoxetine is scarce, a study determined a seven-day halflife for this drug [42]. Previous investigations also indicated a great agreement between nominal and exposure concentrations over 24 h for fluoxetine and three days for norfluoxetine, with recovery rates mostly above 85% (Table S3, Supplementary material) [29,[43][44][45]. Given these data and the experimental conditions in our work, the nominal and exposure concentrations are expected to be close.
The zebrafish embryotoxicity assays carried out fulfilled the validation criterium of Organisation for Economic Co-operation and Development (OECD) guidelines 212 and 236, which indicates as a requirement hatching rates ≥80% for control groups. Few implications of SSRI exposure on survival were found at the tested concentrations. The mortality rate was low at 8 hpf but increased at 32 hpf in the FL and NF exposed embryos. This could be related to the fact that the monoaminergic systems through which the test drugs exert their therapeutic action were not active at 8 hpf. The formation of these systems is known to start at 24 hpf in zebrafish embryos [14,46]. This is supported by the results of our preliminary assays, where similar mortality levels were obtained at 32 hpf when the exposures were started with 24 hpf embryos, compared to the present work. A recent study using zebrafish [16] reported that both FL and NF exposure caused a dose-dependent inhibition of locomotor activity in larvae, and increased the embryo malformation rates. High mortality rates for both compounds were also observed at the highest tested concentrations (89.9 µM for FL and 60 µM for NF) [16]. However, it is important to notice that the concentrations associated with high mortality rates in [16] are much higher than the ones tested in the present work.
The antidepressant drugs investigated in this study elicited an augmentation of embryonic malformations at 32 hpf for FL and 80 hpf for NF. This is of concern for NF, once that a significant increase in relation to solvent control was found for an environmentally relevant concentration (NF1) at 80 hpf. The main anomalies found for both compounds were in pigmentation. This may be due to the interaction of FL and NF with the em-bryonic adrenergic receptors of melanophores [47]. The melanophores are pigment cells that contain eumelanin. They confer black or dark-brown pigmentation to larvae through the tyrosine kinase signalling pathway. Because these cells originate from neural crest tissue, and somites are crucial for patterning neural crest migration, these results highlight the need for future investigation of possible somite defects in FL-and NF-treated embryos [48]. Pigmentation anomalies were previously detected in zebrafish embryos and larvae treated with venlafaxine [27]. An interesting result of the present work is also the time frame for the appearance of pigmentation anomalies, which was different for both compounds, i.e., 32 hpf for FL and 80 hpf for NF. Such an outcome in the developmental window could be related to FL and NF metabolism and/or to differences in potency and ratio of their enantiomers [49]. According to the literature [50], fish can accumulate norfluoxetine over time, reaching values equal or even higher than those of the parent compound [30]. Differences in bioaccumulation and/or MoA could explain the higher frequency of anomalies observed in later phases of development for norfluoxetine-exposed larvae. A recent study investigating the effects of fluoxetine and norfluoxetine (racemate and its S-enantiomer) in the visual motor response of zebrafish larvae also showed that fluoxetine caused impairments at lower concentrations than the metabolite (both when tested as racemate or S-enantiomer) [28]. Moreover, the authors reported additive effects of FL and NF on the embryonic visual motor response. The effects they observed were dependent on the concentrations of the tested equimolar mixtures and appear to be in agreement with the results of the investigations previously done in rats [26]. In contrast, the investigation of Chai et al. [29] indicated that in zebrafish larvae, racemic NF elicited more severe alterations in the cardiac structure, and cardiac fibrosis, than racemic FL. However, racemic FL induced a more severe immune cell infiltration in cardiac tissue than the racemic NF. These authors also showed that racemic FL is biotransformed into NF, with a stronger enrichment in S-norfluoxetine than R-norfluoxetine [29]. Overall, the results reported are complex and dependent on the effects measured, in addition to the test concentrations and interactive effects of FL and NF.

Gene Expression
Gene expression analysis revealed differing global patterns of response to the two SSRIs. Some variability among replicates was found, probably owing to slight differences in the genetic background as wildtype strain animals was used, but overall, the results were in line with previous literature reports on the effects of FL and NF on the expression of these genes in zebrafish larvae [18,19,27]. Mainly significant downregulation was observed in FL-exposed larvae, against a tendency for no change or upregulation relative to the solvent control in those exposed to NF. Nevertheless, for NF larvae, significant differences (downregulation) were only observed for Cu/Zn sod (~1-fold, compared to the solvent control). Obtained results also denoted that the variation in gene expression was globally more accentuated in FL larvae, compared to NF. Data from the literature [50] indicate that NF is more potent than its parent compound. Moreover, it was reported that NF is more toxic (up to 10 times) than FL to aquatic crustaceans and protists such as Tetrahymena thermophila, Spirostomum ambiguum and Thamnocephalus platyurus [51,52]. However, our results contradict those findings, at least for zebrafish in the early stages of development. The results of this study indicate that at similar exposure levels (NF2 and FL1), NF larvae experienced lesser alterations in gene expression and lower mortality than FL larvae. Comparing these two treatments, it is observable that both caused significant anomalies in the embryonic development of zebrafish. At the molecular level, however, for FL, more significant differences were detected in the expression of detoxification genes and of those involved in SSRI MoA than for NF.
Results from the cluster analysis identified eight groups of related genes. Of these, clusters e and h were the most relevant to distinguish NF and FL at the molecular level, while cluster g was useful to distinguish the two fluoxetine concentrations evaluated. From cluster e, the genes dat, mao and 5-ht1a, which were the most indicative genes of differences among group treatments, are involved in the mechanism of action of both compounds. The genes dat, mao and 5-ht1a are present in the monoaminergic systems where both SSRI compounds act. The action of NF and FL in these genes is, however, completely opposite to each other: downregulation for FL and upregulation for NF. The SSRI MoA occurs through inhibition of reuptake by the monoamine presynaptic receptors, mainly the serotonin receptors [3,10]. Considering this, we can hypothesise FL is exerting an expected effect on genes present in monoaminergic systems, causing their downregulation. In contrast, NF-exposed larvae showed a distinct response, with a predominance of upregulation or mild effects. This, once again, is in agreement with previous findings about the potency of norfluoxetine, compared to fluoxetine, in the visual motor response of zebrafish larvae [28].
Cluster g comprised the genes net, serta, pparb and ahr2. Globally, we observed downregulation for the lowest FL concentration and upregulation for the highest one. FL acts by inhibiting monoaminergic systems, which leads to an expected inhibition of net and serta transporters. In the present data, such expected inhibition occurred for serta in FL1-exposed larvae. The opposite expression levels elicited by exposure to FL1 and FL2 also point out a non-monotonic response of this cluster of genes. This kind of response was already described for non-target organisms exposed to SSRI compounds, supporting the above-mentioned assumption [18,53,54]. Furthermore, in a recent study conducted with human intestine cellular lines, the authors found serotonin could induce cyp1a1 expression via the ahr through a process involving its uptake into the cell by sert (serotonin transporter) [55]. Although no evidence of this was described for fish, once both serta and ahr have similar patterns of expression and were comprised in the same cluster, it is possible that FL1 could have caused some disruption in the serotonin-mediated CYP pathway of drug metabolism. However, in-depth follow-up studies would be needed to elucidate this question.
Cluster h was formed by the genes raraa, rxrga, gst, 5-ht2c, rxrab, rxraa, rarga, rarab, drd1b and abcc2. Of these, rarga, rarab and rxraa were the most indicative, characterising the expression patterns and differences among group treatments of the whole group. Those three genes belong to the families of retinoid x and retinoid acid nuclear receptors. Retinoid x receptors have the particularity of forming homodimers, as well as heterodimers with other nuclear receptors (e.g., retinoic acid receptors), to regulate various physiological processes [56]. Differences between the compounds could lead to distinct impact on pathways where those genes are involved, creating potential opposing responses, including expression levels.
Overall, not many studies about effects of NF in fish species are available in the literature, making comparisons difficult. Nevertheless, the results reported herein are generally consistent with previous findings, showing similar up-or downregulation patterns were elicited by these substances [18,19,27,57].

Molecular Biomarkers of Exposure to Psychopharmaceuticals in Fish
Molecular biomarkers can be defined as a group of biomarkers that can be discovered using standard genomics and proteomics methodologies, as is the case of qPCR [58]. Zebrafish is a well-known test model for the study of behavioural neuroscience, depression and several mental diseases [59]. During the past few years, research about new molecular biomarkers of exposure to psychopharmaceuticals has been increasing, mainly for antidepressants ( Table 2).  Mianserin, a tetracyclic antidepressant that has antihistaminic and hypnosedative effects in humans, was reported to exert endocrine disruption on zebrafish [60]. Treatment with this antidepressant induced estrogenicity biomarkers, such as vitellogenin 1 and zona pellucida proteins. These could be defined as biomarkers of exposure to mianserin. Exposure of zebrafish to fluoxetine and sertraline produced a specific expression profile for both antidepressants [61]. However, the authors also reported a set of eight genes that showed similar patterns of expression for both compounds and mentioned the expression pattern of the gene encoding for FK506 binding protein 5 as a possible biomarker of exposure to SSRI expression. Amitriptyline, a tetracyclic antidepressant, was reported to cause changes in the adrenocorticotropic hormone, oxidative stress enzymes, antioxidant defences, nitric oxide production and the total activity of nitric oxide synthase [62]. It is important to notice that, similar to our study, the authors reported U-shaped dose-response relationships of the studied endpoints. In a previous study, Wu and colleagues [63] performed a complete transcriptomic analysis of zebrafish to evaluate embryonic exposure effects of amitriptyline, fluoxetine and mianserin. For that, the authors performed RNAseq (high-throughput RNA sequencing) analysis, followed by qPCR for the most promising set of genes resulting from sequencing. The genes egr1 and egr4 (early growth response 1 and 4) were under-expressed for all tested antidepressants, making it a possible molecular biomarker of exposure. Later on, exposure of zebrafish to FL for 96 hpf indicated the upregulation of serotonin and GABA transporters (slc6a4a, slc6a4b and slc6a11) as of interest as biomarkers of exposure [58]. A transgenerational study investigating the effects of early developmental exposure to FL [17] reported decreased cortisol levels associated with reduced exploratory behaviour of adult males for three generations. The authors concluded that FL exposure of 3 hpf embryos for 6 d caused disruption of the stress axis as manifested by a reduction of the basal and stress-induced cortisol levels, following an acute net handling stressor. The authors also found altered expression levels of genes controlling cortisol synthesis and the pathways of cholesterol or steroid synthesis. These FL-induced low cortisol levels were more prominent in males and were associated with significantly reduced exploratory behaviours for two generations. Comparative transcriptomic analysis of zebrafish larvae exposed to fluoxetine or paroxetine identified upregulation of genes involved in mitochondrial and neurodevelopmental processes as also of interest as markers of exposure [64]. Besides the identified malformations in cardiac development, Chai and colleagues [29], also detected abnormal expression (mostly upregulation) of five genes (gata4, tbx5, tgf-β, ctgf and cyp2k6) in 9 d larvae exposed to racemic FL or NF and S-norfluoxetine, which may provide useful biomarkers of exposure. In a transcriptomic study with sea bass juveniles exposed chronically (21 days) to antidepressants FL and venlafaxine, Costa et al. [65] also identified the serotonin receptor 5-ht3a as a promising molecular biomarker of exposure to the studied compounds, in addition to other genes involved in metabolism and neurotransmission (Table 2).
Our study analysed a vast set of genes using qPCR, unveiling drd2b, 5-ht2c and abcc2 as markers of exposure to FL. Furthermore, different profiles of expression for low and high FL concentrations were identified, which could be used as molecular biomarkers of exposure in zebrafish. On the other hand, although the results obtained in this study show differences in the gene expression patterns between FL and NF, no clear differences could be depicted between NF-exposed larvae and controls. This makes it difficult to determine a possible molecular profile that could be used as biomarker of exposure to the metabolite. Future research should focus on other types of methodologies and endpoints (e.g., biochemical analysis, transcriptomics and proteomic studies) to discover other promising molecular biomarkers of exposure to these compounds.

Conclusions
In conclusion, FL caused higher mortality levels than NF. In contrast, both compounds increased the frequency of abnormalities during zebrafish embryonic development; pigmentation anomalies were the most frequently found. At the molecular level, two clusters of genes with altered expression were useful to distinguish the two SSRIs. One cluster was linked to the adrenergic pathway eliciting pigmentation anomalies. The other was associated mostly with nuclear receptors. Interestingly, NF and FL were mostly found to elicit opposing genomic expression in early zebrafish larvae, possibly originating from differences in their intrinsic pharmacokinetic properties that could lead to different MoA. Overall, embryo pigmentation and the expression of several genes appear as promising biomarkers of fish exposure to antidepressant drugs. Future research should focus on conducting further in-depth investigations comparing the toxicity of parent psychopharmaceuticals and their metabolites, single and in mixtures as detected in wastewater effluents, to improve monitoring and risk assessment routines.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/w14030417/s1, Table S1: Genebank accession numbers, name, function and detailed primer information for the studied target and reference genes; Table S2: Results of the analysis of variance (ANOVA) performed for the normalized expression of each of the genes of the clusters with significant results in the MANOVA (cluster e, g and h) to investigating differences of treatments relative to the control solvent; Table S3: Nominal and exposure concentrations, or recovery%, reported for fluoxetine and norfluoxetine in previous works. Funding: This study was funded by the EU and FCT/UEFISCDI/FORMAS, through projects REWA-TER (ERA-NET Cofund WaterWorks2015, Water JPI) and by the project Ocean3R (NORTE-01-0145-FEDER-000064), supported by the North Portugal Regional Operational Programme (NORTE2020), under the PORTUGAL 2020 Partnership Agreement and through the European Regional Development Fund (ERDF). G.L. was supported by Strategic Funding UIDB/04423/2020 and UIDP/04423/2020 (FCT, ERDF). R.P. was supported by a fellowship from FCT (SRFH/BD/134518/2017).

Informed Consent Statement: Not applicable.
Data Availability Statement: Raw data supporting this study results can be found in [36].

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

Appendix A. Identification of Amplicons Procedure
To evaluate if the identity of the sequences was correct, qPCR (polymerase chain reaction) was performed in a Biometra thermocycler (T Gradient, Biometra GmbH, Goettingen, Germany). qPCR reactions were performed with the following components, volumes, and concentrations, in a 20 µL volume per reaction: 2 µL MgCl 2 (2.5 mM), 4 µL of 5× buffer, 1 µL of forward primer and 1 µL of reverse primer (1 µM), 0.4 µL of DNTP's (0.2 mM), 0.1 µL of TaqPolimerase (Promega) (2.5U), 9.5 µL water, and 2 µL of cDNA template (obtained from polls of~40 zebrafish embryos at 80 hpf stage). The reaction was performed with the following protocol: 2 min of denaturation at 94 • C; 40 cycles with a denaturation period of 30 s, followed by 30 s of annealing or hybridization at 51 • C, 54 • C, 55 • C (54 • C for all genes except for vmat (51 • C), 5-ht2c and drd1b (55 • C)) and a final 30 s period of polymerisation at 72 • C; and 10 min at 72 • C for a further elongation. Band size was assessed on a 2% agarose gel buffered with 1µLof TAE1x Gel Red, under UV light. Bands with the amplicon expected size were cut from and purified using the illustra GFX PCR DNA and Tm Gel Band Purification Kit (GE Healthcare). The fragments were inserted into the vector pGEM (pGEM(R) -T Easy Vector Systems -Promega) and introduced into E. coli with NovaBlue Singles Competent Cells (Novagen). Selected colonies were developed on solid medium for 10 h (35 g/L of LB Broth, ampicillin 0.1 mg/mL, X-gal 100 mM and IPTG 0.1 mM) at 37 • C. Plasmids were isolated from 5 mL of culture medium with 5 µL ampicillin and incubated overnight at 37 • C, with continuous stirring. DNA extraction was performed according to the Wizard Kit Plus SV Minipreps DNA Purification System (Promega), manufacturer instructions. Final products were sequenced by Stabvida (Portugal) and the sequence's identity was verified with the Alignment Basic Local Search Tool (Blast, blastn) at the National Centre for Biotechnology Information (NCBI).