Introduction

Fear is an emotional state elicited by perceived threat. Behaviorally, fear manifests as stereotypical responses that have a striking similarity across a range of taxa (Adolphs, 2013; Darwin, 1872). These responses can be classified into two broad categories: active (e.g., fighting or fleeing) and passive (e.g., freezing or immobility) with the decision of the right behavior at the right time being critical to survival. For example, evasive behaviors may be ideal when a predator is actively striking whereas freezing (i.e., a lack of movement) may be best when a predator is lurking. Recalling where fearful stimuli have been experienced is key to avoiding predation, enabling animals to anticipate danger and optimize survival strategies (Mobbs et al., 2015). This place-fear learning, termed contextual fear conditioning, has been widely studied, however, the extent to which individuals categorically differ in fear expression and how such differences are represented in the brain, remains poorly understood despite being relevant to a variety of psychiatric disorders like PTSD and anxiety (Mancini et al., 2025; Maren et al., 2013; Ressler et al., 2022).

One of the most reliable and ethologically relevant elicitors of fear in animals are alarm pheromones. These are small molecules released by animals as a form of communication, signaling the presence of a threat (Tirindelli et al., 2009). The existence of alarm pheromones is well established in many species of insects, fish, and mammals (Verheggen et al., 2010) and is being increasingly recognized in humans (de Groot and Smeets, 2017). In fish, the first description of alarm pheromone was by Karl von Frisch (von Frisch, 1938) in minnows (Phoxinus laevis), where he described how an extract of epithelial cells resulted in fear-related behaviors like freezing and erratic swimming. Responses were also species specific; fish preferentially respond to conspecific alarm substance (CAS) while ignoring substances from heterospecifics. He also found that fish formed memories of where they previously experienced exposure to CAS, and that individuals sometimes varied in their behavioral responses. More recently, variability in responsiveness to CAS has been reported in a range of fish species (Pfeiffer, 1977) including zebrafish (Danio rerio) (Speedie and Gerlai, 2008; Suboski et al., 1990). These findings suggest the presence of consistent individual behavioral differences in fear reactivity, but they have not been systematically examined.

Given the ethological relevance of CAS as a stimulus, it holds considerable potential for use in understanding individual variation in associative fear learning in zebrafish. In fear conditioning, a conditioned stimulus (CS) is paired with an intrinsically aversive unconditioned stimulus (US). Later presentation of the CS alone results in a conditioned response that is a measure of the strength of the memory. The CS can be either a discrete cue, like a light or sound, or a diffuse stimulus, like a tank or context. CAS has been successfully used as a US in Pavlovian fear learning assays in zebrafish where it has been associated with lights, odors, and contexts (Baker and Wong, 2019; Hall and Suboski, 1995; Maximino et al., 2018; Ruhl et al., 2017). These data suggest that fear learning using CAS is a robust phenomenon in zebrafish, making it ideal for determining if individual animals vary in fear memory expression and identifying the corresponding neural substrates.

To identify the presence of individual differences in fear behavior, we collected data from over 300 adult zebrafish across four inbred strains and both sexes. To increase throughput, we automated behavioral data analysis by combining DeepLabCut (Mathis et al., 2018) with a random forest supervised machine learning model (James et al., 2013). We then used unsupervised machine learning to identify the presence of individual differences in fear memory associated behaviors. Lastly, to identify the neural substrates of fear memory and variation in its expression, we mapped the whole-brain activity of animals during fear memory recall.

Results

Automated behavioral identification

To assess contextual fear memory in zebrafish we designed a three-day behavioral paradigm (Figure 1A). On the first day, animals were given five minutes to acclimate to the conditioning tank. On the second day, animals were placed in the tank for 10 minutes. Baseline behavior was captured during the first five minutes (i.e., pre-exposure). CAS (derived from a mix of donors from all four strains and both sexes) or vehicle was then added to the tank for the last five minutes. On the third day, fish were placed back in the same tank in the absence of CAS to assess fear memory. To automate the identification of different swimming behaviors, we trained a DeepLabCut (Mathis et al., 2018) network to track three points on the fish: the head, trunk, and tail (Figure 1B). We then used the output from this tracking to train a random forest model to identify five different behaviors: straight swimming, normal turns, freezing, burst swimming, and erratic movements. We used 34 parameters derived from the tracking to capture the movement of the animals over time (Table S1-1). These features were calculated using a sliding window of 45 frames (750 ms) to capture the temporal dynamics of movement. Approximately 14,000 frames were labelled for each of five different behaviors, with a 80/20 split for training/testing. We optimized the number of trees and number of features per tree, finding that 2,000 trees and 10 variables per tree optimized the out of bag error rate (2.88%). The resulting model had at least 94% accuracy for each behavior (Figure 1C). Notably, most of the misclassifications occurred when distinguishing between straight swimming versus normal turns and burst swimming versus erratic movements. Because these labels capture the same behavioral concept, and to simplify the subsequent analysis, we use “normal swimming” to refer to either straight swimming or normal turns and “evasive behavior” to capture erratic movement or burst swimming.

Figure 1.

Contextual fear learning

Given prior reports of behavioral variation in response to CAS, we hypothesized that this variation may be due, in part, to the background strain and/or sex of the animals. When looking at how fish respond to CAS itself, we found that almost all strains and sexes had an increase in freezing compared to vehicle treated animals (Figure 2A, top), and all groups were elevated relative to pre-exposure (Figure 2B, top). The only exception to this was female TL fish, where the vehicle treated animals also increased their freezing behavior, likely as a response to the small disturbance that occurs upon the addition of vehicle to the tank. We analyzed evasive behaviors as percent of time active (i.e., not freezing) since many fish spend a large proportion of their time freezing. This gives us a better estimate of evasion that is less influenced by the more dominant freezing behavior. Evasive behaviors in response to CAS had more variation.

Figure 2.

Compared to vehicle treated animals, CAS treated males in two out of four strains increased evasive behavior (except WIKs and TUs that approached significance; Figure 2A, bottom), whereas there were no changes in females. Compared to pre-exposure, CAS resulted in an increase in evasive behaviors in all groups except male WIKs (Figure 2B, bottom).

We next asked whether fish formed a long-term associative memory between the tank and CAS. To test this, fish were placed back into the tank one day after CAS exposure with no CAS present (i.e., memory day). On memory day, there was elevated freezing compared to vehicle treated animals and relative to pre-exposure in nearly all strains and sexes, with the exception of female TU fish that did not differ from their vehicle treated counterparts and male WIKs and ABs that only approached significance (Figures 2A and 2B, top). As before, evasive behaviors were more variable on memory day. When compared to vehicle treated animals, evasive behavior was only elevated in female TUs (Figure 2A, bottom). Compared to pre-exposure, no strains/sexes had an increase in evasive behaviors (Figure 2B, bottom).

To better understand the dynamics of freezing and evasive behavior, we examined temporal profiles of the behavior during the different stages of the paradigm (Figure 2C). During pre-exposure, there was an initial increase in freezing behavior that quickly decreased to ∼20%.

Evasive behavior remained consistently low throughout the trial. During exposure, CAS-treated animals had an initial burst in evasive behavior that rapidly gave way to persistent freezing. On memory day, we saw an initial increase in freezing behavior in both vehicle and CAS-treated animals, but while the freezing in vehicle treated animals quickly declined, the freezing in CAS treated animals was higher and persisted for longer (Figure 2C). Evasive behavior during memory day was largely similar between vehicle and CAS-treated animals, with only a slight increase in CAS treated animals at the beginning of the trial.

Transitions between behavioral states

The aggregated data we have examined so far does not yield insight into how fish switch between behavioral states. To examine these state changes, we constructed transition matrices binning behavior into 750 ms bins that matched the sliding window used for behavioral identification (Figure 3). We found that at all stages, freezing is the most stable behavior (i.e., the behavior is most likely to follow itself), even during pre-exposure when it is a relatively small proportion of the overall behavior (Figure 3A). Amongst the normal swimming behaviors of straight swimming and normal turns, transitions are most likely between these states across all stages of behavior. A similar pattern is seen for evasive swimming (i.e., erratic movements and burst swimming); however, these behaviors are less likely to occur after normal swimming than normal swimming is to occur after evasive behavior. Overall, there are no large changes in the patterns of transitions across the different stages of the behavioral paradigm. These data suggest that exposure to CAS does not result in wholesale reorganization of behavior, but instead a change in how much of each behavior is selected.

Figure 3.

To capture more subtle behavioral changes, we examined state transitions in response to CAS and at memory day in more detail (Figure 3B). To increase clarity, we collapsed straight and normal turns into “normal swimming” and burst swimming and erratic movements into “evasive behavior”. Here, we see that during both CAS and vehicle treatment, there is a slight destabilization of normal behavior such that normal behavior becomes less likely following itself; instead, in the vehicle group normal behavior is more likely to be followed by evasive behavior, and in CAS treated animals, normal behavior is more likely to be followed by freezing (Figure 3B, top row). This increase in evasion in the vehicle group likely represents a response to the water disturbance that occurs when solution is added to the tank. We also see that CAS causes freezing behavior to become more stable (Figure 3B, right column): compared to pre-exposure, freezing is more likely to occur following normal, evasive, and freezing at both exposure and memory day. Evasive behavior increases in likelihood following itself during CAS treatment, with little to no change following normal or freezing behavior (Figure 3B, middle column). This suggests that CAS results in longer bouts of evasive behavior without increasing its likelihood of occurring. In contrast, the CAS-induced elevation in freezing behavior during exposure and memory appears to be due to an increase in both the duration and number of freezing bouts.

Individual differences in fear memory expression

Given the large amount of variation in behavioral responses to CAS, especially during memory day (Figure 2B), we hypothesized the presence of distinct fear memory behavioral types. To uncover these types, we applied Louvain clustering to a series of k-nearest neighbor networks and identified the k-value that yielded the most robust and best clustering (k = 74; Figure S4-1) (Rajput et al., 2022). Distances between nodes (i.e., individual fish) in the network were calculated using percent evasion, percent evasion of active swimming, and percent freezing. Percent evasion of active swimming was included because fish that exhibit high levels of freezing will necessarily be low on non-freezing behaviors. By calculating evasion as a percent of active (i.e., non-freezing) behavior, we better capture what fish choose to do when not freezing. From this analysis, we identified four behavioral clusters (Figure 4A). We named these clusters (Figure 4B): (1) “non-reactive”: fish low in evasion and freezing; (2) “evaders”: fish low in freezing and high in overall levels of evasive behavior; (3) “evading freezers”: fish high in freezing and evasive behaviors; (4) “freezers”: fish high in freezing and low in evasive behaviors.

Figure 4.

Temporal profiles of behavior also reveal distinct differences (Figure 4C). For example, both freezer groups overlapped in the percent freezing behavior over time but differ in their evasive behavior, where the evading freezers have both higher overall levels of evasion and a greater proportion of active behavior spent evading throughout the entire memory trial. Likewise, the temporal profile of freezing for the evaders and non-reactive fish are equally low, but the evaders slowly increase their evasive behaviors over time whereas non-reactive fish remain low in evasion throughout the trial.

Finally, we asked whether behavioral types are affected by the strain and/or sex of the fish (Figure 4D). We found that in AB and WIK fish there was neither over nor underrepresentation in any of the four behavioral types. In contrast, TU fish were overrepresented in the non-reactive group, and TL males were overrepresented in the evading freezer group. This suggests that TU fish may be a poor choice for studying CAS-induced fear learning.

Behavioral clusters over time

We next asked if behavior during exposure to CAS is predictive of behavior at memory day. To do this, we identified which cluster each fish fell into at each behavioral stage (Figure 5A). In the vehicle treated group, animals were largely identified as either non-reactive or evaders at all three time points. In CAS-treated fish, only 14% of animals belonged to one of the freezing groups at pre-exposure. During CAS exposure, this proportion increased to 66%. On memory day, 48% of CAS treated fish were in one of the two freezing groups, about equally split between evading freezers (28%) and freezers (20%). We also found that the behavioral cluster of the fish on exposure day roughly predicted behavior on memory day (Figure 5B, diagonal).

Figure 5.

Most fish in the non-reactive and evader groups on memory day fell into the same category during exposure (54-56%) or transitioned between evaders and the non-reactive groups (29%). Evaders and non-reactive fish rarely transitioned to freezing groups on memory day (11-16%). In contrast, a substantial minority of animals in the freezing groups transitioned to non-freezing groups (30-41%) with the majority remaining in one of the high freezing clusters (58-69%).

Taken together, these data suggest that fish behavior is largely consistent between exposure and memory days, although there remains a sizable minority of fish that appear to forget and switch from freezing to non-freezing behaviors at memory day.

Brain-behavior covariation

To determine the brain regions underlying different behaviors, we performed whole-brain activity mapping followed by partial least squares (PLS) analysis. We chose the AB strain of fish because they showed approximately equal distribution of the different behavioral types (Figure 4D). For brain mapping, animals were euthanized 15 minutes after fear memory recall and processed for tissue clearing using iDISCO. Neural activity was captured using in situ hybridization chain reaction (HCR) against cfos mRNA. Brains were imaged using light-sheet microscopy and then registered to the adult zebrafish brain atlas (AZBA) (Kenney et al., 2021; Rajput et al., 2025). We used PLS analysis (Krishnan et al., 2011) to determine how brain activity covaries with behavior during fear memory recall (Figure 6). This approach uncovers latent vectors that capture the maximal covariance between neural and behavioral states. The first behavioral contrast identified was significant (P = 0.0090), capturing freezing in the positive direction and evasion/normal swimming in the negative (Figure 6A, left); the second contrast approached, but did not reach significance (P = 0.107; Figure 6A, right). However, because this second contrast uniquely captured evasive behavior, we analyzed it with the caveat that these findings are not as strongly supported as those for the first contrast.

Figure 6.

To determine which brain regions consistently contributed to the first behavioral contrast, we calculated bootstrap ratios (BSRs, i.e., salience scores). BSRs capture a combination of the magnitude and stability of regional contributions to a contrast (McIntosh and Lobaugh, 2004).

For the first contrast (Figures 6B and C, S6-1A), we found numerous positive contributions [BSRs > 2.5; i.e., neural activity covaried with freezing. Although PLS is multivariate, we have provided example univariate relationships between cfos activity and behavioral contrasts to give an intuition as to what relationships are being captured (Figure S6-2)]. This activity pattern suggests that freezing behavior, despite the lack of movement, is not passive at the neural level but instead involves a coordinated brain-wide activity. The cerebellum in particular stands out with several regions with very high BSRs: the CC (cerebellar crest; BSR = 4.7), CCe-m (molecular layer of the corpus cerebelli; BSR = 4.0), MON (medial octavolateralis nucleus; BSR = 3.6), and LCa (caudal lobe of the cerebellum; BSR = 6.5). The MON is one of the primary sites of lateral line input to the brain (Valera et al., 2021; Wullimann and Grothe, 2013) and interacts with the rest of the cerebellum (Dohaku et al., 2019; Pose-Méndez et al., 2023). Thus, this data suggests that freezing involves constant monitoring for movement in the water, possibly to maintain posture and minimize drifting. Maintenance of freezing is likely driven by regions of the reticular formation, which projects to the spinal musculature (Chong and Drapeau, 2007; Huang et al., 2013), where we also see strong associations with contrast 1: the SRF (superior reticular formation; BSR = 5.0), IMRF (intermediate reticular formation; BSR = 4.0), and SRN (superior reticular nucleus; BSR = 3.4). Finally, several regions of the telencephalon positively covaried with contrast 1, such as the nuclei of the lateral olfactory tract (nLOTa: BSR = 3.8, nLOTp: BSR = 3.0), parts of the bed nucleus of the stria terminalis (BSTpd: BSR = 3.8, BSTa: BSR = 3.22), the central and lateral portions of the ventral telencephalon (Vc: BSR = 3.1; Vl: BSR = 3.2), and the posterior and central portions of the dorsal telencephalon (Dp: BSR = 3.0; Dc: BSR = 2.8). These regions are likely working to scan and evaluate external chemosensory threats. The nLOTs and Dp regions receive input from the olfactory bulbs (von Bartheld et al., 1984; Yáñez et al., 2022). These regions likely work in concert with the BSTs, which are a part of the extended amygdala, to assign valence to chemosensory stimuli (Porter and Mueller, 2020). The BSTs likely interact with the Vc, which is proposed to be a part of the central amygdala (Porter and Mueller, 2020), a region important for the expression of emotional learning like fear conditioning (Ciocchi et al., 2010; Wilensky et al., 2006). Finally, the Dc is thought to be homologous to mammalian iscortex (Mueller, 2022; Mueller et al., 2011), integrating input from other parts of the brain, including the telencephalon, while projecting to midbrain regions to organize brain-wide activity (Yáñez et al., 2022).

For the second contrast (Figures 6D, 6E and S6-1B), there were only two regions with positive BSRs, and thus covaried with evasive behavior: the nucleus isthmi (NI: BSR = 2.6) and the subcomissural organ (SC: BSR = 3.1). The NI is important for visual attention and modulating retinotectal input, particularly in response to looming stimuli that signal danger (Fernandes et al., 2021; Gruberg et al., 2006; Henriques et al., 2019). Thus, fish with more evasive behavior may be driven by greater attentiveness to potential visual threats. The SCO has recently been found to secrete glycoproteins that form Reissner fiber’s in the spinal cord that help detect spinal curvature (Cantaut-Belarif et al., 2018; Orts-Del’Immagine et al., 2020), and so this region may be more active in order to repair Reissner fibers that come under stress due to the rapid changes in spinal curvature that occur during evasive movements. Regions related to foraging behavior and mechanosensory feedback negatively covaried with contrast 2 and thus are more engaged during normal swimming. For example, with respect to foraging, we have negative BSRs for regions important for spatial navigation like the periventricular gray zone of the optic tectum (PGZ: BSR = -2.8) and the anterior thalamic nucleus (A: BSR = -2.5), which are involved in processing a wide range of visual information (Baier and Scott, 2024; Baier and Wullimann, 2021; Cheng et al., 2017) and the dorsal tegmental area (DTN: BSR = - 2.7), which is critical for monitoring head direction (Petrucco et al., 2023). The central nucleus of the inferior lobe (CIL: BSR = -3.0) is part of a region implicated in feeding behavior (Muto et al., 2017; Wee et al., 2019), but may also play a role in integrating sensory input (Bloch et al., 2019). Finally, several regions involved in sensory feedback are shut down during evasive behavior, likely in anticipation of the large inputs that would occur during evasive movements. These regions include those important for auditory processing in the nucleus of the lateral lemniscus (NLL: BSR = -3.1) and central nucleus of the semicircular torus (TSc: BSR = -3.0) (Yáñez et al., 2024), and vestibular information like the nucleus lateralis valvulae (NLV: BSR = - 3.6), and the caudal lobe of the cerebellum (LCa: BSR = -2.9) (Pose-Méndez et al., 2023).

Brain networks underlying individual differences in behavior

The brain is best thought of as a network of interacting regions whose complex interactions give rise to behavior. To determine the patterns of activity that are present during different types of fear memory expression, we generated networks from correlated cfos activity between brain regions for each experimental group (Figure 7). The overall pattern of correlations (Figures 7A & B) is similar in the groups exposed to the behavioral tank (i.e., non-reactive, evading freezer, and freezer groups; similarity of matrices: r’s=0.41-0.43) while clearly differing from the habituation group that was not exposed to the behavioral tank (r’s=0.18-0.24). To identify which regional interconnections were significant for network generation, we applied a threshold (FDR < 0.001; Figure 7C). This identified fewer significant correlations in the non-reactive group (306) than the evading freezer (588) and freezer groups (655), suggesting these latter two groups form a more coherent brain network.

Figure 7.

To generate insight into how the patterns of brain activity in the evading freezer and freezer groups differed, we generated functional networks (Figure 8). In these networks, the nodes represent brain regions and edges represent the presence of a suprathreshold correlation (Figure 7C). To better visualize how the networks in these two behavioral groups differed, we identified edges whose confidence intervals did not overlap with those of the other behavioral network or the non-reactive network (Figure 8A, solid lines). Significant edges/correlations that had overlapping confidence intervals are depicted as dashed lines.

Figure 8.

All networks (freezer, evading freezer, non-reactive, and habituation; Figures 8A and S8-1, S8-2) have strong interconnectivity in the telencephalon (blue), including high degree centrality both the lateral zone of the dorsal telencephalon [Dl (degree: evading freezers = 30, freezers = 27, non-reactive = 25, and vehicle = 23)] and the medial zone [Dm (degree: evading freezers = 23, freezers = 27, non-reactive = 20, and vehicle = 27)]. These regions are notable for being the putative homologs to the mammalian hippocampus (Dl) basolateral amygdala (Dm) (Anneser et al., 2024; Ganz et al., 2012; Mueller et al., 2011; Porter and Mueller, 2020), both of which are well known to be critical for contextual fear conditioning in mammals (Tovote et al., 2015). This finding is particularly interesting because it suggests that these regions are activated somewhat non-specifically simply by being placed back into the conditioning tanks.

High inter-connectivity of these regions does not occur following an initial exposure to a novel tank (Rajput et al., 2025), and thus this pattern of activity may represent a general response to being placed back into a previously experienced context.

The evading freezer network is notable for its numerous interactions between the subpallium (ventral telencephalon) and hypothalamus (orange). In particular, there is high degree centrality in the mammalian homologs of the central and medial amygdala (Vd-vd, Vc, and Vp), striatum (Vd), and septum (Vv) (Porter and Mueller, 2020; Rink and Wullimann, 2004). Thus, in evading freezers it appears that the Vd-vd, Vc, and Vp may help evaluate the threat while vacillation between active and passive responses is likely regulated by the septal homolog (Vv; Besnard et al., 2019; Rizzi-Wise and Wang, 2021) with the evasive behavior program likely initiated by the striatal homolog (Vd). Notably, the Vv, is less engaged in the freezer network (degree = 9) where there is less switching between freezing and evasive behaviors. The Vv also has extensive projections to the preoptic area and hypothalamus (Rink and Wullimann, 2004; Rizzi-Wise and Wang, 2021), which contains several regions overrepresented in the evading freezer group [i.e., the dorsal and ventral hypothalamus (Hd and Hv) and the anterior and posterior preoptic areas (PPp and PPa)]. These preoptic and hypothalamic regions are key regulators of autonomic functions homologous to their mammalian counterparts (Biran et al., 2015). This suggests that the evading freezer behavioral type may be driven by an elevated stress response relative to freezers.

In contrast, the freezer functional network is notable for interactions between the telencephalon and the cerebellum, particularly the CC which has several unique connections with the dorsal telencephalon and has much higher degree than in the other two networks (degree: freezers = 25, evading freezers = 2, non-reactive = 0). The other notable aspect of the freezer network is the cluster that contains several unique connections between the pretectum and posterior tuberculum. Several parts of the posterior tuberculum stand out as being highly interconnected in freezers, like the medial and lateral preglomerular nuclei (PGm and PGl, respectively), and the posterior thalamic nucleus (P). These nuclei are thought to serve as key sensory processing regions for vision, gustation, and olfaction that then project to the pallium, not dissimilar in function to the mammalian thalamus (Bloch et al., 2020; Mueller, 2012; Northcutt, 2006). Thus, in freezers, enhanced sensitivity to sensory input may be modulating the pallium to signal to the cerebellum to control and maintain the freezing posture and suppress evasive behavior.

Discussion

The present study is the first to unambiguously identify the presence of individual differences in fear response and memory expression in zebrafish by building on prior work that found fish exhibit a mix of active and passive responses to CAS (Hall and Suboski, 1995; Maximino et al., 2018; Speedie and Gerlai, 2008; von Frisch, 1938). We can now add that individuals express fear memories in different ways: some individuals have high levels of freezing interspersed with normal swimming whereas others intersperse freezing with evasive behavior. Finally, we uncovered patterns of brain activity that differentiate these behavioral types, finding that engagement of some regions in the telencephalon are common across behavioral types, but differentiated by regions like the cerebellum, preglomerular nuclei, preoptic area, hypothalamus, and pretectum.

We found that fish respond to CAS with an increase in both freezing and evasive behaviors, with freezing being the more robust and persistent response that is less influenced by strain and sex. These responses to CAS itself are consistent with previous reports of active and passive responses in fish (Hall and Suboski, 1995; Maximino et al., 2018; Speedie and Gerlai, 2008; von Frisch, 1938). We also found that strain had a modest effect on behavior, with the most prominent finding being that TU fish were overrepresented in the non-reactive cluster during fear memory expression. This echoes prior work where the TU strain was found to have a faster rate of contextual fear memory extinction than the AB and TL strains (Kenney et al., 2017). This could potentially be due to a reduction in learning as the TU fish did not increase their freezing response to CAS as much as the other strains (Figure 2A). We also saw a sex difference with respect to evasive behaviors, where males are more likely to increase their evasive behaviors than females when CAS is present. Freezing, however, is largely similar across sexes. This is consistent with a recent study that found both male and female fish increase freezing in response to CAS (Akinrinade et al., 2023). However, the prior study found no increase in evasive behaviors in response to CAS, and no sex differences. We suspect that this may be because in the present study we assessed evasive behavior as a percentage of time not freezing whereas Akinrinade et al (2023) examined levels of evasive behavior as a proportion of total time. We took our approach because many fish spend more than half the time freezing in response to CAS, and so the opportunity for other behaviors is limited. By examining evasive behavior as the percentage of time spent not freezing (i.e., active time), we likely have greater sensitivity to uncover subtle behavioral differences.

A key finding from the present study was the identification of individual differences in how fish behave during fear memory recall. We found that fish behavior falls into one of four groups: non-reactive, evaders, evading freezers, and freezers. The evading freezer and freezer groups are the most clearly associated with fear memory. This is due to the substantial number of animals that fell into the evader category in the vehicle treated fish and the CAS treated fish during pre-exposure (Figure 5A). This suggests that an increase in evasion in the absence of freezing may be a behavioral response to netting and placement in the tank rather than to the alarm substance (Shishis et al., 2023). We suspect this elevated evasive behavior may also be due to the small size of the tank used in this study, which we chose to limit the amount of CAS needed to elicit a fear response.

Our findings on individual differences in fear memory expression are distinct from prior work examining such differences in rodents and fish that classify animals as being either high or low in fear expression. For example, using an active avoidance paradigm with CAS, Maximino and colleagues (2018) classified zebrafish as high or low in avoidance during memory expression, finding that high avoiders had higher levels of freezing and lower levels of erratic movements during CAS exposure. A similar approach has been taken in rats, where individual differences have been identified by differentiating animals based solely on their freezing responses to conditioned stimuli (Graham and Richardson, 2016; Walters et al., 2016).

Interestingly, rats have been found to exhibit a mix of active (darting) and passive (freezing) fear behaviors during fear conditioning to an electric shock where females exhibit more darting than males (Gruene et al., 2015; Mitchell et al., 2022). Elevated darting was also associated with less freezing in rats. These findings mirror some of the findings in the present study, such as the presence of both active (evasive) and passive (freezing) fear behaviors. However, we do not find as clear of a sex difference in active and passive behaviors during fear memory expression, and we find that an increase in evasive behavior does not necessarily lead to less freezing. These differences may be due to the use of a different unconditioned stimulus (CAS versus electric shock) and/or species.

We find that the cerebellum stands out as being strongly associated with freezing behavior. This conclusion is supported by both the PLS and network analysis. PLS analysis found that neural activity of several cerebellar regions (CC, CCe-m, LCa, and MON) strongly contributed to the covariation of brain activity with freezing behavior. The network analysis also identified the CC as a hub with many unique interactions with the pallium and subpallium.

Relatedly, the inferior olive (IO) was also identified as important by the PLS and network analysis. The IO is a key source of instructive input to the Purkinje cells of the cerebellum that process information and serve as its primary output (De Zeeuw et al., 1998; Dohaku et al., 2019; Pose-Méndez et al., 2023). A role for the cerebellum in fear-induced freezing is consistent with prior work in goldfish (Yoshida et al., 2004), and recent views of the cerebellum playing an important role in socio-emotional behavior alongside its more well established involvement in motor coordination (Strata, 2015). In looking closer at the functional connectivity of the cerebellum in freezers, we find connections between the CC, IO, and Dm where the Dm is thought to be the teleostan homolog of the basolateral amygdala (Mueller, 2022; Porter and Mueller, 2020) which has been implicated in fear memory learning in both zebrafish and rodents (Fanselow and LeDoux, 1999; Lal et al., 2018; Sun et al., 2020).

Another collection of brain regions that stands out as key to freezing behavior are the posterior tuberculum and pretectum. Within the posterior tuberculum, the preglomerular nuclei are amongst the top contributors to freezing in the PLS and network analysis (particularly the PGm which is supported by both sets of analyses). Preglomerular nuclei receive a variety of sensory inputs and project to the pallium, making it functionally analogous to the mammalian dorsal thalamus (Mueller, 2012; Wullimann, 2020). In particular, the PGm in goldfish, a related teleost species, receives chemosensory input (Northcutt, 2006). Thus, the PGm may work together with other chemosensory regions to recall the CAS memory, like the nLOTa, nLOTp, and Dp that positively co-vary with freezing behavior. Pretectal nuclei, like the PSp and PPv, were also identified as important by both the PLS and network analysis. The pretectum is key to visual processing as these regions receive retinal input and have reciprocal interactions with the optic tectum (Yáñez et al., 2018). Interestingly, several pretectal nuclei, including the PPv, receive or send projections to the CC of the cerebellum (Yáñez et al., 2018). Taken together, these data suggest that the freezing behavior in the absence of elevated evasive behavior, may be driven by a bias towards generating memories from combined visual and chemosensory inputs. These sensory memories then engage a cerebellar-telencephalic circuit to initiate the coordination of freezing behavior. However, further anatomical and functional manipulations would be required to more directly draw these conclusions.

In the functional brain network of evading freezers, we see greater engagement of the hypothalamus and preoptic areas than in freezers. For example, the network analysis identified the SC, PPp, and PPa, parts of the preoptic area, as high in degree and for having a higher proportion of unique edges than other brain regions. These regions formed a cluster that has many edges to hypothalamic and subpallial regions (Figure 8A). The preoptic area of zebrafish has been found to be a neurosecretory area, roughly homologous to the paraventricular nucleus in mammals (Herget et al., 2014). This region, along with other hypothalamic nuclei, modulate a broad spectrum of behaviors and physiological processes, including the hypothalamic-pituitary-interrenal (HPI) axis that regulates the release of cortisol, a central stress hormone (Biran et al., 2015; Wendelaar Bonga, 1997). Relatedly, CAS has been found to increase cortisol in zebrafish (Abreu et al., 2017); and in the frillfin goby (Bathygobius soporator) cortisol modulates the behavioral effect of CAS (Barreto et al., 2014). Thus, it may be the case that individual variation in the status of the HPI axis in zebrafish may determine if zebrafish exhibit a combination of evasive and freezing behaviors (i.e., evading freezers) or freezing interspersed with normal swimming (i.e., freezers).

In conclusion, we have identified the presence of individual behavioral variation in how zebrafish express fear memories. We find that the strain and sex of fish affect fear responses and that the behavioral response during learning is predictive of memory expression. Behavioral variation is associated with differences in the activity of brain regions like the cerebellum, preglomerular nuclei, preoptic area and hypothalamus. Future work will be needed to probe the function of these circuits and determine how they may be differentially wired to give rise to variation in behavior.

Materials and Methods

Subjects

Animals were female and male zebrafish from the AB, TU, TL and WIK strains 6 to 9 months of age. Fish were bred in-house and raised at Wayne State University and were within two generations of fish obtained from the Zebrafish International Resource Center (ZIRC - University of Oregon). Fish were housed in a re-circulating system (Aquarius Fish System, Aquatic Enterprises) with a 14/10-hour light/dark cycle (lights on at 8:00 a.m.; pH 7.5 ± 0.2 °C; water conductivity 500 ± 10 μS; temperature 27.5 ± 0.5 C). Fish were fed two times a day with dry feed in the morning (Gemma 300, Skretting) and brine shrimp in the afternoon (Artemia salina; Brine Shrimp Direct). All the procedures were approved by the Wayne State University Institutional Animal Care and Use Committee (Protocol ID: 21-02-3238).

For the behavioral experiments, animals were housed in male/female pairs (2 female/male pairs per tank separated by a transparent divider) in a 2L tank for one week before behavioral manipulations. The pair-housing system minimizes potential stress induced by social isolation (Parker et al., 2012), and maintains a consistent environment while allowing for the identification of individual fish across days. Behavioral experiments were performed between 11:00 and 15:00 h. Fish were moved to the behavioral room 1 hour before experiments to allow animals to adapt to the new environment. The fear conditioning task consisted of 3 days of behavioral testing: day 1: habituation; day 2: conditioning; day 3: memory recall (Figure 1A).

After behavioral recording, animals remained in the room for 1 hour before being moved back to the re-circulating system. Sex was determined using three secondary sex characteristics: body shape, color, and presence of tubercles on the pectoral fins (McMillan et al., 2015). At the end of the experiment, we confirmed sex after euthanization by dissection to identify the presence of eggs. To ensure data reliability, eight independent batches of offspring for each strain were used for behavioral testing (n = 12 – 16 per strain in each batch). All behavioral testing was carried out using a counterbalanced design across sex and placement in tank (front/back).

Fear conditioning

The fear conditioning protocol was adapted from prior work (Kenney et al., 2017) to use the Zantiks automated system (Zantiks Ltd) and conspecific alarm substance (CAS) as the unconditioned stimulus (Maximino et al., 2018). Two tanks (20 x 7 x 9 cm) were recorded simultaneously during the experiments at 60 fps. A pilot study was performed to determine the best CAS concentration starting with previous work in which a concentration of 3.5mL/L elicited fear-related responses up to 7 days after exposure (Maximino et al., 2018). We found a 1:10 dilution of this dosage (0.35 mL/L) elicited a robust freezing response with some variation on memory day (Figure S2-1). Fish were tracked using DeepLabCut (Mathis et al., 2018) by labelling three points (head, trunk and tail). We used ResNet101 network and trained a deep neural network on 200 manually labelled frames (divided across the four strains and both sexes), the network was then refined and improved by correcting outliers to include 578 additional frames.

CAS extraction

CAS is a substance released by fish when their epithelial cells are damaged, typically in response to a strike from a predator. We extracted CAS as previously described (Egan et al., 2009; Quadros et al., 2016). Briefly, CAS was obtained from donor fish that were euthanized via rapid cooling in 4 °C water. Fish were placed in a 10 cm petri dish and epidermal cells were damaged via 10 shallow incisions on both sides of the body using a surgical blade on ice without drawing blood. 10 mL of distilled water was then added to the dish and it was gently agitated to ensure complete coverage of the lacerated areas on the fish body. Because 4 different strains were used for the behavioral experiments, we prepared stock solutions of alarm substance by mixing extract from 1 fish per sex/strain and then aliquoting before storing at -20 °C for no longer than 4 weeks.

Behavioral analysis

To automate capture of behavior we developed a random forest machine learning model. After fish were tracked using DeepLabCut, we calculated postural information for each individual, extracting 34 parameters (Table S1-1) over a sliding window of 750 ms (45 frames). We manually labelled 14,199 frames per behavior for training the following behaviors: (1) freezing, the absence of movement with fast opercular movements, (2) erratic movements, zig-zagging movement of at least 3 quick C-bends, (3) burst swimming, a rapid increase of speed in one direction, 4) straight swimming, normal speed swimming in one direction and, (5) normal turns, (Figure 1B). We used an 80/20 (train/test) split for training and testing the model.

Optimization of the model was done by testing a range of variables per level (1-15) and trees. We found that 10 variables and 2,000 trees yielded the lowest out of bag error rate (2.88%). Results of testing on 20% of unseen data also yielded excellent results with errors of no more than 6% and typically much lower (Figure 1C).

To analyze the temporal organization of zebrafish behavior, we modeled behavioral state transitions using discrete-time Markov chains. Behavioral sequences were determined by calculating the modal behavior in 750 ms time windows. This time window matched the length of the sliding window used in the random forest classifier. Transition matrices were calculated using maximum likelihood estimation using the markovchain package in R (Spedicato, 2017).

Matrices were fit using data pooled across all animals within each experimental group and day. Ninety-five percent confidence intervals were computed and used for determining significance.

Identification of behavioral clusters

Clusters were identified by creating a k-nearest neighbor network and then applying the Louvain community detection algorithm (Blondel et al., 2008). We only used behavior on memory day from fish exposed to CAS. To calculate distances for the k-nearest neighbor network, we first standardized the behavioral data (% freezing, % evasion, and % evasion as a proportion of active time). We then calculated similarity scores between individual fish:

where D is the Euclidean distance between each fish in three-dimensional behavioral space. We selected the optimal k value by calculating internal clustering metrics using the clusterCrit package (version 1.3.0) at each value of k: Calinski–Harabsz index (Caliński and Harabasz, 1974), silhouette index (Rousseeuw, 1987), and Davies–Bouldin index (Davies and Bouldin, 1979). Based on these data, we chose k = 74, which fell within a range that optimized internal clustering metrics and showed robustness to small changes in k (Figure S4-1). Clustered data were visualized using a uniform manifold approximation (UMAP) (McInnes et al., 2020) as implemented in the UMAP R package (version 0.2.10).

To identify clusters in new data, such as during pre-exposure and exposure or in fish used for brain imaging analysis, we first standardized the new data using the parameters from the initial clustering described above. We then assigned clusters to each fish based on the number of connections to its 33 nearest neighbors in the initial network. Thirty-three was chosen because it is half the size of the smallest cluster identified.

Tissue clearing and staining

Brains were collected 15 min after the end of fear memory expression. Fish were rapidly euthanized via submersion in ice cold water and heads were removed and fixed in 4% paraformaldehyde (PFA) overnight. Animals that were assigned the wrong sex (determined post-mortem as described above) were removed from data analysis, as well as its paired fish (< 2%). After overnight fixation, brains were carefully dissected in cold PBS. Brains damaged during the dissections were discarded. Tissue staining and clearing was performed using iDISCO+ (Renier et al., 2016) as previously described (Rajput et al., 2025). After dissections, samples were washed three times with PBS (30 min washes) at room temperature followed by dehydration using methanol/water (20, 40, 60, 80 and 100% methanol for 30 min each).

Samples were placed in 100% methanol for one hour and chilled in the fridge (4° C). They were then incubated at 4° C overnight in 5% hydrogen peroxide in methanol. The next day, samples were rehydrated in methanol/water (80, 60, 40 and 20%, 30 min each) followed by two PBS washes, one PTx.2 (PBS with 0.2% TritonX-100) of one hour each and an overnight incubation in 5xSSCT buffer (sodium chloride sodium citrate and 0.2% TritonX-100). The following day, an acetylation step was performed with 0.25% v/v acetic anhydride for 30 min followed by three washes with ultrapure water for 5 min each. Samples were then incubated in hybridization buffer (Molecular Instruments) for 15 min at room temperature followed by one hour incubation at 37° C with the same buffer. A solution of hybridization buffer and cfos probes (2 pmol per sample) was then prepared and samples were incubated with the probes for 3 days. Samples were then washed three times using probe wash buffer (Molecular Instruments) followed by two 5xSSCT washes of one hour each. Finally, samples were incubated for 1 h with amplification buffer (5×SSC, 0.1% Tween 20, 10% dextran sulfate) and then incubated in the dark for 3 days in 125 µL of amplification buffer and 2.5 µL of each B1 hairpin (7.5 pmol). After amplification, samples were kept in the dark and washed five times in 5xSSCT buffer and left overnight. Clearing then began with dehydration in methanol/water mixtures (20, 40, 60, 80, 100, and 100%) for one hour each, followed by three hours incubation in dichloromethane/methanol (66%/33%) solution. Samples were then incubated two times for 15 min with dichloromethane. Lastly, the solution from the tubes was removed and dibenzyl ether was added to the samples in which they remained until imaging.

Whole-brain imaging

Cleared samples were imaged via light-sheet microscopy on an UltraMicroscope II (Miltenyi Biotec). Samples were immersed in dibenzyl ether and affixed to a pedestal using a refractive index matched ultraviolet cured resin (NOA 61, Norland Products). Images were acquired at 5× magnification. Illumination of samples was from the right laser using dynamic horizontal focusing and a step size of 4 µm. Autofluorescence images were acquired using a 480 nm laser and used for registration. Cfos images were acquired using an illumination laser of 641 nm. The exposure time was 200 ms, sheet width of 50%, and a sheet numerical aperture of 0.156. Images were acquired as 1 × 2 tiles and stitched using TeraStitcher (Bria and Iannello, 2012).

Automated cell detection

Detection of cfos positive cells was done using CellFinder, version 0.4.20 (Tyson et al., 2021) as previously described (Rajput et al., 2025). Briefly, cell candidate detection was first performed and those cell candidates that were closer than 9 µm to one another were removed. We then applied the CellFinder artificial neural net supervised machine learning algorithm to differentiate cells from non-cells. To train the CellFinder network, we labelled 25,460 cells and non-cells and trained the network to an accuracy of 97.7%. For one brain, cell counting was unsuccessful in the telencephalon and so these regions were excluded from the analysis.

Registration to AZBA

Images were registered to the adult zebrafish brain atlas (AZBA) using advanced normalization tools (version 2.4.4; Avants et al., 2011) to automatically parcellate brains into individual brain regions as previously described (Rajput et al., 2025). We initially generated a template from autofluorescence images that was a combination of a previously generated template (Rajput et al., 2025) and 10 brains from the present study. The autofluorescence image from AZBA was then registered to the template and the transformation from each individual brain to the template was used to bring the AZBA autofluorescence image into the space of individual brains. We had poor registration of the olfactory bulbs in three brains and so these regions were excluded from the analysis for these samples.

Partial least squares (PLS) analysis

We implemented PLS analysis (Krishnan et al., 2011) in R to identify covariation between behavior and regional brain activity. Because of batch-to-batch variation in the brightness of cfos images, cell counts were normalized based on the median overall cfos count for each cohort. This normalized data was then used for PLS (n = 87). Only fish exposed to CAS during training day were included in the analysis. Using mean-centered behavioral PLS correlation, we applied a singular value decomposition to a matrix of correlations between normalized cfos counts and four behaviors: % normal swimming, % freezing, % evasive swimming, and % evasive swimming while active (i.e., not freezing). Because the behavioral data is measured as percentages, some of which were skewed, we used an arcsine square root transformation prior to the PLS analysis:

To determine which latent variables to consider, we applied permutation tests (resampled 10,000 times without replacement) and used this to calculate p-values associated with each singular value. To determine the reliability of brain region contributions to a specific latent variable, we used bootstrapping (resampled 10,000 times with replacement). Bootstrap ratios were calculated by dividing the salience associated with each brain region by its standard error calculated from bootstrapping. Bootstrap ratios are similar to z-scores and thus we consider a threshold of ± 2.5 as significant, which is roughly equivalent to a two-tailed p-value of 0.01.

Brain network analysis

Brain networks were derived from Pearson correlation matrices capturing covariation between cfos counts per brain region across animals in each experimental group. Matrices were thresholded using a false discover rate (FDR) ≤ 0.001 (Benjamini and Hochberg, 1995). For calculating the initial p-values from Pearson’s correlations, we used the lowest number of observations (n = 16) from the groups used for making networks. We made this conservative correction because the p-value is highly sensitive to the number of observations, which would cause a large skew in the number of significant edges identified for each group. To capture edges that were unique to each network, we identified edges whose 95% confidence intervals did not overlap between networks. Only those edges that were significant (i.e., FDR ≤ 0.001) and did not have overlapping confidence intervals (CI) were labelled as ‘unique’. Those with only one overlap were labelled as ‘1 overlapping CI’.

Clustering of graphs was done by identifying the partition that maximized modularity. Relative degree scores for evading freezer and freezer networks were calculated using the following equation:

Where D_i,Fr and D_i,EF are the degree of the ith node for the freezer and evading freezer groups, respectively. Networks were visualized using Cytoscape version 3.10.2 (Shannon et al., 2003; Su et al., 2014).

Coding and statistics

Analysis was performed using R version 4.3.3. To identify differences at specific time points during behavior we FDR (Benjamini and Hochberg, 1995) corrected independent samples t-tests and compared CAS and vehicle treated groups. One sample t-tests were used to compare difference scores to zero with p-values also corrected using the FDR. The markovchain package (version 0.9.5) was used to create the Markov chains for generating transition matrices to identify behavioral transitions (Spedicato, 2017). Missing values for brain activity were imputed via nuclear-norm regularized singular value decomposition as implemented in the softImpute package (version 1.4.3; Hastie et al., 2015; Mazumder et al., 2010). The regularization penalty (λ = 25) was determined empirically by identifying the inflection point in the scree plot of the singular values. Permutation resampling without replacement (10,000 times) was used to assess significance of strain and sex for each behavioral cluster. We used the igraph (Csardi and Nepusz, 2006; version 2.02) package for network analysis and the cccd (version 1.6) package to build the k-nearest neighbor networks.

Data availability

All data needed to evaluate the conclusion in the paper are present in the paper and the supplemental material. Supplemental files are available at Dryad: https://doi.org/10.5061/dryad.sqv9s4nfh:

File S01. Statistics for Figure 2A.

File S02. Statistics for Figure 2B.

File S03. Behavioral data for Figure 2.

File S04. Neural cfos and behavioral data for Figure 6.

File S05. Bootstrap ratios for Figure 6.

File S06. Neural cfos data for network generation for Figure 7.

File S07. Correlations and confidence intervals for freezer, evading freezer, and non-reactive networks.

File S08. Graphml file for evading freezer network in Figure 8A.

File S09. Graphml file for freezer network in Figure 8A.

File S10. Graphml file for non-reactive network in Figure S8-1.

File S11. Graphml file for vehicle network in Figure S8-2.

File S12. Cytoscape file containing all four networks (evading freezer, freezer, non-reactive, and vehicle).

Acknowledgements

We thank Ian Lee for copy editing of the manuscript.

Additional information

Funding

Funding was provided by the National Institute of General Medical Sciences (R35GM142566, JWK).

Author contributions

Conceptualization: BDF and JWK; Methodology: BDF, JH, and JWK; Investigation: BDF, JH, NR, DK and DL; Formal analysis: BDF, NR and JWK; Writing: BDF and JWK.

Funding

HHS | NIH | National Institute of General Medical Sciences (NIGMS) (R35GM142566)

• Justin W Kenney

Additional files

Supplemental figures. Figure S2-1. Freezing behavior in response to different concentrations of CAS. Figure S4-1. Internal clustering metrics for K-nearest neighbor partitions. Figure S6-1. Bootstrap ratios for PLS analysis. Figure S6-2. Example covariations between cfos activity in individual brain regions and behavioral contrasts. Figure S8-1. Functional network for non-reactive zebrafish. Figure S8-2. Functional network for vehicle zebrafish.

Supplemental tables. Table S1-1. Parameters used for training random forest machine learning model. Table S6-1. Brain region abbreviations.