Introduction

A central function of the brain is to learn and change with experience. These adaptations can reflect attempts to identify and attend preferentially to salient stimuli. For example, identifying the smell of a predator or prey may be crucial, while identifying that my home still smells like my kin is not. This ability to suppress responses to continuous non-salient stimuli is known as habituation, a process generally considered to be the simplest form of learning and memory (Rankin et al., 2009). Habituation is conserved across all animals, and like other forms of plasticity, exists in at least two mechanistically distinct forms: transient short-term habituation, and protein-synthesis dependent long-term habituation. Here we focus on long-term habituation, which serves as a pragmatic model to dissect plasticity processes in neural circuits.

Work on long-term habituation in various species and paradigms has led to significant insights into the adaptations underlying this process (Cooke and Ramaswami, 2020; McDiarmid et al., 2019b), nonetheless a consensus model on the general principles underlying habituation is yet to emerge. Physiological and genetic work in Aplysia, and C. elegans were consistent with a model in which homosynaptic depression of excitatory synapses drives habit-uation (Bailey and Chen, 1983; Rose et al., 2003) (although see (Glanzman, 2009)). In contrast, work in the Drosophila olfactory and gustatory systems indicate that the potentiation of inhibitory neurons drives habituation rather than depression of excitatory connections (Das et al., 2011; Paranjpe et al., 2012; Trisal et al., 2022), and habituation to specific orientations of visual cues in mice is associated with the potentiation of neuronal activity and synapses in the visual cortex (Cooke et al., 2015), which requires GABAergic interneurons (Kaplan et al., 2016; Hayden et al., 2021). These studies are more consistent with a model in which the potentiation of inhibition, rather than depression of excitation, drives habituation (Cooke and Ramaswami, 2020).

Recently, we found that long-term habituation of the response of larval zebrafish to sudden pulses of whole-field darkness, or dark flashes (DFs), involves multiple molecularly independent plasticity processes that act to suppress different components of the behavioural response (Randlett et al., 2019). Similar behavioural, pharmacological, and genetic experiments have led to comparable conclusions in acoustic short-term habituation (Nelson et al., 2023), and habituation in C. elegans (McDiarmid et al., 2019a,b), indicating that habituation generally acts via multiple modular plasticity processes. These modules act to mute or shift behavioural responses to repeated stimuli. How and where these processes are implemented in the circuit, and how conserved or derived these processes are across species or paradigms remains to be determined. Here we have used a combination of high-throughput behavioural analyses, pharmacology and Ca²⁺ imaging to dissect DF habituation. Our results are consistent with a model in which habituation results from a multidimensional and distributed plasticity process, involving multiple independent molecular mechanisms. We propose that GABAergic inhibition is central to DF habituation, but how individual cell types and molecular events lead to behavioural adaptations during habituation will require targeted genetic and cellular approaches.

Results

Volumetric 2-photon Ca²⁺ imaging of habituation learning

When stimulated with a dark flash (DF), larval zebrafish execute an O-bend response (Figure 1A). The O-bend is characterised by a strong body bend and a large turn that forms part of the phototactic strategy of larval zebrafish, helping them navigate towards lit environments (Burgess and Granato, 2007; Chen and Engert, 2014). When presented with repeated DFs, larvae habituate and reduce their responsiveness, remaining hypo-responsive for multiple hours (Figure 1B), (Randlett et al., 2019).

Figure 1.

Figure 1—figure supplement 1. Validation of motion analysis based on image artifacts during 2-photon imaging

To explore the circuit mechanisms leading to this form of habituation, we asked how individual neurons within the DF responsive circuit adapt to repeated dark flashes. We used a head-fixed paradigm to perform 2-photon Ca²⁺ imaging in larvae expressing nuclear-targeted GCaMP7f pan-neuronally. Imaging was performed with a resonant scanner and piezo objective, enabling us to cover a volume of ≈ 600 × 300 × 120 μm (x,y,z) sampled at 0.6 × 0.6 × 10 μm resolution, leading to the detection of 30890±3235 ROIs per larvae (±SD, Figure 1C-E). ROIs were aligned to the Z-Brain atlas coordinates (Randlett et al., 2015), demonstrating that this volume spans the majority of the midbrain, hindbrain, pretectum and thalamus (Figure 1C-E).

We focused on a single training block of 60 DFs to identify neuronal adaptations that occur during the initial phase of learning (Figure 1Bi). This paradigm induced strong Ca²⁺ activity in neurons (Figure 1F), some of which were clearly associated with the DF stimuli. Ca²⁺ transients in response to DFs generally decreased across the 60 stimuli, though this pattern was not seen in all neurons, and substantial heterogeneity in their adaptations were observed. Strong correlated patterns were also seen in large groupings of neurons, predominantly in the hindbrain, which were associated with movement events through their correlation with motion artifacts in the imaging data (Figure 1-figure Supplement 1).

To explore the spatial patterns in these data we used a 2-dimensional lookup table to visualize tuning with regressors representing either DF stimuli or movement (Figure 1G, H). This revealed segregated populations of neurons coding for the DFs (pink) and movement (green/teal). As expected, DF-tuned neurons were located predominantly in visual sensory areas of the midbrain (tectum) and the diencephalon (pretectum and thalamus). Motor-coding neurons dominated in the hindbrain, with the exception of the cerebellum and inferior olive, which was predominantly tuned to the sensory stimulus. Some neurons did show approximately equal correlation values to both stimuli, as evidenced by the blue-ish hues. Finally, some areas of the brain appeared to contain mixtures of neurons with different coding properties, including the ventral diencephalon and midbrain.

To determine if there was any spatial logic to how different neurons adapt their responsiveness to DFs during imaging, we plotted the ROIs using a lookup table highlighting the preference of for either the first three DFs (pink, naive response), or last three DFs (green, trained response). Strong preferences for the naive stimuli reflects a depressing response profile (Figure 1I, J). While most neurons did show tuning consistent with strong depression, there were neurons that showed an equal preference for naive and trained stimuli, or even stronger preference for the latter, indicating stable or potentiating response profiles. These non-depressing neurons were mostly contained in the dorsal regions of the brain, including the torus longitudinalis, cerebellum and dorsal hindbrain. These results demonstrate that while the majority of neurons across the brain depress their responsiveness during habituation, a smaller population of neurons exists that show the opposite pattern.

Functional classification and anatomical localization of neuronal types observed during habituation learning

To explore the functional heterogeneity within the DF-tuned neurons we used affinity propagation clustering. This method has the advantage that cluster numbers do not need to be defined beforehand, and instead attempts to identify the most representative response profiles (Förster et al., 2020). This identified 12 clusters that differed both in their adaptation to repeated DFs, as well as the shape of their response to the DF (Figure 2A,B).

Figure 2.

We therefore use these two aspects of the response to label the clusters:

Adaptation Profile

No Adaptation = ^noA : Cluster 1, 9, 10

Weak Depression = ^weakD : Cluster 5, 6, 11

Medium Depression = ^{med D} : Cluster 2, 3, 7

Strong Depression = ^strgD : Cluster 4, 8

Potentiation = ^Pot : Cluster 12

Response Shape

On-response = _On : Cluster 1, 2

Long/sustained response = _L : Cluster 3, 4

Medium-length response = _M : Cluster 5, 6, 9

Short/transient response = _S : Cluster 7, 8, 10, 11

Yielding clusters: , and

While these results indicate the presence of a dozen functionally distinct neuron types, such clustering analyses will force categories upon the data irrespective of if such categories actually exist. To determine if our cluster analyses identified genuine neuron types, we analyzed their anatomical localization (Figure 2C-E). Since our clustering was based purely on functional responses, we reasoned that anatomical segregation of these clusters would be consistent with the presence of truly distinct types of neurons. Indeed, we observed considerable heterogeneity both within and across brain regions. For example: was mostly restricted to medial positions within the optic tectum; and were more prevalent within motor-related regions of the brain including the tegmentum and hindbrain rhombomeres; was the most prominent cluster in the torus longitudinalis, consistent with the presence of non-depressing signals in the area (Figure 1I,J).

We then quantified the similarity in the spatial relationships among the clusters by looking at the correlations in the positions of the ROIs in the Z-Brain (Figure 2E). This revealed similar hierarchical relationships to those identified functionally (Figure 2B), especially with respect to Response Shape, indicating that physical location is associated with functional response type.

Finally, since our functional analysis was performed purely based on correlations with the DF stimuli, we asked to what extent neurons belonging to each cluster were correlated with motor output (Figure 2F). This identified as the most strongly correlated to motor output, consistent with its strong habituation profile and its localization within motor-regions of the hindbrain. This indicates that neurons likely occupy the most downstream positions within the sensory-motor network.

These results highlight a diversity of functional neuronal classes active during DF habituation. Whether there are indeed 12 classes of neurons, or if this is an over- or under-estimate, awaits a full molecular characterization. Independent of the precise number of neuronal classes, we proceed under the hypothesis that these clusters define neurons that play distinct roles in the DF response and/or its modulation during habituation learning.

Pharmacological screening to identify habituation modulators

We next used a pharmacological screening approach to both identify molecular mechanisms of habituation and to further probe the habituating circuit. For this we screened 1953 small molecule compounds with known targets (Figure 3-source data 1)), in conjunction with the high-throughput assay we previously established, which has a maximum throughput of 600 larvae/day (Figure 3A, (Randlett et al., 2019)). As we aimed to identify modulators specific for habituation, we included additional behavioural assays as controls, including the response to acoustic stimuli, the optomotor response, and the spontaneous swimming behaviour of the fish in the absence of stimulation (Figure 3B,C). In each 300-well plate, 40 groups of 6 larvae were treated in individual wells, and compared to 60 vehicle treated controls (Figure 3A). We chose these numbers based on a sub-sampling analysis that determined these numbers were sufficient to identify the effect of a known modulator of habituation (haloperidol (Randlett et al., 2019)) at a false-negative rate of less than 0.05 (not shown), while allowing us to screen 80 compounds per experiment across 2 plates.

Figure 3.

We were able to collect the full behavioural record of 1761 compounds (Figure 3D, Figure 3-source data 2)), indicating that the fish survived the treatment and maintained their ability to swim. Behavioural records for fish treated with each compound were compressed into a fingerprint (Rihel et al., 2010) – a vector representing the strictly standardised mean difference (SSMD) across 47 aspects of behaviour (see Methods). For measurements related to dark-flash habituation behaviour, responses were time-averaged across three epochs chosen to highlight changes in habituation: the naive response (first 5 dark flashes), the response during the remaining training flashes, and the re-test block 5 hrs after training (Figure 3B). This was done across 10 different components of the dark flash response (Probability of Response, Latency, Displacement, etc.).

We found that 176 compounds significantly altered at least one aspect of measured behaviour, yielding a 9% hit rate (hit threshold of |SSMD| ≥ 2). While the average effect was to suppress behavioural output , which could reflect non-specific toxicity or a generalized inhibition of motor output, most small molecules induced both positive and negative changes in behavioural output, indicating that toxicity is not the primary phenotypic driver. While the false negative rate is difficult to determine since so little is known about the pharmacology of the system, we note that of the three small molecules we previously established to alter dark flash habituation that were included in the screen, Clozapine, Haloperidol and Pimozide (Randlett et al., 2019), the first two were identified among our hits while Pimozide was lethal at the 10μM screening concentration.

Correlational structure in the pharmaco-behavioural space

To explore the pharmaco-behavioural space in our dataset we clustered the hits based on their behavioural phenotypes (Figure 4A). This strategy can identify compounds that share common pharmacological targets, or perhaps distinct pharmacological targets that result in convergent behavioural effects (Bruni et al., 2016; Rihel et al., 2010). Indeed, compounds known to target the same molecular pathways often showed similar behavioural fingerprints lying proximal on the linkage tree, indicating that our dataset contains sufficient signal-to-noise to recover consistent pharmaco-behaviour relationships.

Figure 4.

Alternatively, compounds can be considered as tools to manipulate different aspects of brain function agnostic to their molecular mechanisms. Consequently, using similarities and differences among the induced alterations should uncover molecular and neural linkages among different behavioural outputs. Following this logic, the ability of a compound to co-modify different aspects of behaviour would reflect molecular and/or circuit-level dependencies. For example, visual behaviours that all depend upon photoreceptors should be similarly affected by any compounds that modulate phototransduction in these photoreceptors. We quantified these relationships by calculating the correlated effects on our different behavioural measurements across compounds (Figure 4B).

Consistent with our previous results highlighting uncorrelated learning across the behavioural components of the O-bend response during habituation (Randlett et al., 2019), we found that different aspects of the response were independently affected pharmacologically, resulting in distinctive correlated groupings within the correlogram. While we previously found that O-Bend response Probability and Latency habituate independently in individual fish (Randlett et al., 2019), in our small molecule screen data these appear to be tightly coupled (Figure 4B). The performance of the animals in the OMR assay under different treatments was also associated with O-bend Probability and Latency, suggesting that pharmacological modulation of vision or arousal could drive these correlations within the small molecule screen dataset.

These analyses confirm habituation behaviour manifests from multiple distinct molecular mechanisms that independently modulate different behavioural outputs.

Modulation of habituation by GABA, Melatonin and Estrogen signaling

For the remainder of the analyses we decided to focus on the mechanisms leading to the habituation of response probability, as this is the criterion for which it is easiest to identify the link between neural activity and behavior, providing the best entry point for studying the circuit mechanisms of long-term habituation. To identify the most promising hits, we sought to identify compounds that:

1. Have minimal effects on the naive response to DFs, but strong effects during the training and/or memory-retention periods. This would prioritize pathways that affect habituation, rather than simply DF responsiveness.

2. Have minimal effects on other aspects of behaviour, in order to exclude compounds that would alter generalized arousal, movement ability/paralysis, or visual impairment. Such compounds would strongly influence DF responsiveness, but likely independently of pathways related to habituation.

3. Show similar behavioural effects to other compounds tested that target the same molecular pathway. Such relationships can be used to cross validate, yet we note that our library choice was very broad, and target coverage is non-uniform. Therefore a lack of multiple hits targeting the same pathway should not be taken as strong evidence of a false positive.

This manual prioritization led to the identification of the GABA_A/C Receptor antagonists Bicuculline, Amoxapine, and Picrotoxinin (PTX). PTX treatment had the strongest effects, with increased responsiveness to DFs during the training and test periods, indicative of defects in habituation (Figure 4Ci). Dose-response experiments confirmed a strong effect of PTX on inhibiting the progressive decrease in responsiveness during habituation learning at 1-10μM doses (Figure 5A). Importantly, like the naive dark-flash response, the probability of responding to an acoustic stimulus and the optomotor response were not inhibited (Figure 5-figure Supplement 1A). While strong GABA_A/CR inhibition results in epileptic activity in larval zebrafish, we did not observe evidence of seizure-like behaviour at these doses, consistent with a partial GABA_A/CR in our experiments and previous results (Bandara et al., 2020). Therefore, we conclude that partial antagonism of GABA_AR and/or GABA_CR is sufficient to strongly suppress habituation but not generalized behavioural excitability, indicating that GABA plays a very prominent role in habituation. This is consistent with a model in which the potentiation of inhibition actively silences sensory-induced activity during habituation to suppress motor output (Cooke and Ramaswami, 2020; Ramaswami, 2014).

Figure 5.

We next turned our attention to the upper portion of the clustered behavioural fingerprint graph (Figure 4A), where strong and relatively specific inhibition of responsiveness during training and testing were observed, indicative of enhanced habituation (Figure 4Cii, iii). Among the hits observed here were multiple agonists of both Melatonin and Estrogen receptors, indicating that hormonal signaling may play a prominent role in habituation. Dose response studies with Melatonin confirmed strong potentiation of habituation (Figure 5B). Melatonin did cause a decrease in spontaneous movement behaviour, consistent with its role in arousal/sleep regulation in zebrafish and other vertebrates (Gandhi et al., 2015; Zhdanova et al., 2001), yet Melatonin did not inhibit the naive response to dark flashes, the responsiveness to acoustic stimuli or OMR performance (Figure 5B, Figure 5-figure Supplement 1B). Melatonin’s effect on habituation was also most prominent for the Probability of response, and did not strongly alter habituation for Displacement Figure 5-figure Supplement 1F, indicating it does not cause generalized sedation but modulates specific aspects of behaviour at these doses, including increasing habituation of the Probability of response.

We similarly validated that the Estrogen Receptor agonists Ethinyl Estradiol and Hexestrol, potentiated habituation at 5-100μM and 1-10μM doses, respectively (Figure 5C,D). Ethinyl Estradiol strongly suppressed movement rates at these doses, and both treatments suppressed acoustic responsiveness and OMR performance at doses >=10μM (Figure 5 - figure Supplement 1C,D). Thus, it is less clear how specific or generalized Estrogen Receptor agonism is on behaviour, although the effective doses of Hexestrol for influencing habituation (1-5uM) were lower than those that significantly affected other behaviours (10uM). Nevertheless we decided to focus on PTX and Melatonin for the remaining experiments.

Our screening approach identified both expected (GABA) and unexpected (Melatonin, Estrogen) pathways that strongly modulate habituation of responsiveness. We also implicated other pathways in habituation, including Progesterone and PPARα (Figure 4C), and identified compounds that strongly modify other aspects of behaviour (OMR, acoustic and spontaneous behaviour). These hits can be mined for future projects investigating the molecular basis of behaviour.

Pharmacological manipulations of functional circuit properties during habituation

Our Ca²⁺ imaging experiments identified 12 distinct functional classes of neurons during habituation learning, but it is unclear how these might be organized in a circuit. Based on the diversity of functional response profiles identified, it is clear that solving this circuit will take considerable further work. As as starting point in this long-term effort, we used the pharmacological manipulations as these treatments provide us with tools to ask how treatments that potently alter habituation behaviour also alter the functional properties of neurons. We compared the Ca²⁺ activity patterns after treatment with vehicle (0.1% DMSO), PTX, or Melatonin (Figure 6). At the behavioural level, we found a trend indicating that we were able to manipulate habituation pharmacologically in our tethered imaging assay, though this was very subtle (Figure 6A). This discrepancy relative to the very strong behavioural effects in freely-swimming animals (Figure 5) likely result from the head-restrained protocol, which itself strongly inhibits behavioural output. Yet, since we did observe a trend in behavioural data, we proceeded under the assumption that the compounds were having the desired effects.

Figure 6.

As PTX and Melatonin have opposing effects on habituation behaviour, we reasoned that these two treatments should have opposite effects in the circuit, with PTX inhibiting depression and Melatonin promoting depression. Indeed Melatonin has been found to increase the effects of GABA, and so such a relationship could be direct (Cheng et al., 2012; Niles et al., 1987). In contrast to this straightforward hypothesis, what we observed was considerably more complex. We did not observe alterations of the average response profiles of individual neuronal classes, which remained indistinguishable after the treatments (Figure 6-figure Supplement 1C-K). Instead, the proportion of neurons that belonged to the different classes was altered (Figure 6B-D). Therefore, the pharmacological manipulations did not alter the activity of neurons in such a way as to alter the average activity states of populations, but instead the proportion of neurons belonging to different populations changed. This may point to fixed and relatively inflexible processing strategies that the brain is using in the context of dark-flash habituation which constrain the possible functional response types.

The effect of PTX on cluster reassignment generally tended towards weaker depression, increasing the proportion of cells falling into the weaker depressing classes at the expense of strongly depressing classes for a given response shape (Figure 6D). This pattern was most clear in the classes with “Short” and “Long” Response Shapes, which are those that included the most strongly depressed classes of neurons.

Based on the hypothesis that Melatonin and GABA cooperate during habituation, we expected PTX and Mela-tonin to have opposite effects. This clearly does not fit with our observations: for example, the size of the neuron population was increased by both PTX and Melatonin (Figure 5C). While habituation of the Probability of response is oppositely modulated by PTX and Melatonin, this is not true of behaviour globally – the behavioural fingerprints of Melatonin and GABA are not opposites (Figure 5A,B) and opposing effects are not seen for the habituation of Displacement (Figure 5-figure Supplement 1E,F). Therefore, a lack of coherent shifts across the entire neural population when applying these treatments is expected. However, opposite effects of PTX and Melatonin were observed for neurons (Figure 6C), which we found to be most strongly correlated with motor output (Figure 2F). Therefore, this class might be most critical for habituation of response Probability.

Combined, these experiments reveal that pharmacological manipulations that affect habituation behaviour manifest in complex functional alterations in the circuit. These effects can not be captured by a simple model, and considerable additional knowledge of the circuit, including the connectivity and signalling capacity of different neurons will be necessary to understand these dynamics.

Identification of GABAergic neurons classes in the habituating circuit

Since our pharmacological experiments point to the importance of GABAergic inhibition in habituation, we asked which functional classes of neurons are GABAergic? An obvious model for habituation would assign a GABAergic identity to the neurons that potentiate their responses, and thus could act to progressively depress the responses of other neuronal classes. We began with virtual co-localization analyses with 3D atlases to identify candidate molecular markers for functionally identified neurons. Such a strategy can be powerful to generate hypotheses from brain-wide imaging data, provided sufficient stereotypy exists in the positioning of neurons, and the relevant marker exists in the atlas (Dunn et al., 2016; Randlett et al., 2015). Therefore, we analyzed the spatial correlations for markers contained in the Z-Brain (Randlett et al., 2015), Zebrafish Brain Browser (Gupta et al., 2018; Marquart et al., 2017; Tabor et al., 2018), and mapZebrain atlases (Kunst et al., 2018; Shainer et al., 2022). We identified markers showing the highest spatial correlations with any of our functional clusters (corr. > 0.15, n=89 of 752 markers), and organized these hierarchically (Figure 7A). GABAergic reporter lines based on the gad1b promoter were located in a region of the hierarchy showing greatest spatial similarity with and (Figure 7B-E). An enrichment along the medial tectum is common to markers in this region of the hierarchy, where the highest density of GABAergic neurons within the tectum reside.

Figure 7.

To confirm that and classes are GABAergic, we imaged the response of neurons in Tg(Gad1b:DsRed); Tg(elavl3:H2B-GCaMP6s) double transgenic larvae, and classified neurons as gad1b-positive or -negative based on DsRed/GCaMP levels (Figure 7F,G). Indeed we saw a heterogeneous distribution of gad1b-positive neurons across functional clusters, including a significant enrichment in not only and , but also the other two clusters with the “Short” Response Shape ( and ). The remaining clusters either showed no significant bias, indicating that they contain mixed populations, or a significant depletion of gad1b-positive cells, suggesting that they comprise mostly of excitatory or neuromodulatory neurons ( and ).

These experiments indicate that GABAergic neurons in the habituating circuit are not characterized by their Adaptation Profile (other than non-potentiating), and instead have a characteristic “Short” Response Shape, perhaps reflecting a transient bursting style of activity relative to other neuronal types that exhibit more sustained firing patterns. This lack of coherence in adaptation profile may explain why global manipulations of GABAergic signaling through PTX have complex manifestations in the functional properties of neurons (Figure 6D).

Discussion

Molecular mechanisms of DF habituation

To explore the molecular mechanisms of habituation, we performed a small molecule screen testing for effects on DF habituation behaviour. Analyses of the correlated effects of drugs across different aspects of behaviour (Figure 4) are consistent with our previous results indicating that habituation results from multiple molecularly independent plasticity processes which act to adapt different aspects of the DF response during habituation (Randlett et al., 2019). Here we focused our analysis on those pharmacological agents and pathways that strongly and relatively specifically modulated habituation when measuring response Probability. We found that inhibition of GABA_A/C Receptors using PTX reduced habituation learning. GABA is the main inhibitory neurotransmitter in the zebrafish brain, and deficits in GABA signaling lead to epileptic phenotypes (Baraban et al., 2005). We were fortunate that our screening concentration (10μM) did not cause seizures, but was still sufficient to inhibit habituation. This implies that the habituation circuit is exquisitely sensitive to changes in GABA signaling at levels well below the threshold required to drastically change excitatory-inhibitory balances. We cannot rule out the possibility that off-targets of PTX, or subtle non-specific changes in excitatory/inhibitory balance alter habituation behaviour. However, the lack of strong modulation of other behaviours, including the response to acoustic stimuli or the optomotor response (Figure 5-figure Supplement 1A), suggests that GABAergic inhibition plays a crucial role in the process of DF habituation.

A critical role for GABA in habituation is also consistent with data from Drosophila, where both olfactory and gustatory habituation have been linked to GABAergic interneurons (Das et al., 2011; Paranjpe et al., 2012; Trisal et al., 2022). Therefore, this circuit motif of increasing inhibition to drive habituation may be a conserved feature of habituation, which would allow for a straightforward mechanism for habituation override during dis-habituation via dis-inhibition (Cooke and Ramaswami, 2020; Trisal et al., 2022).

Our screen also identified that neuro-hormonal signaling is critical for habituation, where Melatonin and Estrogen receptor agonists potently increase habituation learning rate. The role of Estrogens in learning and memory is well established (Luine et al., 1998; Nilsson and Gustafsson, 2002). Though its role in habituation is less well explored, it has previously been shown to increase memory retention for olfactory habituation in mice (Dillon et al., 2013). To our knowledge, Melatonin has not previously been implicated in habituation, though it has been implicated in other learning paradigms (El-Sherif et al., 2003; Jilg et al., 2019). Notably, Melatonin was shown to block operant learning at night in adult zebrafish (Rawashdeh et al., 2007), and therefore Melatonin appears to be able to both promote or inhibit plasticity in zebrafish, depending on the paradigm.

While Melatonin and Estrogen were not strong candidates for involvement in DF habituation plasticity before our screen, their previous associations with learning and memory reinforce the idea that these molecules play critical roles in plasticity processes. In support of this idea, we have previously shown that habituation is regulated in a circadian-dependent manner (Randlett et al., 2019), and both Melatonin and Estrogen levels fluctuate across the circadian cycle (Alvord et al., 2022; Gandhi et al., 2015; Zhdanova et al., 2001), suggesting that either or both of these pathways may act to couple the circadian rhythm with learning performance.

Finally, approximately 2% of the US population use Melatonin as a sleep-aid (Li et al., 2022), and a substantial proportion of US women take Estrogen as part of either oral contraceptives or hormone replacement therapy. Therefore, understanding the roles these molecules play in neuroplasticity is a clear public health concern.

Circuit mechanisms of DF habituation

Based on behavioural experiments, we previously postulated that multiple plasticity loci cooperate in the habituating dark-flash circuit, arranged in both parallel and series within the circuit (Randlett et al., 2019). Here, our Ca²⁺ imaging experiments identified a diverse range of neuronal Adaptation Profiles, including non-adapting and potentiating neurons spread throughout sensory- and motor-related areas of the brain. Thus, non-habituated signals are transmitted throughout the brain, consistent with a distributed learning process. Such a model is further supported with brain-wide imaging data for short-term habituation to looming stimuli, where distributed neurons were identified that showed differential rates of habituation (Marquez-Legorreta et al., 2022). It is important to point out that Marquez-Legorreta et al. did not observe non-adapting or potentiating neurons in their experiments. This may be due to differences in analysis methods, or could highlight a difference between short- and long-term habituation circuit mechanisms, the latter of which may rely on more complex circuit mechanisms involving both potentiation and suppression of neuronal responses.

We also observed classes exhibiting an On-response profile ( and ). These neurons fire at the ramping increase in luminance after the DF, making it unlikely that they play a role in aspects of acute DF behaviour we measured here. These neurons exist in both non-adapting and depressing forms suggesting a yet unidentified role in behavioural adaptation to repeated DFs.

While we have insufficient anatomical data to constrain circuit connectivity models that drive DF habituation, here we demonstrate the use of pharmacology, functional imaging and neurotransmitter classifications to constrain our models. Specifically, pharmacology indicated a role for GABA and Melatonin in habituation, and our functional imaging identified distinct classes of neuronal types in the DF circuit, including potentiating neurons . These results point to a model where neurons are GABAergic and thus progressively inhibit the other neuronal classes, and that perhaps this effect is bolstered by Melatonin. However, in silico co-localization analyses and double transgenic Ca²⁺ imaging identified neurons as predominantly non-GABAergic, inconsistent with this simple model. Instead, we found that the GABAergic neurons in the circuit are characterized by their short burst of activity to the stimulus onset. If the GABAergic neurons are not increasing in their firing rates but do drive habituation, then per-haps it is the potentiation of GABAergic synapses that drives habituation. This is a somewhat unexpected model, as studies of long-term synaptic plasticity (e.g. LTP and LTD) have overwhelmingly focused on plasticity at excitatory synapses. Although a functional link to behaviour is less well established, long-term inhibitory synaptic plasticity has been well documented, including inhibitory (i)-LTP and i-LTD (Castillo et al., 2011). Alternatively, there may be a key minority subset of neurons that are GABAergic and exert a strong influence over the rest of the circuit driving depression and habituation.

We also found that the same pharmacological treatments that result in strong alterations to habituation behaviour in freely swimming larvae (Figure 5), resulted in relatively subtle and complex functional alterations in the circuit (Figure 6). Making direct comparisons between freely-swimming behaviour and head-fixed Ca²⁺ imaging is always challenging due to the differences in behaviour observed in the two contexts, and therefore our failure to identify a clear logic in these experiments may have technical explanations that will require approaches to measure neural activity from unrestrained and freely-behaving animals to resolve (Kim et al., 2017). Alternatively, these results are again consistent with the idea that habituation is a multidimensional and perhaps highly non-linear phenomenon in the circuit, which cannot be captured by a simple model.

Circuit loci of DF habituation

Where in the brain does habituation take place? As discussed above and previously, our data is inconsistent with a single-locus of plasticity (Randlett et al., 2019). Instead, we propose that plasticity is distributed throughout the circuit. Since PTX inhibits most aspects of habituation learning (Figure 5Ai), these all may all involve GABAergic motifs. Moreover, the different functional classes of neurons are distributed through sensory- and motor-related areas of the brain, consistent with the notion that habituation plasticity occurs in a very distributed manner. While distributed, there are clear associations between anatomical location and functional neuron type (Figure 2A-E), indicating that there is some degree of regional logic to the localization of Adaptation Profiles. For example, and are the most prevalent in the pretectum, and mostly absent from the tegmentum and posterior hindbrain, whereas and are numerous in tegmentum and posterior hindbrain, and thus likely occupy more downstream positions in the sensori-motor circuit.

The tectum is one of the largest brain areas in larval zebrafish, and is directly innervated by nearly all retinal ganglion cells (Robles et al., 2014). Therefore, the tectum is a prime candidate for implementing DF habituation for anatomical reasons. In further support of this notion, the neurons we have identified as GABAergic and propose to be driving habituation ( and ) are concentrated in the tectum (Figure 2C,D). The tectum contains multiple anatomically distinct types of GABAergic neurons, most of which are locally projecting interneu-rons (SINs, ITNs, PVINs), although GABAergic projection neurons have been observed with axons projecting to the anterior hindbrain (Gebhardt et al., 2019; Martin et al., 2022; Nevin et al., 2010; Robles et al., 2011). Therefore, we expect that our GABAergic classes correspond to subsets of these GABAergic tectal neurons, which is testable using genetic approaches based on marker co-expression and/or single cell morphometric and transcriptomic analyses.

Beyond the tectum, conspicuous neuronal clustering was observed in the inferior olive and cerebellum, which have been implicated in motor-related learning behaviours in larval zebrafish (Ahrens et al., 2012; Lin et al., 2020; Markov et al., 2021). Both structures contained many stimulus-tuned neurons (Figure 1I), and non-adapting ( and ), and potentiating neurons were among the most concentrated in the cerebellum (Figure 2C,D). Non-adapting neurons were also prominent in the torus longitudinalis. The torus longitudinalis has recently been implicated in the binocular integration of luminance cues (Tesmer et al., 2022), and therefore is ideally placed to influence habituation to whole-field stimuli like DFs.

Collectively, our brain-wide imaging data indicate that the adaptations underling habituation span many regions of the brain, and therefore a comprehensive model will need to span many regions of the brain in order to explain the neural and behavioural dynamics underlying habituation learning

Conclusion

Habituation is the simplest form of learning, yet despite its presumed simplicity a model of how this process is regulated in the vertebrate brain is still emerging. Here we have combined two methods offered by the larval zebrafish model: whole brain functional imaging and high-throughput behavioural screening. By applying these methods to long-term habituation, we identified and validated pharmacological agents that strongly modulate habituation learning, and distinct classes of neurons that are activated by DFs and adapt their activity during learning. The systematic datasets we generated contain large amounts of additional information that await future validation and integration into our understanding of DF habituation. Nonetheless, the diversity of molecular pathways and functional neuronal types we have identified here indicate that considerable biological complexity exists that awaits discovery within the “simplest” form of learning.

Methods

Animals

All experiments were performed on larval zebrafish at 5 days post fertilization (dpf), raised at a density of ≈1 larvae/mL of E3 media in a 14:10h light/dark cycle at 28-29°C. Wild type zebrafish were of the TLF strain (ZDB-GENO-990623-2). Transgenic larvae used were of the following genotypes: Tg(elavl3:H2B-GCaMP7f)^jf90 (Yang et al., 2021), Tg(elavl3:H2B-GCaMP6s)^jf5 (Freeman et al., 2014), and Tg(gad1b:DsRed)^nns26 (Satou et al., 2013). Zebrafish were housed, cared for, and bred at the Harvard MCB, UPenn CDB, and Lyon PRECI zebrafish facilities. All experiments were done in accordance with relevant approval from local ethical committees at Harvard University, the University of Pennsylvania, and the University of Lyon.

High-throughput screening setup and protocol

Larvae were assayed for behaviour in 300-well plates using the apparatus described previously (Randlett et al., 2019). Briefly, each well is 8mm in diameter and 6mm deep, yielding a water volume of ≈300uL. Behaviour plates are suspended below a water bath kept at 31°C, which acts as a heated lid to prevent condensation and maintains the water temperature in the well at 29°C. Behaviour was tracked using a Mikrotron CXP-4 camera, Bitflow CTN-CX4 frame grabber, illuminated with IR LEDs (TSHF5410, digikey.com). Visual stimuli were delivered via a ring of 155 WS2812B RGB LEDs (144LED/M, aliexpress.com). For a dark flash stimulus, the LEDs were turned off for 1s, and then the light intensity was increased linearly to the original brightness over 20s. The optomotor response was induced by illuminating every 8th LED along the top and bottom of the plate, and progressively shifting the illuminated

LED down the strip resulting in an approximately sinusoidal stimulus, 5.5 cm peak to peak, translating at 5.5 cm per second. Direction of motion was switched every 30 s, for a total testing period of 1 hour, and performance was scored as the average change in heading direction towards the direction of motion during these 30s epocs. Acoustic tap stimuli were delivered using a Solenoid (ROB-10391, Sparkfun). The behavioural paradigm was designed to be symmetrical such that 1hr worth of stimulation was followed by 1hr worth of rest (Figure 1B), allowing us to alternate the view of the camera between two plates using 45-degree incidence hot mirrors (43-958, Edmund Optics) mounted on stepper motors (Figure 1A, ROB-09238, Sparkfun), driven by an EasyDriver (ROB-12779, Sparkfun).

Apparatus were controlled using arduino microcontrollers (Teensy 2.0 and 3.2, PJRC) interfaced with custom written software (Multi-Fish-Tracker), available here: github.com/haesemeyer/MultiTracker.

The protocol for assessing behaviour (Figure 1B, Figure 3B) consisted of dark flashes repeated at 1-minute intervals, delivered in 4 training blocks of 60 stimuli, separated by 1hr of rest (from 0:00-8:00, hr:min of the protocol). For analyses, this epoch is separated into periods reflective of the Naive response (first 5 stimuli), and the remaining 235 stimuli during training. From 8:00-8:30, no stimuli are delivered and fish are monitored for spontaneous behaviour. From 8:30-9:00 fish are given acoustic stimuli, and from 10:00 - 11:00 fish are assayed for the optomotor response and turning towards the direction of motion (light blue). Finally, at 12:00-13:00, larvae are given 60 additional dark flashes during the test period (red).

Behavioural analyses

The behaviour of the fish was tracked online at 28 hz, and 1-second long videos at 560 hz were recorded in response to DF and Acoustic Tap stimuli. Offline tracking on recorded videos was performed in MATLAB (Mathworks) using the script “TrackMultiTrackerTiffStacks_ParallelOnFrames.m”, as described previously, to track larval posture (Randlett et al., 2019). Tracks were then analyzed using Python. Analysis code available here: github.com/owenrandlett/lamire_2022.

Responses to DFs and to taps were identified as movement events that had a bend amplitude greater than 3rad and 1rad, respectively. Behavioural fingerprints were created by first calculating the average value for each fish reflecting either the DF response during the specified time period (Naive = DFs 1-5, Training = DFs 6-240, Test = DFs 241-300), or the average response during the entire stimulus period (Acoustic Taps, OMR, Free Swimming). Periods where the tracking data was incomplete were excluded from the analysis. DFs where larvae did not respond were excluded from the behavioural components other than the Probability of Response. The Strictly Standardized Mean Difference was then calculated for each of these average fish values for the compound-treated larvae relative to the vehicle (DMSO) control (Figure 3C). The threshold for determining hit compounds was set at |SSMD| ≥ 2. These analyses were performed using: Analyze_MultiTracker_TwoMeasures.py.

Hierarchical clustering (Figure 3D, Figure 4A-C) was performed using SciPy (Virtanen et al., 2020). Correlations across different behavioural measures (Figure 4B) was calculated computing all pairwise comparisons for each behavioural measure in the SSMD fingerprint across the 176 hit compounds.

Further details and code for the analyses used to create the figure panels are in the following notebook: 2022_LamireEtAl_BehavFigs.ipynb. Analyses made use of open-source Python packages, including: NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), matplotlib (Hunter, 2007), seaborn (Waskom, 2021), and open-cv (Bradski, 2000).

Pharmacology

Compounds were prepared as 1000x frozen stock solutions in DMSO. Stock solutions were initially diluted 1:100 in E3, yielding a 10x solution. 30uL of this solution was then pipetted into the wells, yielding a 1x compound solution in 0.1% DMSO (Sigma). Vehicle treatment followed the same protocol, using pure DMSO. Larvae were incubated in compound solution for between 30 to 90 minutes prior to behavioural testing.

The small molecule compound library (Selleckchem Bioactive: FDA-approved/FDA-like small molecules, Figure 1-source data 1) was obtained from the UPenn High-Throughput Screening Core. The library concentration was 10mM, and thus all compounds were screened at approximately 10μM. For subsequent pharmacological experiments chemicals were obtained from: Picrotoxinin: Sigma, P-8390; Melatonin: Cayman, 14427; Sigma, M5250; Ethinyl Estradiol: Cayman, 10006486; Hexestrol: Sigma, H7753

Microscopy

Imaging was performed on 5dpf larvae, mounted tail-free in 2% LMP agarose (Sigma A9414) in E3, using a 20x 1.0NA water dipping objective (Olympus). Volumetric Imaging (Figure 1,Figure 2, Figure 6) was performed at 930 nm on a Bruker Ultima microscope at the CIQLE imaging platform (Lyon, LYMIC), using a resonant scanner resonant scanner over a rectangular region of 1024×512 pixels (0.6μm x/y resolution) and piezo objective mount for fast z-scanning. Imaging sessions began by taking an “Anatomy Stack” consisting of 150 slices at 1μm z-steps, summed over 128 repeats (imaging time ≈11 minutes). This served as the reference stack used for alignment to the Z-Brain atlas, and to detect Z-drift in the imaging session (see below). The functional stack consisted of 12 slices separated at 10μm steps, thus covering 120μm in the brain acquired at 1.98 hz). To image Tg(elavl3:H2B-GCaMP6s);Tg(gad1b:DsRed) double transgenic larvae (Figure 7), we used a custom built 2-photon microscope (Haesemeyer et al., 2018), imaging 512×512 images at (0.98 μm x/y resolution) at 1.05 hz. The anatomy stack was taken at 2 μm step sizes for both the green and red channels in the dark. Functional imaging was performed only on the green/GCaMP channel since the red stimulus LED was incompatible with DsRed imaging.

When developing this protocol we determined that substantial shifts of more than a cell-body diameter (5uM) in the Z-plane are common during the ≈1.2 hrs of imaging. We determined this by comparing the sum of the functional image planes during 5 equally sized time epochs (1540 frames per epoch), aligned to the “Anatomy Stack”, using “phase_cross_correlation” in the scikit-image library (van der Walt et al., 2014). This allowed us to quantify shifts in the imaging plane as shifts in this alignment. These tended to occur within the first hour of imaging, therefore we performed an hour of imaging of this functional stack before beginning the DF stimulation protocol to allow the preparation to settle under imaging conditions. Dark flashes were delivered using a 3mm red LED mounted above the fish, controlled by an Arduino Nano connected to the microscope GPIO board and the Prairie View software to deliver pulses of darkness consisting of 1 sec light off, 20 sec linear ramp back to light on, delivered at 60 second intervals.

Even with this pre-imaging protocol, z-shifts were still observed in a considerable number of fish. Since our habituation-based analysis is focused on how individual neurons change their responses over time, shifts in the z-plane are extremely problematic as they are not correctable post-acquisition and can result in different neurons being imaged at individual voxels. This could easily be confused for changes in functional responses over time during habituation. Therefore, any fish showing a z-drift of greater than 3μm was excluded from our analysis. Stable z-positioning was further confirmed by manual inspection of the eigan images in the imaging timecourse using “View registration metrics” in suite2 to confirm that these do not reflect z-drift. Of 56 larvae imagined total, 22 were excluded, leaving 34 included. Larvae were treated with 0.1% DMSO, Picrotoxinin (PTX, 10uM), or Melatonin (1uM), from approximately 1hr before imaging. These fish were analyzed as a single population (Figure 1,Figure 2) and separately to determine the effects of the treatments (Figure 6).

To quantify responses to the dark flash stimuli we used motion artifacts in the imaging data to identify frames associated with movements (Figure 1-figure Supplement 1). Motion artifact was quantified using the “corrXY” parameter from suite2p, which reflects the peak of phase correlation comparing each acquired frame and reference image used for motion correction. The “motion power” was quantified as the standard deviation of a 3-frame rolling window, which was smoothed in time using a Savitzky-Golay filter (window length = 15 frames, polyorder = 2). A response to a dark flash was defined as a “motion power” signal greater than 3 (z-score) occurring within 10-seconds of the dark-flash onset, and was used to quantify habituation in the head-embedded preparation (Figure 6A).

Ca²⁺ imaging analysis

ROIs were identified using suite2p (Pachitariu et al., 2017) using the parameters outlined in RunSuite2p_BrukerData_ScreenPaper.py and RunSuite2p_MartinPhotonData_ScreenPaper.py scripts for the data from the Bruker Ultima microscope (Figure 5-Figure 7), and custom built 2-photon microscope (Figure 7F,G), respectively. These ROIs mostly reflected individual neuronal nuclei/soma. The clustered heatmap image of neural activity ((Figure 3F) was generated using the suite2p GUI using the “Visualize selected cells” function, and sorting the neurons using the rastermap algorithm (Pachitariu et al., 2017), github.com/MouseLand/rastermap). The imaging planes were then aligned to the anatomical stack taken before functional imaging using “phase_cross_correlation” in the scikit-image library (van der Walt et al., 2014). For the volumetric data, the anatomical stack was then aligned to the Z-Brain atlas coordinates using CMTK, and ROI coordinates were transformed into Z-Brain coordinates using streamxform in CMTK. These steps were performed using Bruker2p_AnalyzePlanesAndRegister.py.

To identify ROIs that were correlated with the stimulus we use a regression-based approach (Miri et al., 2011), where we identified ROIs that were correlated with vectors representing the time course of the DF stimuli convolved with a kernel approximating the slowed H2B-GCaMP time course with respect to neuronal activity. These regressors reflected either the entire 21 second dark flash stimulus, or only the onset of the flash, and either the first 3, last 3, or all 60 flashes (6 regressors in total). To identify neurons correlated to motor output, we took advantage of the plane-based registration statistics calculated by suite2p. Specifically, the “ops[‘corrXY’]” metric, which reflects the correlation of each registered image frame with the reference image. We reasoned that movements would cause image artifacts and distortions that would be reflected as a transient drop in these correlations. Indeed, we confirmed this association by imaging the tail using an infrared camera, and compared the motion index calculated through tail tracking, and that which we calculated based on the motion artifacts, which showed good overall agreement in predicted movement events and average correlation of 0.4, demonstrating that these image-based artifacts can be used as reliable proxies of tail movements (Figure 1-figure Supplement 1). Therefore, regressors based on these motion indices were used to identify neurons correlated with motor output.

Images for the functional tuning of individual neurons (Figure 1G-J) were computed using the the Hue Saturation Value (HSV) colorscheme, with the maximal correlation value to either regressor mapped to saturation, and the hue value reflecting the linear preference for either regressor. Clustering of functional response types (Figure 2) was done by first selecting all those ROIs that showed a correlation of 0.25 or greater with any of the 6 stimulus regressors across all imaged fish. Then among these ROIs we removed any ROIs that did not show a correlation of 0.3 or greater with at least 5 ROIs imaged in a different larvae. This filtered out ROIs that were unique in any individual fish, allowing us to focus on those neuron types that were most consistent across individuals. We then used the Affinity Propagation clustering from scikit-learn (Pedregosa et al., 2011), with “affinity” computed as the Pearson product-moment correlation coefficients (corrcoef in NumPy (Harris et al., 2020)), preference=-9, and damping=0.9, and clustered using Hierarchical clustering (cluster.hierarchy in SciPy (Virtanen et al., 2020)). Cluster number was assigned based on the ordering of the hierarchical clustering tree.

To generate the final cluster assignments we re-scanned all the ROIs calculating their correlation with the mean-response vectors for each of the identified 12 functional clusters, selecting those with a correlation value of 0.3 or greater, which were then assigned to the cluster with which they had the highest correlation. To determine the cluster assignments for the data from Tg(Gad1b:DsRed);Tg(elavl3:H2B-GCaMP6s) double transgenic larvae (Figure 7F,G) data were realigned and interpolated to match the frame rate of the clustered data, and assigned to the 12 clusters as above.

To compare the spatial relationships between the neuronal positions of different functional clusters (Figure 2E), and between the functional clusters and reference brain labels (Figure 7A-E), image volumes were cropped to the imaged coordinates (Figure 1E), downsampled to isometric 10 um³ voxels, and linearized to calculate the Pearson’s correlation coefficient between the image sub-volumes.

Analyses made use of multiple open-source Python packages, including: suite2p (Pachitariu et al., 2017) NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), scikit-learn (Pedregosa et al., 2011), scikit-image (van der Walt et al., 2014), numba (Lam et al., 2015), matplotlib (Hunter, 2007), seaborn (Waskom, 2021), and open-cv (Bradski, 2000). Details of the analyses used to create the figure panels are in the following notebook: 2022_LamireEtAl_FunctionalFigs.ipynb

Data and Code Availability

Code for data analysis and for generating the figure panels is available here: github.com/owenrandlett/lamire_2022 Data are available here: doi.org/10.5061/dryad.jdfn2z3fc

Acknowledgements

We thank the Randlett, Granato and Engert group members for helpful advice regarding the manuscript and work, Gregory Forkin for help in zebrafish husbandry, the Burgess and Baier groups for sharing anatomical data from the Zebrafish Brain Browser and mapZebrain atlases, and Armin Bahl for sharing the mapZebrain to Z-Brain transformation. This work was supported by funding from the ATIP-Avenir program of the CNRS and Inserm (O.R.), a Fondation Fyssen research grant (O.R.), the IDEX-Impulsion initiative of the University of Lyon (O.R.), the NIH grants MH109498 (M.G.), NS123887 (M.H.), U19NS104653 and 1R01NS124017 (F.E.), the National Science Foundation grant IIS-1912293 (F.E.), and the Simons Foundation grant SCGB 542973 (F.E.)

Author Contributions

Conceptualization: OR, MG

Methodology: OR

Software: OR, MH

Formal Analysis: OR

Investigation: L-AL, OR

Supervision, Resources and Funding Acquisition: FE, MG, OR

Writing – Original Draft Preparation: OR

Writing – Review & Editing: MH, FE, MG, LA-L, OR