Jnk1 and downstream signalling hubs regulate anxiety-like behaviours in a zebrafish larvae phenotypic screen

Current treatments for anxiety and depression show limited efficacy in many patients, indicating the need for further research into the underlying mechanisms. JNK1 has been shown to regulate anxiety- and depressive-like behaviours in mice, however the effectors downstream of JNK1 are not known. Here we compare the phosphoproteomes from wild-type and Jnk1-/- mouse brains and identify JNK1-regulated signalling hubs. We next employ a zebrafish (Danio rerio) larvae behavioural assay to identify an antidepressant- and anxiolytic-like (AA) phenotype based on 2759 measured stereotypic responses to clinically proven antidepressant and anxiolytic (AA) drugs. Employing machine learning, we classify an AA phenotype from extracted features measured during and after a startle battery in fish exposed to AA drugs. Using this classifier, we demonstrate that structurally independent JNK inhibitors replicate the AA phenotype with high accuracy, consistent with findings in mice. Furthermore, pharmacological targeting of JNK1-regulated signalling hubs identifies AKT, GSK-3, 14–3-3 ζ/ε and PKCε as downstream hubs that phenocopy clinically proven AA drugs. This study identifies AKT and related signalling molecules as mediators of JNK1-regulated antidepressant- and anxiolytic-like behaviours. Moreover, the assay shows promise for early phase screening of compounds with anti-stress-axis properties and for mode of action analysis.


Distance, turning, pausing and spurting are altered by AA drugs during the startle period
We next examined the behavioural syllables that make up the motility response.Thus, in addition to distance travelled, we also measured turning (sum of turning angles), pausing (total time spent pausing), spurt velocity, time thigmotaxis and distance thigmotaxis near the well border (Fig. 2A).We characterised the effect of AA drugs on this expanded set of features using 411 zebrafish larvae behaviour tracks (Fig. 2B).Fluoxetine, imipramine, LiCl and diazepam elicited strikingly similar behavioural syllables during the startle period (Fig. 2B).All clinical drugs tested increased turning behaviour and decreased thigmotaxis distance and time, and spurting behaviours at lower doses.Only diazepam (at 1 µM or higher) diverged somewhat in that thigmotaxis distance was increased.Ketamine mimicked MK801 at high dose (100 µM) with regards to thigmotaxis (which increased), whereas turning increased with 10 µM ketamine or 1 µM MK801, consistent with the effect being mediated by the NMDA receptor for which MK801 has higher affinity.Notably, increased turning was a prominent feature of the AA-drug treated fish (Fig. 2B).In addition, MK801 induced a large increase in motility at 10 µM consistent with the known hyperlocomotor effect of NMDAR inhibition in rodents and fish 42,[44][45][46] , however the effect of ketamine on distance moved was minor in comparison during the startle period (Fig. 2B and Supplementary Fig. 5A, B and 6A, B).
In order to validate that the test battery was sensitive to a HPI stressor, we investigated the response to elevated saline, which is known to activate the HPI axis and increase cortisol in zebrafish larvae 47,48 .Moreover, it has been shown that NaCl or white light-induced motility responses require the glucocorticoid receptor in zebrafish larvae 48 .Interestingly, NaCl stress (which freely crosses the blood brain barrier) is also a classical means with which to activate JNK in the nervous system 28 .We treated fish with 100 mM NaCl, or with E3 medium alone, and measured the motility over 30 min (Fig. 2C).Saline treatment increased motility as previously 46 , and altered distance, spurting and pausing in a JNK-dependent manner, as responses were reversed upon JNK inhibitor treatment (SP600125, 10 µM for 1 h) (Fig. 2D).We next tested the effect of saline stress in the test battery.NaCl-induced behaviour syllables that were generally in the opposite direction of those elicited by anxiolytic www.nature.com/scientificreports/and antidepressant drugs (Fig. 2E).This indicates that the battery test replicates stress responses that are JNK and cortisol-dependent.

Identification of JNK1 pathway signalling hubs
Having established phenotypic responses to clinically used AA drug classes in zebrafish larvae, we next turned our attention to the JNK1 pathway.As JNK1 regulates many physiological processes in brain, we were interested to identify downstream signalling hubs that could serve as alternative, possibly more specific AA drug targets than JNK1 itself.Initially we identified JNK1-regulated phosphoproteins by comparing the phosphoproteome from wild-type (WT) and Jnk1-/-mouse brains.Phosphoproteins that differed significantly by more than 1.5 fold in Jnk1 knockout brain are presented as a circular array (Fig. 3A).We next predicted the most highly connected phosphoproteins to the JNK1-regulated ones using the GeneMANIA physical interaction network database.The predicted proteins that physically interact with the JNK1-regulated phosphoproteins are displayed at the centre of the circle (Fig. 3A).From the identified interacting proteins, we determined 29 candidate hubs from the JNK1-regulated phosphoproteome based on interaction counts and interaction confidence levels.To assess the significance of these interactions, we compared interaction counts for JNK1-regulated phosphoproteins to 1000 randomly generated datasets derived from the entire brain phosphoproteome (Fig. 3B).These hubs demonstrated significantly higher connectivity in Jnk1-/-mouse brains compared to WT brains.From among the hubs, www.nature.com/scientificreports/we further selected those for which small molecule drugs were available for testing in the zebrafish behavioural screen.The selected hubs included YWHAZ/E (14-3-3s), AKT1, GSK3B, BRAF, and PKCE/G.

JNK1 pathway hub compounds induce AA-like stereotypic behaviours during the startle period
We next screened JNK1 pathway drugs during the startle period.JNK inhibitors SP600125 at 1 µM, and JNK-IN-8 at 10 nM mimicked the AA drug profiles (Fig. 4A).Downstream of JNK1, we targeted signalling hubs 14-3-3 and AKT using R18, an inhibitor of 14-3-3 client binding 49 and SC79 (an AKT activator).This also induced a profile similar to that obtained with clinically used AA drugs.In contrast AKTi, an inhibitor of AKT, showed an opposite effect on distance and pausing and less effect on turning (Fig. 4A).GSK-3 inhibitor (SB216763) at 1 µM and LiCl at 1 mM, both known to inhibit GSK-3, produced similar responses.We also screened PKC-targeted drugs.PKC-epsilon activator (FR236924) 50 at 0.1 µM increased turning and pausing while decreasing distance (Fig. 4A).Bryostatin-1, a mixed PKC agonist, also increased turning, but differed from FR236924 in that it reduced spurting.Also, at 10 nM (a dose that activates PKCα and PKC δ isoforms, but antagonizes PKCγ), it increased thigmotaxis and decreased pausing.Conversely, the pan-PKC inhibitor 3,4-Bis(3-indoyl)maleimide (BIM) increased turning and decreased thigmotaxis at 100 nM, a dose that inhibits PKCα more effectively than PKCε 51 .

Drug behaviour profiles during the post-startle period
The post-startle period is also relevant for stress responses and habituation 52 , therefore we examined zebrafish behaviour during this period.Clinically validated drugs produced behaviours similar to those during the startle period but with a more pronounced reduction in distance moved, increased pausing, and milder effects on thigmotaxis (Fig. 5A).Turning behaviour during the post-startle response period increased more linearly with dose, and at 10 µM, ketamine's behaviour profile resembled classical AA drugs in this period.In contrast, saline stress strongly increased distance travelled and decreased turning and spurting at 100 mM NaCl (Fig. 5B).These behaviour profiles are consistent with a block of the HPI axis stress response by the clinical drugs, as expected.Among the JNK1 pathway hub compounds, several induced AA-like profiles during the post-startle period (Fig. 5C).For example, the 14-3-3 inhibitor R18 (0.1 µM) mimicked classical AA drug profiles, reducing distance moved while increasing turning and pausing.The selective PKCε activator FR236924 (0.1 µM) showed a similar effect, decreasing distance and thigmotaxis while increasing turning and pausing.Signalling hubs derived from the Jnk1-/-mouse brain phosphoproteome (labelled using gene names) represent the most highly connected phosphoproteins from among all phosphoproteins that are significantly altered in Jnk1-/-brain verses wild-type.For each hub, the distribution of physical interaction count from the 1000 background networks is represented with a boxplot.The number of physical interactions for the same hub in the Jnk1-/-mouse brain phosphoproteome (i.e. the phosphoproteins that are significantly altered in Jnk1-/-mouse brain phosphoproteome), is depicted with a red cross.The most highly ranked JNK1regulated signalling hubs are shown.www.nature.com/scientificreports/JNK inhibitor SP600125 (1 µM) reduced distance, decreased thigmotaxis, and increased pausing, while the less well characterised JNK-IN-8 at 0.01 µM, had no significant effect.The AKT activator SC79 (2 µM) produced a response similar to JNK inhibitor and AA drugs, while the AKT inhibitor AKTi (50 µM) induced an opposite response, possibly representing anxiogenic-like behaviours.
GSK-3 inhibitor SB216763 mimicked the response of AKT activator SC79.The PKC activators FR236924 (0.1 µM) and PMA evoked similar responses, although Bryostatin-1, which activates some PKC isoforms, differed substantially in profile 51 .Notably, we found that NMDA receptor antagonists MK801 and high-dose ketamine elicited similar effects to each other in zebrafish, as previously reported 53 .Interestingly, the behavioural profiles evoked by ketamine and MK801 after the startle period (rather than during) (Fig. 5), aligned more closely with those of the AA drugs.AKTi, BIM, and Bryostatin-1 did not fit the AA phenotype, as expected, whereas haloperidol fit both "AA" and "other" classification.Finally, saline stress (NaCl) induced responses during the post-startle period were measured.These were notably similar to those exerted by AKTi (Fig. 5B), consistent with an anxiogenic profile.

Machine learning identifies compounds with antidepressant/anxiolytic classification
In addition to traditional statistical analysis, we employed machine learning to create a comprehensive AA phenotype classification.For this we used Receiver Operating Characteristic (ROC) curves to assess the "fit" of test drugs to this classification.This approach has the advantage of overcoming the limitations and biases of regular statistical analysis and enhances the interpretation of complex data.We tested three supervised learning models: Random Forest classification-regression, GLMNET linear-regression, and SVM non-linear-regression.The SVM non-linear regression model achieved the highest training accuracy (> 0.99) for both the startle and post-startle response periods with AA drugs, making it the chosen model to validate the test compounds.
When using post-startle response data, the accuracy of predictions for test drugs was generally higher.Thus, SC79, R18, SB216763, PMA, and FR236924 scored highly for an AA phenotype in the post-startle period, even surpassing the scores of JNK inhibitors, suggesting that greater functional specificity is achieved by targeting www.nature.com/scientificreports/downstream of JNK1.A summary of the multiparametric zebrafish larvae behaviour test and the impact of signalling hubs downstream of JNK on the AA phenotype are outlined in Fig. 6C.

Discussion
In this study, we establish a multiparameter behavioural screen to classify a common AA phenotype in zebrafish larvae using clinically proven AA drugs.We use a range of classical antidepressant and anxiolytic drugs; fluoxetine, representing a commonly used serotonin reuptake inhibitor; imipramine, a tricyclic antidepressant, and ketamine, a fast-acting antidepressant, as well as the mood stabilizing drug lithium chloride and diazepam, a commonly used anxiolytic of the benzodiazepine class.Despite having distinct pharmacological profiles, these drugs elicit similar effects on zebrafish larvae motility patterns.All AA drugs reduce distance moved, thigmotaxis www.nature.com/scientificreports/and spurting, while they increase turning.This differs from the behavioural responses to psychosis-inducing drugs MK801 and high dose ketamine NMDA receptor antagonist drugs and anti-psychotic haloperidol, and saline-induced HPI axis stress, all of which produce opposing behaviours in the same test battery.We employ machine learning to characterise an "AA phenotype" using all of these behavioural features, and go on to use the assay to evaluate the signalling hubs downstream of JNK1.The results uncover new regulators downstream of JNK1 that align with the classical AA-drug phenotype performing even better than JNK inhibitors, suggesting increased functional specificity.Dysregulation of the hypothalamic-pituitary-adrenal axis leading to hypercortisolemia is one of the main circuit irregularities contributing to anxiety disorders and major depression [54][55][56][57] .This hormonal change can trigger a cascade of events including the activation of JNKs.Accordingly, JNK is activated by cortisol and induces pro-inflammatory cytokines as well as amygdaloid dendritic hypertrophy 32,[58][59][60] .Additionally, JNK plays a role in the homeostatic regulation of glucocorticoid receptor transcriptional activity 28,31 .In contrast, when JNK is inhibited in mouse brain a reduction in anxiety-and depressive-like behaviours is observed 61 , as is the case with JNK inhibitor treatment or genetic deletion 36,[62][63][64] .In contrast, chimeric JNK activation induces depressive-like behaviour and impaired assessment of risk/reward benefit 64,65 .In zebrafish brain at the larval stage, JNK1 orthologues (mapk8a and mapk8b) are expressed 66 , and consistent with behavioural findings in rodents 60 , we find that treatment with structurally independent JNK1 inhibitors SP600125 and JNK-IN-8, induce an AA phenotype with A.U.C.s of 97.9% and 98.3% respectively in the startle period.Thus, JNK inhibitors elicit a behavioural response in zebrafish larvae that matches closely to that of clinically used AA drugs.
A main goal of this study was to identify signalling hubs downstream of JNK1 influencing anxiety-and depressive-like behaviours.We discover 14-3-3ζ (YWAZ) and 14-3-3ε (YWHAE), members of the 14-3-3 client binding family 67 , as key hubs in this process.This is interesting as JNK is known to phosphorylate 14-3-3ζ leading to client protein release 68 .Consistent with this, in Jnk1-/-brain, we detect an increase in the interaction count for ζ and ε 14-3-3 isoforms, indicating that JNK1 regulates 14-3-3 binding dynamics in the CNS.In zebrafish larvae, we show that R18 peptide which disrupts 14-3-3 client binding 69 , recapitulates the AA-drug phenotype with A.U.C. of 93.5% and 97.2% in the startle and post-startle period respectively, consistent with the involvement of ζ or ε isoforms 67,70 .14-3-3ε has also been shown to mediate chronic stress-induced depression 71 .Interestingly, fluoxetine rescues the behavioural deficit in zebrafish lacking 14-3-3 ζ suggesting a shared mechanism with fluoxetine 72 .
We also identify AKT as a potential signalling hub downstream of JNK1 based on mouse brain phosphoproteome interaction network analysis.AKT and 14-3-3 are mechanistically linked, thus inhibition of 14-3-3 client binding activates AKT 73,74 , potentially placing AKT regulation by JNK1 downstream of 14-3-3.We find that the AKT activator SC79 evokes an AA phenotype with A.U.C.s of 99.4% and 98.3% for the startle and post-startle response periods in the zebrafish larva screen.Conversely, the allosteric inhibitor of AKT (AKTi) induced an opposite response (A.U.C. 60.3% and 60.4% respectively).That AKT activation induces an AA phenotype fits with the well-documented antidepressant action of AKT in animal models and depressed subjects, and its regulation of BDNF 75,76 .Consistent with this, we also identify that GSK-3β is a hub downstream of JNK1 in mouse brain.This is not surprising as GSK-3α/β isoforms are phosphorylated directly by AKT on S21 and S9 (α and β isoforms respectively), leading to GSK-3 inhibition 77 .In the zebrafish larva assay, the GSK-3α/β inhibitor SB216763, elicits an AA phenotype with A.U.C.s of 97.6% and 99.0% during the startle and post-startle response phases.
B-Raf was another signalling hub identified downstream of JNK1.In the interaction screen, B-Raf displayed an increase in protein interactions (from "0" to "10") in the phosphoproteome from Jnk1-/-mouse brain compared to wild-type.However, when we tested the B-Raf inhibitor SB590885, the resulting phenotype did not classify as "AA" either during the startle or post-startle period.Interestingly, B-Raf is also known to interact with 14-3-3, which functions to maintain its kinase domain in an inactive state 78 .The expectation therefore would be that activation of B-Raf upon disruption of 14-3-3 binding may allow it to contribute to an AA-like state, as observed with the 14-3-3 client binding inhibitor R18.However, it was not possible for us to test B-Raf further due to lack of activators for this molecule.The oncogenic nature of B-Raf activating mutations contraindicates the use of such molecules.Finally, we tested PKC activators and inhibitors in the zebrafish screen to explore their impact as possible effector hubs downstream of JNK1.We used the PKCε-activator FR236924 and the pan-PKC activator PMA, both of which evoked a clear anxiolytic-like effect (A.U.C.s of 97.8% and 95.7% respectively).Consistent with this, the broad specificity PKC inhibitor (BIM) had the opposite effect.Notably however, in some contrast with our findings in zebrafish larvae, BIM treatment attenuates corticotrophin releasing factor (CRF) facilitation of the acoustic startle response in mice, although it had no effect in the absence of CRF 79 .Moreover, PKCγ-/-mice exhibit low anxiety 80 .Interestingly also, that PKCε activator (FR236924) promotes an AA phenotype in zebrafish is supported by the response to byrostatin-1 at 100 nM, (a dose that activates PKCε 81 ) during the startle response, even though when analysed at the full dose response in our study the overall phenotype is "other".While these studies suggest isoform-specific roles for PKC variants in anxiety-like behaviours, we identify that PKCε represents a signalling hub downstream from JNK1 which controls an AA phenotype in zebrafish larvae, indicating that further analysis of this isoform in the context of stress-related behaviours is warranted.
Our observation that distinct clinical drugs, namely fluoxetine, imipramine, diazepam and LiCl, elicit overlapping behaviours in this short timeframe assay merits some discussion.Although these drugs share some pharmacological endpoints, their primary mode of action is classically seen to differ.For example, fluoxetine and imipramine block serotonin uptake, whereas imipramine also blocks norepinephrine uptake and antagonises dopamine D2 and acetylcholine receptors, and the primary target of diazepam is GABA receptors 82 .Moreover, these drugs elicit effects on neurotransmitter levels both in responders and non-responders, giving rise to the theory that other effects such as neuroplasticity and neurogenesis changes may be needed for the therapeutic effect 82,83 .A recent groundbreaking study identified the Tropomyosin receptor kinase (TrkB) as a common binding site for antidepressant drugs including fluoxetine, imipramine and ketamine.The authors demonstrated that binding of these drugs to TrkB occurs within minutes and is required for the antidepressant-like effect 84 .They rationalized that neuroplasticity effects of BDNF, mediated via its receptor TrkB, could explain the antidepressant effect and its timeline.In our experiments, many lines of evidence support a common role for TrkB in the observed responses.Firstly, it is feasible, as zebrafish larvae express TrkB orthologs (trkB1 and trkB2) from 6 days post fertilization onwards 85 .Secondly TrkB activation upon BDNF binding leads to rapid (within minutes) activation of AKT 75 , and the behaviors measured here occur after only 60 min drug exposure.Thirdly, we show that JNK1-regulated signalling hubs identified here converge on the AKT pathway which incorporates 14-3-3 ζ/ε and GSK-3.Finally, the AKT activator (SC79) evokes an AA phenotype that closely matches the clinically used drugs with an A.U.C. of 99.4%.Thus, AKT pathway activation offers the most likely explanation for the temporal and phenotypic overlap observed in response to the clinically used and JNK1-regulated hub drugs.
Animal models are limited when it comes to modeling complex disorders such as anxiety and depression.They do not account for subjective experience, or complex genetic and environmental background, and the models that do exist in mice require large numbers to attain statistical power 86,87 .Here we present an assay that is more straightforward and scalable than rodent models.It measures stereotypic responses to acoustic and light stress, based on larval behavioural sequelae (distance, turning, pausing, spurting, thigomataxis distance, thigmotaxis time, dose and size effect).This may lack face and construct validity, and it is not designed to model a specific disease mechanism.Nonetheless, it does measure HPI axis stress responses and the neuroendocrine and monoamine control systems underlying the stress axis is conserved from chordates to mammals 88 , moreover locomotor activity in zebrafish larvae is inherently cortisol dependent 89 .A similar approach has been used to identify novel anti-psychotic compounds in zebrafish larvae 16 , and to investigate the paradoxical excitatory effects of GABA and serotonin 17 , for example.Interestingly, here we also observe a similar hyperactivity effect with the GABAR agonist diazepam when fish are exposed to violet light and acoustic stimuli together (Fig. 1H).Our study in addition exploits machine learning to classify the desired phenotype, enabling use of complex combinations of features.The result is high accuracy prediction scores for test drugs that are known to be anxiolytic/antidepressant in mice e.g.AKT activator and GSK-3 inhibitor 75,76,90 .This highlights the potential of the zebrafish larvae screen to potentially identify compounds that mimic AA drugs, for further testing in rodent models, which is rather the goal of a phenotypic screen.
A possible explanation for the high prediction score could be linked to the dependence of the zebrafish larvae locomotor response on the HPI-axis 48 , and the drugs used regulate the HPA-axis 91 .Although the stereotypic behaviours measured here may not directly relate to AA conditions, JNK is activated by cortisol 37,92 leading to synapse withdrawal 32 , which may influence circuit behaviour and long-term mental state.Notably also, the phenotype of AA-drugs differs from those of psychotic compounds MK801 and ketamine.Similarly, the antipsychotic D2 antagonist haloperidol shows dual classification which is consistent with its known paradoxical side effect, whereby some patients can respond with worsening of symptoms such as increased agitation, insomnia and hallucinations in patients 89 .Consistent with this in the zebrafish larvae, we find that haloperidol shows a distinct hyper-motility phenotype, which possibly accounts for its classification as "other"; yet it also increases turning and diminishes spurting, consistent with the "AA" phenotype (Fig. 2; Supplementary Fig. 7).This paradoxical effect to haloperidol was also previously reported in zebrafish larvae 16 .A long-term advantage of the screen is the potential reduction in the number of mice needed for drug discovery, while expanding the utility of zebrafish larvae in neuropsychopharmacology studies.

Materials and methods
Zebrafish husbandry-WT zebrafish (Danio rerio) larvae were derived from a cross between the high-performance AB strain and the Tu short-fin wild type strain 22 .Fish were maintained in the Turku Bioscience Zebrafish Core facility at 28.5 °C, under a 12 h day/night cycle.Embryos were cultured in E3 medium (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl 2 , 0.33 mM MgSO 4 ) with 30 eggs per 10 cm diameter dish.Larvae of indeterminate sex were used from 3 to 8 dpf as indicated.Experiments were carried out in the afternoon to maximise consistency of locomotor behaviour 93 .Fish were bred in two tanks.At least two plates of fish were tested for every drug.There were minor differences in baseline motilities between plates therefore data was normalised to internal control fish measured on the same plate.After the experiment, zebrafish larvae are anaesthetised with 200 mg/L Tricaine for 5min, and then euthanized by transferring to 1% chloramine ice-water solution.Behavioural experiments were carried out following the guidelines that were approved by the National Animal Experiment Board (ELLA) and the Regional State Administrative Office of Southern Finland (ESAVI) under licenses ESAVI/16343/2019 and ESAVI/30686/2022, which is the authority in Finland that ensures that animal experimentation is carried out following the law and thus follows EU ethical guidelines for animal experimentation.All methods are performed and reported in accordance with the ARRIVE guidelines (https:// arriv eguid elines.org).
Behavioural battery test-7 days post fertilization (dpf) zebrafish larvae were transferred to square well 96-well plates and drug dosages were administered following which the plate was directly transferred to the DanioVi-sion™ device (Noldus Information Technology, BV, Wageningen, Netherlands) and acclimatized for 1 h in the dark at 28.5 °C.Larvae then underwent the test battery consisting of 5 cycles of a 60 s string of acoustic/vibratory (tapping) and light stimuli programmed using Ethovision13 XT software (Noldus Information Technology, BV, Wageningen, Netherlands).The 60 s stimuli were as follows; high (H) and low (L) intensity patterns (1.6 & 4.2 W): 1 s pause, 6 × L-tap (0.5 s break between individual taps), 4 s pause, 4 × H-tap, 1 s pause, 8 × L-tap, 8 s pause, 2 × L-tap, 10 × H-tap, 7 s pause, 2 × H-tap, 13 s pause, 14 × H-tap, 3 s pause.Visual stimuli evoked using red (635 nm, 192 Lm), blue (470 nm, 240 Lm), purple (red and blue together) and white LEDs, switched colour at 10 s intervals in the following order: darkness, red, blue, purple (red & blue), white, and sequential flickering light of all colours sequentially at 0.1 to 0.2 s intervals.After the stimuli period there was a 29 min intermission period.This 30 min sequence was repeated 4 more times.Digital video tracking was performed using the Daniovision imaging unit using a pre-programmed protocol with the included Ethovision XT software.Acquisition was 60 fps and exported as csv file for automated analysis using customized software in R. The data analyst was agnostic to the drug properties.
Stress response assay-5 or 7 dpf larvae (as indicated) were transferred to a 48-well plate with 1 fish/well and acclimatized for 1 h following which they were exposed to stress: either with white light illumination (10,000 Lux) for 50 s or addition of 100 mM NaCl.Larvae were tracked for a further 30 min following stressor in the dark using infrared light.NaCl experiments were performed using ambient level (500 Lux) white light.
Drug treatments-Larvae at the indicated ages were treated with increasing doses of drugs and compounds as follows: fluoxetine, diazepam, FR236924 and PMA (Phorbol 12-myristate 13-acetate) and ketamine were from Tocris Bioscience, Abingdon, UK; haloperidol, imipramine, LiCl, SP600125, SB590885, Bryostatin-1 and MK801 from MERCK Sigma Aldrich Solutions, Darmstadt, DE; R18 from Enzo Life Sciences distributed by AH Diagnostics, Helsinki, FI; SB216763, JNK-in-8 and AKTi from Selleck Chemicals LLC, Houston, Texas; bisindolylmaleimide I (BIM) from Santa Cruz Biotechnology Inc., Dallas, Texas and SC79 from R&D systems, McKinley Place NE, Minneapolis or NaCl in E3 medium.Treatments were for 1h at 28.5°C before commencing the startle battery.Control fish were treated with equivalent volume of carrier.
Behaviour data analysis-Motility: Raw motility tracking data was exported in csv format, where each file contained data from 1 experimental cycle, and each worksheet contained time series data from 1 fish.Relative motility index (M.I.R ) (plots of left side) represents the motility during a certain condition (Mc x ) relative to the motility during the control condition (Mc 0 ).The relative motility index was calculated as follows from all 5 cycles: "Relative Peak Analysis" (shown in Fig. 1F-K right side panels) represents the motility during peaks P1, P2, P3, P4, P5, and P6 of a test condition relative to the same peaks from control fish in the same plate.The data scientist handling the analysis was agnostic to the drug treatments.Analysis was done in R.
Analysis of other features: X and Y positions, heading, and distance features from the time series data were utilized to calculate behavioural features: distance, turning, pausing, spurt velocity, and thigmotaxis by distance and time.Distance -Distance was summed from all time points per fish and mean total distance travelled per treatment dose was calculated.Turning -Turning angles between consecutive time points were calculated from the heading parameter of the last and current time points, and converted to magnitude representing the smaller angle in 360 degrees.Relative turning angles were calculated by summing all time points for each fish, then normalizing to distance travelled, and presented as average values according to drug dose.Pausing-Time points with value 0 for distance were summed to give a value for total pause time for each fish.Spurt velocity-Time bouts between pauses were identified and the velocity for each time point was calculated from distance feature.The maximum velocity from each time bout was recorded and the mean maximum velocities from all time bouts in a single fish were calculated.Thigmotaxis-The well border coordinates from the 96-well plate were recorded from Ethovision graphic user interface (GUI).An inner border was (virtually) created that was 40% of the width and length of the original border, and placed in the centre of each well.The coordinates of these inner borders served as thresholds for thigmotactic behaviour.Time points where a fish travelled outside the inner border were identified and the ratio of near outer border travel versus total travel for each fish was calculated for time and distance.Analysis was done in R.
Machine learning analysis-Machine learning algorithms were applied to zebrafish larvae behaviour data during the 1 min startle period and during the 10 min of the recovery period.Specifically, the quantitative variables total pause duration, total distance travelled, relative turning angle, average spurt speed, total thigmotaxis time and total thigmotaxis distance were normalized against control (or zero-dosage fish) data as percentages.These normalized metrics, along with dosage level variable, were incorporated as features for the machine learning models.The models were tasked to differentiate between "antidepressant or anxiolytic" drugs and "other" drugs.
Vol.:(0123456789) were categorized as "antidepressant/anxiolytic" class, while ketamine and MK801 treatments were categorized as "other".Generalized Linear Model via penalized maximum likelihood (Glmnet), Support Vector Machine (SVM) with a polynomial kernel, and Random Forest algorithms were utilized to build 3 separate models each for startle response behaviour (SRB) data during the 1 min startle battery and post-startle response behaviour (PSRB) data during the 10 min recovery period.The constructed models were subsequently employed to predict the class category of other drug treatments.Machine learning was performed with "caret" package within the R environment.
Statistics-All behavioural features were compared to control fish from the same experiment.For statistical testing, feature values were collected into a distribution, and this data according to drug dose was compared to respective control distributions using Wilcoxon signed-ranked test or one-way or two-way ANOVA with Tukey or Dunnett's multiple comparison test as indicated.Calculations and plots were made using a customized pipeline in R programming language.Data analysis was automated and the data scientist was agnostic to the treatment drugs.
Preparation of samples for mass spectrometry, SDS-PAGE and in-gel digestion -All the chemicals for digestion and mass spectrometry analysis were from Sigma-Aldrich (Stockholm, Sweden).Whole brain was isolated from WT and Jnk1-/-C57B6J mice 94 at embryonic day 15 (E15), post-natal day 0 (P0), post-natal day 21 (P21), and 8 months (Adult).We extracted brains from three mice for each age and genotyped and snap froze them in liquid N 2 .Later tissue was pulverized in a Micro-dismembrator II (Braun Biotech International, Melsungen, Germany).We weighed the powder and stored aliquots at -80°C before use.SDS (1%), supplemented with protease and phosphatase inhibitors (Sigma-Aldrich, St. Louis, USA), was added to the brain-powder.Mechanical disruption followed using a syringe twice.We centrifuged samples for 1 h at 4 °C and collected the supernatant.We quantified protein concentration using the Total Protein Kit, Micro Lowry, Peterson's Modification (Sigma-Aldrich, St. Louis, USA).
A total protein amount of 100 μg for each sample was separated on 12% Criterion gels (Bio-Rad Laboratories, Hercules, Ca, USA), following the manufacturer's instructions.We washed the gel in milliQ, stained for 30 min with GelCode (Bio-Rad Laboratories) and washed overnight in milliQ before manually slicing each lane into 5 equal slices.We de-stained gel slices by washing 3 times in 25 mM ammonium bicarbonate (AMBIC)/50% acetonitrile and dried in a vacuum centrifuge (Speedvac).Samples were reduced (10 mM DTT/100 mM AMBIC, 1 h at 56 °C), alkylated (55 mM iodoacetamide (IAA)/100 mM AMBIC, 45 min at room temperature in the dark), washed 2 times with 100 mM AMBIC and dehydrated with ACN before being dried in the Speedvac.We rehydrated slices with 12.5 μg/ml modified porcine trypsin (Promega, Madison, WI, USA) in 50 mM AMBIC and digested overnight at 37 °C.Peptides were eluted 2 times in 75% ACN/5% FA, dried 1 h in the Speedvac and dissolved in 0.1% formic acid.
TiO 2 phospho-peptide enrichment -we incubated samples for 10 min with Buffer A (1 M GA, 80% ACN, 5% TFA) followed by incubation with the TiO2 Mag Sepharose (GE Healthcare Life Sciences) for 30 min in an Eppendorf Thermomixer at 800 rpm.Liquid was removed and washing steps were done using Buffer A and 2 times using Buffer B (80% ACN, 1% TFA).Phospho-peptides enriched were eluted 2 times using 100 mM ammonium hydroxide in water and dried 1 h in the Speedvac.
LC-MS/MS analysis-Dried peptides were resuspended in 0.1% formic acid and separated using an Eksigent nano-LC 2D system (Eksigent, Dublin, CA, USA) consisting of a solvent degasser, a nano-flow pump and a cooled auto-sampler.Peptide concentrations were measured by Nanodrop (ThermoFisher, Stockholm) and the concentration adjusted to the same amount in all samples.Eight μl of sample was loaded and washed for 15 min onto a pre-column (Acclaim PepMap 100, C18, 3 μm particle size, 50 μm diameter, Thermo Fisher Scientific, Hägersten, Sweden) at a constant flow of 5 μl/min solvent B (0.1% FA in ACN).We ran three biological replicates and two technical repeats for each of the genotypes across four time points.Phosphopeptides were loaded into a RP analytical column (10 μm fused silica emitter, 75 μm ID column, 16 cm Pico Tip, New Objective) packed in-house with C18 material ReproSil-Pur and separated using an eighty-minute gradient from 3 to 40% solvent B at a flow rate of 300 nl/min.The gradient was followed by 20 min column washing with 90% ACN, 0.1% FA and 15 min re-equilibration with 3% solvent B. The HPLC system coupled to an LTQ-Orbitrap XL mass spectrometer (ThermoFisher Scientific, Bremen, Germany) operated in Data Dependent Acquisition mode.Spray voltage was set to 1.90 kV and the temperature of the heated capillary was set to 200 °C.The ten most intense ions from the survey scan performed by the Orbitrap were fragmented by collision-induced dissociation (CID) in the LTQ (normalized collision energy 35, parent mass selection window 0.5 Da, 30 ms activation time and minimum signal threshold for MS/MS scans set to 100).We excluded unassigned charge states and singly charged ions from fragmentation.The dynamic exclusion list was limited to a maximum of 500 masses with a retention time window of 2 min and a relative mass window of 10 ppm.The X-calibur software version 2.0.7 (Thermo Scientific) controlled the HPLC, mass spectrometer and data acquisition.
Shotgun data analysis from LTQ Orbitrap -The raw data was uploaded into Progenesis LC-MS (Waters) for analysis.After raw mass spectra were aligned and normalisation across all runs was carried out by total ion current in the analysis area manually defined for the reference HPLC for the time and mass ranges.We submitted the raw files from the mass spectrometer to Thermo Proteome Discoverer (version 1.0.43.0) for protein identification by Mascot.Cysteine carbamidomethylation was set as a fixed modification, methionine oxidation and phosphorylation of serine, threonine, tyrosine (STY) were set as variable modifications.Data were searched against Uniprot Mus musculus database containing 170.598 sequences (release date April 2014).The mass tolerance was set to 10 ppm for the precursor ion and 0.8 Da for the fragment ions.Trypsin was set as the protease and a maximum 2 of missed cleavages were allowed.We ran raw data using an FDR of 5% and a score cutoff of 20.We discarded proteins marked as contaminants, reverse hits and "non-unique peptides".We selected phospho-peptides with a Mascot score higher than 20 and verified using the default Mascot Delta score of 0.1 for comparative analysis between groups.www.nature.com/scientificreports/Data merging and statistical analysis phosphoproteomics data -Downstream analysis of MS intensity outputs from Progenesis LC-MS (Waters) utilized a customized R (v.3.4.0) software pipeline.Phosphopeptide data were refined as follows: (i) all phosphopeptides from 5 gel slices were merged, (ii) identical phosphopeptides, i.e. those with identical sequence, UniProt identifier and number and position of phosphorylated residues, were merged by summing the intensities.(iii) If a peptide entry had > 1 missing value across 3 replicates, all entries were removed.For quality control of data, a distribution plot of phosphopeptide intensities for each unique genotype at given ages was visualized using boxplot and histogram analysis, using ggplot2 95 (Supplementary F1)."Phosphopeptidecentric" and phosphosite-centric analysis was performed as appropriate and is defined in figure legends.We used the Bioconductor package RankProd 96,97 for statistical analysis.We calculated mean ratios of phosphopeptide intensities [Jnk1-/-/WT] per entry (for both phosphopeptide-centric and protein-centric changes).Only phosphorylation intensity changes that passed a threshold Ratio Jnk1-/-/WT > 2, with p-value < 0.05 were included.
Hub analysis -29 candidate hubs (Fig. 3A,B)that were significantly altered in Jnk1-/-mouse brain compared to wildtype brain and were part of the psychiatric disease associated genes according to the MetaCore™ database schizophrenia gene list were selected to analyse as potential signalling hubs acting downstream of JNK1 in regulating depressive/anxiety-like behaviour.The selection criteria included a threshold for phosphoproteins that showed either ≥ 7 physical interactions, or > 3 high weight (weight > 0.04) physical interactions from the physical interaction network built with GeneMANIA from among the significantly altered phosphoproteins derived from the wildtype and Jnk1-/-mouse brain phosphoproteomes.29 proteins fulfilled these criteria.To control for hub protein's inherent functional plasticity variations, 1000 control phosphoprotein networks of the same size as the Jnk1-/-psychiatric disease phosphoprotein network were randomly generated from the entire brain phosphoproteome dataset using GeneMANIA (v.3.4.1).The number of interactions of the 29 proteins from the significantly altered phosphoproteome network in Jnk1-/-mouse brain was each compared to the distribution of interaction count from 1000 randomly generated control networks to determine the significance of the hubs.Statistical significance for the hubs was calculated using the hypergeometric test.P-values were adjusted using the Benjamini-Hochberg procedure.
Receiver Operating Characteristic (ROC) analysis-ROC curves were constructed to visualise the predictive proficiency of the machine learning models.Individual curves were plotted for each test drug from the test dataset, paired with a control drug that was either classified as a "positive control" or a "negative control".Each treatment drug is assigned with a "true" class label that reflected its known properties and molecular interactions.Drugs of the "antidepressant or anxiolytic" class were coupled with the "negative control" drug ketamine, which provided the "true" class label of "other".Conversely, drugs classified under the "other" class were paired with the "positive control" drug fluoxetine for the "true" class label of "Antidepressant or anxiolytic" for the purposes of ROC curve calculation.The construction of the ROC curve was facilitated by the 'pROC' package within the R statistical programming environment.

Figure 2 .
Figure 2. Testing the effect of AA drugs on zebrafish larvae behavioural sequalae during the 1 min startle phase.(A) The behavioural features extracted from zebrafish larvae tracking are shown.Distance, turning, pausing, spurting, time and distance thigmotaxis are extracted using R statistical computing platform.(B) Fish were exposed to AA drugs and ketamine, MK801 or haloperidol as indicated and features (distance, turning, pausing, spurting, time and distance thigmotaxis) were extracted from the entire 1 min startle period.Averaged data indicates change relative to control for each of these behaviours according to colour code.(C) Mean motility profiles of zebrafish larvae (5 dpf) before and after E3 (n = 24) or NaCl (100 mM, n = 24) are shown.(D) Mean distance travelled, turning, pausing, and spurting of zebrafish larvae (5 dpf) before and during the 10 min after 100 mM NaCl (n = 72) or control (n = 67), are plotted as a % change from control.(E) Zebrafish treated with E3 medium or NaCl as indicated underwent the battery test.Mean data on distance travelled, turning, pausing and spurting are shown.Control: 15 (75 observations), NaCl: 32 (156 observations).P-values were calculated by Wilcoxon Rank Sum test and adjusted with Benjamini-Hochberg procedure, where * p-value ≤ 0.05; ** p-value * ≤ 0.01; ***p-value ≤ 0.001; ****p-value ≤ 0.0001.

Figure 3 .
Figure3.Identification of signalling hubs downstream of JNK1 in mouse brain.(A) JNK1 regulated phosphoproteins from Jnk1-/-mouse brain are organised in a large outer circle with connections to the predicted interacting proteins in the centre.Warmer colours indicate a larger number of physical interactions, red being the highest.(B) Signalling hubs derived from the Jnk1-/-mouse brain phosphoproteome (labelled using gene names) represent the most highly connected phosphoproteins from among all phosphoproteins that are significantly altered in Jnk1-/-brain verses wild-type.For each hub, the distribution of physical interaction count from the 1000 background networks is represented with a boxplot.The number of physical interactions for the same hub in the Jnk1-/-mouse brain phosphoproteome (i.e. the phosphoproteins that are significantly altered in Jnk1-/-mouse brain phosphoproteome), is depicted with a red cross.The most highly ranked JNK1regulated signalling hubs are shown.