Introduction

Animals execute sophisticated behavioral escape responses to avoid environmental threats. These movements are coordinated by specific sensorimotor circuits, which receive inputs via specialized sensory neurons that are tuned to respond to environmental cues of ethological relevance. Visual pathways, for example, are tuned to detect the motion or change in size of nearby objects [1], and downstream motor circuits coordinate action sequences that serve to maintain a balance between responding with sufficient strength to achieve the desired goal and avoiding an overreaction that would exhaust the animal and its metabolic reserves [2]. Critically, these sensorimotor processes must be in place and ready to function the moment the animal engages with the world, which in oviparous species can occur very early in development.

Animal survival also depends on sophisticated internal organ responses to environmental threats. For example, cardiac output varies in response to environmental challenges, allowing adaptive oxygen delivery during vigorous behaviors. This regulation occurs at sub-second timescales and requires the activation of neural circuits, which enable fast and bidirectional communication between the brain and internal organs [3,4]. However, the functional development of these circuits remains poorly understood, largely due to technical limitations. In rodents, for instance, monitoring or manipulating heart-brain circuits often requires invasive procedures performed under anesthesia, which disrupts the very homeostatic responses being studied [5,6]. Moreover, in most mammals, including rodents and humans, autonomic innervation of the heart and other viscera begins in utero, precluding longitudinal functional imaging during the initial establishment of heart-brain connectivity [7]. An exception is found in marsupials, where autonomic innervation occurs postnatally [8], offering rare developmental access, but these species are experimentally challenging. Thus, a model that combines ex utero development, optical accessibility, and intact physiology is ideal to understand how autonomic circuits emerge and function during early life.

Larval zebrafish offer, for several reasons, an ideal alternative to studying the functional development of heartbrain neural circuits in intact animals. First, zebrafish larvae are small, relatively transparent, and amenable to developmental, behavioral, and physiological studies using only modest restraint. Second, the structure of the neural circuits that connect internal organs and the brain is conserved in bony fish, reptiles, birds, and mammals [9]. In all these vertebrates, organ regulation is carried out by the motor vagus nerve and the sympathetic system, which send critical control signals from the brain and the spinal cord to the organs. Similarly, in those vertebrates, most viscerosensory signals are carried by the sensory vagus nerve to the brain [10]. Third, cardiac responses can be measured and manipulated in larval zebrafish with optical methods. These advantages are crucial but not sufficient; here, we establish hardware and software tools, as well as data analysis pipelines and transgenic animals that allow us to study the functional development of heart-brain motor and sensory circuits.

Previous studies of the adult zebrafish heart have shown that stimulation of the motor vagus nerve induces bradycardia, and sympathetic nerve stimulation induces tachycardia [11, 12]. Moreover, in vitro preparations of the zebrafish heart have been developed and used to study neural effects caused by pharmacological agents such as vapor anesthetics [13]. Experiments with adult zebrafish, however, have some of the same technical limitations as experiments with rodents, including a lack of optical access, the need for invasive procedures, and the inability to concurrently observe autonomic neural dynamics. Here, we establish larval zebrafish as a model for heart-brain interactions by characterizing the function of motor and sensory circuits of the heart during development, finding the earliest stages in which it is possible to study the sensorimotor control of the heart while retaining optical access to both the heart and brain.

We establish larval zebrafish as an in vivo model for longitudinal, cell-resolved physiology of heart-brain circuits. With custom instrumentation that enables calcium imaging, optogenetics, and quantitative all-optical cardiac physiology, we show that cardiac control assembles in discrete steps. Early in development, a minimal efferent circuit is sufficient for heart rate control. Direct brain-to-heart projections execute early heart rate control even in the absence of intracardiac neurons. At this stage, heart rate is governed exclusively by the motor vagus nerve. Within the hindbrain, we functionally localize a vagal premotor nucleus that drives this early control.

As development proceeds, sympathetic innervation arrives and enhances the amplitude and kinetics of cardiac responses. Neurons in the most anterior sympathetic ganglia acquire the capacity to accelerate the heart and exhibit proportional, integral, and derivative-like relationships to heart rate, indicating controller motifs that shape gain and dynamics. Soon after, vagal sensory neurons (VSNs) innervate the heart. Distinct VSN subsets increase activity when heart rate falls or rises. These VSNs encode heart rate during spontaneous fluctuations, responses to aversive stimuli, and optogenetically evoked cardiac perturbations. In all of those cases, their dynamics are captured by one temporal basis function whose phase-shifted instances tile heart rate trajectories, consistent with a population code for heart rate.

Our findings uncover how the vertebrate nervous system progressively assembles a closed-loop control architecture for cardiac regulation, beginning with descending motor control and culminating in vagal sensory neuron innervation of the heart. Human autonomic disorders, such as paroxysmal orthostatic tachycardia syndrome (POTS) and long COVID-associated dysautonomia, are hypothesized to involve disruptions in feedback control of the heart, including altered sensory gain and defective reflex adaptation [14, 15]. The framework we present enables in vivo testing of these hypotheses at cellular resolution by allowing precise manipulation of individual system components and direct observation of resulting deficits in cardiac dynamics.

Results

The Emergence of a Minimal Functional Circuit for Neurocardiac Control

The two-chambered heart of the zebrafish embryo starts beating 24 hours post-fertilization, with atrial contraction followed by ventricular contraction, driven by its intrinsic pacemaker [16, 17]. To investigate the emergence of a minimal circuit for neurocardiac control in larval zebrafish, we developed a system, as shown in Figure 1A, that enables the quantification of the heart’s physiological function while simultaneously exposing the animal to visual environmental stimuli. In our setup, the larval zebrafish head is immobilized using agarose, while the tail is free to move (Supp. Fig. 1G). We quantified individual atrial or ventricular contractions based on the reflection of near-infrared light, as shown in Figure 1B. Using this approach, we examined the spontaneous heart rate of fish during development and found, consistent with prior work [18, 19], that basal heart rate first decreases from 4 to 5 days post fertilization (dpf) and increases from 5 to 7 dpf (Supp. Fig. 1B).

Figure 1.

To induce transient changes in heart rate, we used a brief UV flash (Fig. 1C), a flash of darkness (Fig. 1D), or looming stimuli (Supp. Fig. 1C-F), which are well-studied visual threats. We subjected larval zebrafish to these stimuli across developmental stages ranging from 4 to 12 dpf. We found that while at 4 dpf, fish are capable of behavioral escapes as marked by tail flicks (Supp. Fig. 1G), their heart rate remained constant (Fig. 1C, D, Supp. Fig. 1C, H). At 5dpf, larval fish start displaying increased heart rate responses to UV flash, dark flash, and looming (Fig. 1C, D, Supp. Fig. 1D, H). These responses change their temporal dynamics and intensity as the fish matures. By 12 dpf, the response is larger, more rapid in onset and adaptation, and more stereotypical than at earlier stages (Fig. 1C, D, Supp. Fig. 1F, H).

To define the minimal neural structure required for heart rate modulation, we examined the spatiotemporal development of cardiac innervation using immunostaining for acetylated alpha-tubulin to label neurites and Myl7 to label the myocardium (Fig. 1E). We observed that cardiac innervation begins at 5 days post-fertilization (dpf), when the first neurites reach the entry of the atrium, known as the sinoatrial plexus (SAP), a region that includes the pacemaker cells [20,21]. As development proceeds, innervation expands in a stereotyped sequence: the atrium is innervated by 7 dpf, followed by the atrioventricular plexus (AVP) near the valve region, and the ventricle by 12 dpf.

Next, we asked whether intracardiac neurons are required for the initiation of heart rate control. To track their development, we used transgenic zebrafish expressing GFP under the control of the phox2bb promoter, which labels all intracardiac neurons (ICNs) [22]. Consistent with our acetylated α-tubulin staining, we observed no cardiac innervation at 4 dpf (Fig. 1F). By 5 dpf, neurites begin to reach the heart, but no intracardiac neurons are yet present (Fig. 1F; Supp. Vid. 1). The first intracardiac neuron migrates into the heart by 7 dpf (Fig. 1F; Supp. Vid. 2), and by 12 dpf, multiple intracardiac neurons are evident (Fig. 1F; Supp. Vid. 3). To visualize neuronal somata more clearly, we performed immunostaining for Phox2bb alongside Myl7 to label the myocardium (Fig. 1G). At 5 dpf, no neuronal cell bodies were observed within the heart. By 12 dpf, we identified 6 neurons near the SAP, 6 near the AVP, and a cluster of 16 neurons in the dorsal-posterior atrium, which we designate the Intracardiac Ganglion (ICG).

To summarize the timing and variability of cardiac innervation and neuronal migration, we quantified developmental progression across animals (Fig. 1H, I): at 4 dpf, 0/12 fish showed cardiac innervation; at 5 dpf, 12/15 fish showed innervation; at 7 dpf, 14/20 fish had at least one intracardiac neuron; and by 12 dpf, all 12/12 fish displayed robust innervation of both cardiac chambers with multiple intracardiac neurons. Together, these results show that functional cardiac innervation precedes the arrival of intracardiac neurons, establishing that early modulation of heart rate occurs through direct descending motor pathways alone.

Formation of the Direct Brain-to-Heart Motor Circuit

We next sought to determine which branches of the autonomic nervous system contribute to progressive innervation of the heart between 5 and 12 dpf. Because the motor vagus nerve is cholinergic [10], we began by labeling it with antibodies against choline acetyltransferase (ChAT), the enzyme responsible for synthesizing acetylcholine, and labeling the myocardium with antibodies against Myl7. This approach allowed us to visualize cholinergic innervation of the heart beginning at 5 dpf (Fig. 2A).

Figure 2.

We also detected ChAT-positive intracardiac neurons (ICNs), including the first neuron to arrive at 7 dpf, suggesting that this cell is a postganglionic neuron of the motor vagus nerve (Fig. 2A–C). Motor vagus postganglionic neurons are known to be cholinergic and receive input from vagal preganglionic fibers and form direct connections with cardiac tissue [23]. By 12 dpf, we observed six ChAT-positive ICNs near the sinoatrial plexus (SAP). We confirmed the cholinergic identity of the first ICN using fluorescent in situ hybridization chain reaction (HCR) for ChAT, together with GFP boosters to confirm overlap with ChaTA-Gal4, UAS-GFP and phox2bb:GFP (Supp. Fig. 2). We did not detect neurons by the atrioventricular plexus (AVP) or elsewhere, suggesting that the other ICNs identified in Fig. 1G may use alternative neurotransmitters.

We then used ChaTA-Gal4, UAS-GFP transgenic fish to label motor vagal fibers selectively. Because red fluorophores in the heart interfered with the detection of thin axons, we instead imaged the myocardium in transmitted light using a photodiode-based detector, while capturing green fluorescence in a separate channel. This strategy enabled reliable identification of even fine cholinergic neurites, based on their anatomical position and rhythmic motion during the heartbeat. Consistent with our immunostaining and phox2bb:GFP imaging, we observed no innervation at 4dpf (Fig. 2D) (Supp. Vid. 4), emerging innervation at 5 dpf (Fig. 2D, Supp. Vid. 5), emergence of the first ICN at 7 dpf (Fig. 2D, Supp. Vid. 6), and further extension to the AVP by 12 dpf (Fig. 2D, Supp. Vid. 7). We did not detect cholinergic innervation of the ventricle at any stage up to 12 dpf.

To quantify the developmental progression of cholinergic cardiac innervation across individual animals, we scored each fish for the presence of cardiac innervation and the number of cholinergic intracardiac neurons (ICNs) observed at 4, 5, 7, and 12 dpf (Fig. 2E). At 4 dpf, the vast majority of fish (14/15) lacked any evidence of cholinergic cardiac innervation. By 5 dpf, most fish (12/14) displayed innervation of the sinoatrial region, but none had ICNs. At 7 dpf, 11 of 17 fish showed a single cholinergic ICN, while the remaining 6 had innervation but no ICNs. By 12 dpf, all fish exhibited robust innervation of the atrium and SAP, with 11 of 14 fish containing between two and six ICNs. These data highlight a gradual and stage-specific acquisition of cholinergic innervation and ICNs, with clear stereotypy in the timing of each step in the maturation of the vagal motor circuit.

To investigate if the first ICN already receives neural inputs, we performed serial-section electron microscopy (ssEM) of the SAP in a 7 dpf fish (Fig. 2G). We identified five cell types—cardiomyocytes, epithelial cells, blood cells, endothelial cells, and one neuron, consistent with our light microscopy observations (Fig. 2G). Vesiclerich synaptic contacts were found both onto the ICN (Fig. 2H) and from the ICN onto cardiomyocytes (Fig. 2I, J), supporting its integration as a functional node in the developing cardiac motor circuit.

A Vagal Premotor Nucleus in the Hindbrain Controls Heart Rate

To explore how this first anatomical evidence of innervation correlates with the onset of physiological control, we used calcium imaging to measure the activity of motor vagus neurons in the brain. Previous anatomical studies have shown that motor vagus neurons in larval zebrafish are cholinergic and located in the hindbrain [24]; thus, we focused on that population. In order to characterize the dynamics of those neurons, we customized a 2-photon imaging rig with an optical path for simultaneous quantification of cardiac function using near-infrared light and brain functional imaging at single-cell resolution using the 2-photon laser (Fig. 3A). Using ChaTA-Gal4 to drive UAS-GCaMP6s expression, we selectively imaged cholinergic hindbrain neurons while concurrently monitoring heart rate. Dark flash stimulation elicited robust tachycardic responses in 5, 7, and 12 dpf fish (Fig. 3C-E) but not in 4dpf fish (Fig.3B), consistent with our results in Fig. 1. Neural responses were diverse and even when large fractions of neurons displayed synchronous increases in activity, those were not consistently accompanied by changes in heart rate (Fig. 3B-E).

Figure 3.

We asked whether motor-vagus activity carries a time-resolved signature of cardiac control. To do so, we fit a generalized linear model (GLM) that decodes heart-rate (HR) dynamics from the population activity, assigning each neuron a causal finite-impulse-response kernel that quantifies its influence across time lags. A static sigmoid nonlinearity captured saturation in HR output. To visualize the temporal motifs, we pooled predictive neurons—those exceeding an Otsu threshold on cross-validated predictive power (Methods, Decoding GLMs)—and plotted their kernels as a shape-sorted heatmap (Fig. 3F). This revealed a single canonical waveform that recurred across cells with opposite polarities: some neurons expressed the positive version and others the negative of the same kernel (Fig. 3F-H).

Model performance was evaluated with blocked 5-fold cross-validation (80% train, 20% test per fold; Methods). Across 31 fish (5, 7, and 12 dpf), only 2 fish at 5 dpf and 2 at 7 dpf fell below CV R² = 0.20; in all others, the model predicted HR trends with CV R² ranging from 0.25 to 0.60 (Fig. 3I). Even for animals near the cohort median (CV R² = 0.35 and 0.43), the predicted traces closely tracked stimulus-evoked HR responses (Fig. 3J). Although predictive neurons (defined by the Otsu threshold on CV performance) comprised <5% of recorded cells at each age, their prevalence tripled at 12 dpf relative to 5 dpf (Fig. 3K). Among encoders, the fraction with negative kernels—the HR-decreasing sign—rose from ~40% at 5 dpf to ~90% at 12 dpf, consistent with the cardioinhibitory role of the motor vagus (Fig. 3K). Spatially, encoders were concentrated medially across ages; positive-kernel encoders were more anterior and dispersed, whereas negative-kernel encoders formed a reproducible medial cluster (Fig. 3L–O).

Given the stereotyped, medial clustering of negative-kernel encoders and the established cardio-inhibitory role of the motor vagus, we hypothesized that activating this population would reduce heart rate. We expressed the light-gated channel CoChR under ChaTA-Gal4 in cholinergic hindbrain neurons and delivered patterned photostimulation to the negative-kernel locus using structured illumination steered by dual galvos (Fig. 3P,Q). At 4 dpf, before anatomical heart innervation is detectable, photoactivation of this locus did not measurably alter heart rate (Fig. 3R). By 7 dpf, once innervation is present, the same stimulation decreased heart rate (Fig. 3S). As a control, optogenetic activation of the anterior spinal cord did not result in changes in heart rate (Supp. Fig. 3). These data indicate that a small, medial subset of cholinergic hindbrain neurons functions as a vagal premotor nucleus for cardiac control.

Formation of Sympathetic Cardiac Circuits

In zebrafish, the sympathetic nervous system consists of a chain of paravertebral ganglia that extends from the visceral cavity to the tail (Fig. 4A). As in other vertebrates, these ganglia can be visualized using antibodies against tyrosine hydroxylase (TH), the rate-limiting enzyme in catecholamine biosynthesis [26]. To determine when the heart first receives sympathetic input, we performed TH immunostaining alongside Myl7 to label myocardial tissue. At 5 dpf, no TH-positive fibers were detected in the heart (Fig. 4C). By 7 dpf, TH-positive projections reached the sinoatrial plexus (SAP), indicating the onset of sympathetic cardiac innervation (Fig. 4C). These fibers originated from the most anterior sympathetic ganglia, which we refer to as the anterior paravertebral ganglia (APG) (Fig. 4C). Over development, the APG expand and adopt a distinctive C-shaped morphology (Fig. 4B). Despite this morphological maturation, TH-positive projections remain restricted to the SAP through 12 dpf, with no evidence of innervation to the atrioventricular plexus or ventricle.

Figure 4.

To quantify the developmental dynamics and interindividual variability of sympathetic innervation, we assessed individual animals at 5, 7, and 12 dpf and recorded the presence or absence of TH-positive fibers at the heart. At 5 dpf, none of the nine animals examined showed sympathetic innervation. By 7 dpf, 6 of 13 animals exhibited TH-positive fibers at the SAP, indicating a variable onset of innervation across individuals. By 12 dpf, all 10 animals displayed consistent TH-positive projections to the SAP (Fig. 4D, E). These findings reveal a gradual and variable onset of sympathetic connectivity to the heart, with full anatomical incorporation of the APG circuit established by 12 dpf.

Emergence of Control Motifs and Function in Sympathetic Cardiac Circuits

To investigate how the anterior paravertebral ganglia (APG) acquire the ability to regulate cardiac function, we expressed GCaMP6s under the TH1-Gal4 driver and we simultaneously imaged sympathetic neuron activity and heart rate dynamics in response to brief UV flashes (Fig. 5A). At 5 dpf, the APG neurons respond to the UV stimulus with transient increases in activity (Fig. 5B), suggesting early engagement of the sympathetic system in fight-or-flight visceral responses even before its anatomical integration with the heart. At 7 dpf, sympathetic neurons respond with apparent faster dynamics (Fig. 5C). By 12 dpf, as the ganglion increased in size, we observed emergent heterogeneity: newly appearing neurons decreased their activity following the same stimulus (Fig. 5D), hinting at the maturation of diverse sympathetic subtypes.

Figure 5.

To further investigate the diversity of neural dynamics and their temporal relationship with heart rate, we trained a generalized linear model (GLM) to predict heart rate from the activity of all recorded APG neurons. This model estimates a temporal kernel for each neuron that captures how its activity at different time lags influences heart rate. Pooling across animals revealed a striking diversity of kernel shapes (Fig. 5E). Unsupervised clustering of kernels identified six reproducible motifs (Fig. 5F, Supp. Fig. 4). Two motifs were biphasic with opposite polarity (C1, C4), resembling derivative filters that emphasize the effect of changes in neural activity. Two were monophasic bumps (C2, C5), consistent with proportional decoding of magnitude. The remaining two were monotonic ramps (C3, C6), consistent with integrator functions that accumulate activity over time. Together, these motifs map onto canonical motor control operations—derivative, proportional, and integrator—that may form the building blocks of cardiac control.

In our decoding model, neural activity is passed through its assigned kernel, summed across neurons, and then transformed by a static sigmoid nonlinearity to predict heart rate (Fig. 5G). Cross-validated performance improved significantly with age, with CV-R² increasing from 0.2 at 5 dpf to ~0.5–0.6 at 12 dpf (Fig. 5H), indicating that APG neurons acquire predictive power as circuits mature. To illustrate the predictive power of the model even with modest CV-R², we display the heart rate over time (in blue) and the heart rate predicted by the model (in orange) (Fig. 5I). Cluster composition also shifted across development: derivative motifs were absent early (except for one neuron) and emerged at 7 and 12 dpf (Fig. 5J), suggesting that the APG progressively establishes the computational repertoire needed for cardiac control.

To investigate whether specific neuron types are spatially organized in the ganglion, we plotted the imaged APGs coloring each neuron by cluster. We found that locations are variable per fish and per side, without a stereotypic location (Fig. 5K-M).

Last, to directly test the causal role of sympathetic neurons in modulating heart rate, we optogenetically activated APG neurons expressing CoChR under the TH1-Gal4 promoter. Using galvo mirrors, we delivered spatially targeted stimulation to the APG while simultaneously recording heart rate (Fig. 5O). At 5 dpf, stimulation failed to evoke any change in cardiac output (Fig. 5P), consistent with the absence of sympathetic innervation. At 7 dpf, however, optogenetic activation induced modest heart rate increases (Fig. 5Q), which became progressively larger by 9 and 12 dpf (Fig. 4R). Importantly, this response was abolished by propranolol, a β-adrenergic receptor antagonist (Supp. Fig. 5), confirming that sympathetic neurons modulate cardiac function via β-adrenergic signaling.

Formation of the Heart-to-Brain Interoceptive Circuit

The somata of the sensory vagus are organized into four ganglia per side. The largest and most posterior is the nodose ganglion. It is located below the lateral line ganglion and innervates the heart and most of the viscera. The three most anterior ganglia are the epibranchial ganglia, which innervate each of the gill arches and are believed to track oxygen and other ventilatory parameters [28] (Fig. 6A). To visualize vagal ganglia, we expressed GFP under the control of the vglut promoter, which labels the cranial sensory neurons [29]. To assess the developmental timing and variability of cardiac innervation by the sensory vagus nerve, we examined the presence of sensory axonal projections to the heart at 5, 7, and 12 dpf. At 5 dpf, all fish (8/8) had formed vagal sensory ganglia but lacked projections to the heart, indicating that ganglion specification precedes target innervation. By 7 dpf, only 2 of 10 animals exhibited cardiac sensory innervation, revealing substantial interindividual variability in the onset of sensory connectivity. By 12 dpf, all fish (11/11) showed robust vagal sensory projections to the heart, demonstrating that sensory innervation becomes fully established at this stage.

Figure 6.

Cardiac encoding in vagal sensory neurons

In order to investigate how vagal sensory neurons may encode heart rate dynamics, we expressed GCaMP6s under the control of the vglut2-Gal4 driver. We selected this driver because in zebrafish, all cranial sensory neurons, including vagal sensory neurons (VSNs), are glutamatergic [29]. We did not restrict our experiments to 12dpf – when vagal sensory innervation of the heart is established – because heart rate information could be available to VSNs from innervation of vasculature, gills, or other indirectly coupled structures. We stimulated fish with UV-flashes while recording calcium activity in the vagal ganglia concomitantly with heart rate (Fig. 7A). We found robust heart rate responses across development, from 5 to 12 dpf, and diverse neural dynamics, including subpopulations of VSNs that displayed temporally locked responses (Fig. 7B-D). Using generalized linear models (GLMs), we asked how well heart-rate fluctuations predict neural activity, assigning each neuron a causal temporal kernel. Neurons exceeding a cross-validated performance threshold (CV R² > 0.2) were classified as cardiac-encoding. The proportion of encoding neurons was comparable at 5, 7, and 12 dpf (Fig. 7E), indicating that subsets of VSNs track heart-rate trajectories even before direct sensory innervation of the heart is fully established.

Figure 7.

Kernel heatmaps aligned by shape revealed that the diversity of response filters could be described as a single temporal motif expressed at systematically shifted phases (Fig. 7F). In other words, rather than multiple distinct kernel families, each encoding neuron appeared to express the same biphasic kernel advanced or delayed in time. Principal-component analysis confirmed this interpretation: the first two PCs accounted for nearly all variance (Supp. Fig. 6) and arranged the kernels along a circular trajectory in phase space (Fig. 7G). This geometry is consistent with a basis set of sinusoidal quadrature components, where neurons are distributed around the circle such that each one samples a different phase of the same underlying waveform (Fig. 7H).

The quality of this representation improved across development. Model performance increased, with crossvalidated (CV) R² values rising from ~0.3 at 5 dpf to >0.5 by 12 dpf (Fig. 7I), indicating that HR became a progressively stronger predictor of neural activity. To illustrate the predictive power of the model even with CVR² values near the median, we display the ΔF/F₀ of a neuron with phase π (top) and a neuron with phase 0 (bottom) over time (in black) and the ΔF/F₀ predicted by the model (in green) (Fig. 7J). In both cases the signal is tracked accurately, with every large fluctuation in measured ΔF/F₀ matched by a large fluctuation in the predicted ΔF/F₀. This is also important to illustrate that neurons with phase π decrease activity when heart rate increases following the dynamics of the canonical kernel, and neurons with phase 0 increase activity when heart rate increases, also following the dynamics established by the canonical kernel.

The distribution of kernel phases also changes systematically during development: at 5 dpf, neurons did not occupy phases uniformly, but instead showed a bias toward 0 and π (Fig. 7K, blue radial histogram). By 7 dpf, these biases strengthened (Fig. 7K, green radial histogram), and by 12 dpf neurons clustered into two modes of representation: a positive biphasic kernel (at phase 0) and a negative biphasic kernel (at phase π) (Fig. 7K red radial histogram). We confirmed this observation using a pairwise permutation test on the kernel phases, which showed that 5dpf fish have a different distribution than fish at 7 and 12 dpf (Fig. 7K inset heatmap), consistent with emergent bimodal encoders with opposite phases.

Anatomical mapping showed no spatial segregation of phase groups: neurons from different ensembles were intermixed within vagal ganglia and across both sides (Fig. 7L–N). Thus, the phase structure of the code is functional rather than anatomical, arising from temporal offsets in kernel phase rather than physical clustering.

Together, these findings demonstrate that vagal sensory neurons implement a distributed, phase-shifted population code for heart rate. By expressing a single canonical kernel at different phase offsets, the ensemble provides a robust and redundant representation of cardiac dynamics, well before direct cardiac sensory innervation is fully stabilized, and, as the fish matures, the encoding shifts to a bimodal representation with kernels of opposite phases.

Stimulus-invariant heart rate encoding in vagal sensory neurons

The observation that a single kernel shape, expressed at shifted phases, accounts for heart rate encoding across fish and developmental stages suggests that this is a core computational property of vagal afferents. To test this hypothesis, we asked whether the same phase-shifted kernel emerges under distinct modes of heartrate modulation: (1) spontaneous fluctuations, which occur without environmental perturbations (Fig. 8A); (2) fluctuations elicited by mosaic optogenetic activation of small subsets of cardiomyocytes (Fig. 8B); and (3) acute tachycardia triggered by brief optogenetic activation of the entire myocardium (Fig. 8C). Spontaneous fluctuations were chosen because, in contrast to fight-or-flight responses, they occur while the external sensory environment is constant and therefore span a broad range of amplitudes, durations, and temporal profiles (Fig. 8D). Optogenetic manipulations were included because they perturb the heart directly and independently of external stimuli. For these experiments, we generated transgenic zebrafish expressing the light-gated cation channel CoChR under the cmlc2 promoter, which labels cardiomyocytes [30]. Mosaic activation produced local perturbations that did not evoke immediate tachycardia but instead generated delayed fluctuations, sometimes tens of seconds after stimulation (Fig. 8E). In contrast, global heart activation produced immediate and robust tachycardia (Fig. 8F). Together, these three conditions provided diverse trajectories of heart rate and neural activity (Fig. 8D–F) without altering the external environment, allowing us to probe encoding mechanisms under controlled and distinct perturbations.

Figure 8.

We next trained GLMs, as in Fig. 7, to identify causal kernels that predict neural activity from heart-rate dynamics. Across all conditions, kernel heatmaps revealed the same smooth green↔magenta gradient (Fig. 8G–I) seen previously under UV stress (Fig. 7F), indicating that most encoding neurons are well described by a single canonical biphasic filter, differing primarily in phase (latency/polarity). PCA of kernel shapes confirmed this structure: in all three conditions, kernels distributed along a circular trajectory in PC space, with neurons clustering near phase 0 and π (Fig. 8J–L), as observed in 12 dpf fish during UV stimulation (Fig. 7K). Encoding strength was similar across conditions, with median CV R² values in the 0.4–0.55 range (Fig. 8M). Moreover, the canonical kernel identified in each case converged to the same shape (Fig. 8N).

As in the neural recordings of responses to UV flashes, anatomical mapping showed no spatial segregation of phase groups: neurons from different ensembles were intermixed within vagal ganglia and across both sides (Fig. 8O–Q). Across fish, we observe a higher likelihood of neurons being located on the ventral side of the nodose ganglion, however, there are also neurons in more dorsal parts of the ganglia (Fig. 7L-N, Fig. 8O-Q).

These results demonstrate that vagal sensory neurons implement a stimulus-invariant, phase-shifted (latency/polarity) population code for heart rate. Whether fluctuations arise from environmental stressors, spontaneously, from localized perturbations, or from optogenetically-induced tachycardia, the dynamics of VSN encoding are well described by a linear time-invariant filter set, with phase (latency/polarity) offset as the primary degree of freedom.

Discussion

Here, we have defined the developmental sequence by which heart–brain functional circuits emerge and mature in a vertebrate. Using larval zebrafish, we combined longitudinal calcium imaging, optogenetics, quantitative optical cardiophysiology, and computational modeling to dissect the limbs of the autonomic circuit. We show that descending vagal motor projections are sufficient for early cardiac control even in the absence of intracardiac neurons; that a small premotor nucleus in the hindbrain provides the first locus for cardioinhibitory drive; that sympathetic inputs from the anterior paravertebral ganglia add modulatory capacity and their dynamics display proportional, integrative, and derivative-like control motifs; and that vagal sensory neurons encode heart rate with a canonical biphasic kernel expressed at distinct phase offsets, yielding a robust, stimulus-invariant population code. Together, these findings reveal how neural regulation of the heart is assembled in stages and how each stage contributes a distinct computational role.

Early brain-to-heart communication via the motor vagus

Our findings demonstrate that the first stage of functional cardiac control in zebrafish is mediated exclusively by descending vagal motor neurons (Fig. 1 and 2). By 5 dpf, cholinergic fibers reach the sinoatrial plexus, and vagal activity becomes capable of altering heart rate even in the absence of intracardiac neurons. Optogenetic activation of hindbrain cholinergic neurons confirmed that these descending projections are sufficient to slow the heart, revealing a minimal efferent circuit. This stage is particularly striking because it separates the earliest onset of neural control from the later incorporation of intracardiac ganglia. In mammals, parasympathetic ganglia embedded in the cardiac plexus are thought to be essential for modulation of pacemaker activity [31, 32], but our results suggest that brain-to-heart projections alone can transiently carry this function. Such a configuration may provide a rapid developmental solution for ensuring that cardiac activity is already under neural control during the critical window when larvae first engage with the external environment.

A premotor nucleus for descending control

In parallel with the onset of vagal motor projections, we identified a small medial cluster of hindbrain neurons whose activity reliably predicts cardiac inhibition. Encoding models revealed a canonical kernel with negative polarity, and targeted optogenetic activation of this locus elicited bradycardia in animals as early as 7 dpf (Fig. 3R). These data indicate that a premotor nucleus in the hindbrain provides a dedicated node for descending cardiac control. The temporal precision and polarity of these kernels suggest that the premotor nucleus exerts a gain-limited inhibitory drive on the heart, consistent with the well-established cardioinhibitory role of the motor vagus. Although the molecular identity of these premotor neurons remains to be fully resolved, their reproducible medial location and robust functional impact indicate that they are a conserved substrate for parasympathetic regulation.

Sympathetic innervation and emergence of control motifs

The second stage of circuit assembly involves the arrival of sympathetic projections from the anterior paravertebral ganglia (Fig. 4). Anatomical innervation emerges and optogenetics show that sympathetic neurons develop the capacity to accelerate the heart at 7dpf. Importantly, encoding analyses of neural activity uncovered kernels with proportional, integrative, and derivative-like temporal motifs. This computational repertoire mirrors classical control architectures, suggesting that the sympathetic system contributes not only to the magnitude but also to the dynamics of cardiac responses. Such motifs may allow the system to respond proportionally to changes in heart rate, integrate sustained deviations, or accentuate rapid fluctuations, thereby shaping the gain and stability of autonomic control. Mammalian baroreflex studies indicate modulation of gain and sensitivity that is consistent with proportional control, and dynamic changes in reflex latency under different physiological conditions, suggesting that components of the control motifs we observe may be conserved [33, 34]. Our findings are the first evidence for heart rate integral- or derivative-like temporal dynamics in sympathetic neurons.

Development of the intracardiac nervous system

Although our primary focus was on the staged emergence of autonomic control, our results also provide new insight into the developmental progression of the intracardiac nervous system (ICNS). By tracking phox2bb+ intracardiac neurons alongside functional assays of heart rate control, we found that cardiac innervation is first detectable at 5 dpf, but intracardiac neurons are absent at this stage. The earliest intracardiac neuron appears by 7 dpf, and by 12 dpf multiple clusters are evident, including a reproducible population in the dorsal–posterior atrium that we designate the intracardiac ganglion. Thus, the arrival of ICNs lags behind the onset of functional vagal efferent control, which is already functional by 5 dpf. This indicates that the initial establishment of parasympathetic influence over heart rate occurs independently of intracardiac ganglia, with these neurons contributing only after the basic circuit is in place.

In mammals, intrinsic cardiac ganglia are thought to be essential for modulation of pacemaker activity and for shaping the spatiotemporal dynamics of cholinergic input to the sinoatrial and atrioventricular nodes [31, 32]. Our findings suggest that in zebrafish, these neurons are dispensable for the earliest phases of functional control. They likely refine or diversify cardiac modulation as they integrate into the circuit. This sequence may reflect a developmental principle: early reliance on extrinsic parasympathetic fibers to ensure inhibitory capacity, followed by later recruitment of intracardiac ganglia for fine-tuned control.

Sensory feedback and the establishment of heart rate encoding

The final stage in the developmental sequence is the engagement of vagal sensory neurons (VSNs). By 12 dpf, sensory axons reliably reach the sinoatrial and atrioventricular plexuses, but encoding analyses revealed that VSNs can track heart-rate fluctuations even at earlier ages, likely through indirect coupling via vasculature, gill innervation or innervation of other tissues with heart rate-coupled dynamics (Fig. 6-8). Across spontaneous fluctuations, responses to visual threat, local optogenetic perturbations, and optogenetically-induced tachycardia, a single biphasic kernel accounted for neural responses. Neurons differed only in phase, tiling the kernel at systematic offsets. This phase-shifted organization sharpened with development, evolving into a bimodal distribution centered at opposite polarities. Such a population code allows continuous coverage of heart-rate dynamics while remaining robust to the source of modulation. The emergence of this stimulusinvariant representation demonstrates that VSNs contribute a stable and generalizable code for cardiac interoception.

We term this computational property of heart rate encoding, chronoreception by analogy to chemoreception and baroreception, and consistent with the use of the word “chronotropic” in reference to changes in heart rate [35]. To our knowledge, such heart rate-responsive cells have not been previously described. Hemodynamically, it makes sense to be able to ascertain heart rate directly because, along with contractility, it is one of the two variables that determine cardiac output. Regardless of whether heart rate is detected at the heart or in other structures, the heart rate encoders are present across conditions and are inputs to the brain.

Integration of efferent and afferent development into a staged program

Together, these results reveal a temporally ordered program for assembling heart–brain circuits. First, descending motor vagus fibers establish a minimal efferent control structure with a hindbrain premotor nucleus providing inhibitory drive (Fig. 9, 5dpf). Second, sympathetic inputs add modulatory power to gain and dynamics (Fig. 9, 7dpf). Third, vagal sensory neurons innervate the heart, they implement a distributed, phaseshifted code for heart rate feedback even before stable heart innervation (Fig. 9, 12dpf). Each stage introduces new computational capabilities without disrupting earlier ones, allowing larvae to maintain functional control as circuit elements mature. This progressive layering provides a framework for understanding how regulation of internal organs is assembled during development.

Figure 9.

Comparisons to mammalian autonomic development

The sequence we describe parallels and opportunities for the study of autonomic development in mammals. In rodents and humans, autonomic innervation of the heart arises prenatally, making it difficult to observe the onset of function in vivo [7]. Fetal studies suggest that vagal tone emerges early and that sympathetic influences strengthen later in gestation [36]. Sensory encoding of blood pressure has been inferred from nodose baroreceptor recordings in mammals, which respond to beat-to-beat pressure oscillations [37, 38], but direct demonstration of heart rate dynamic encoding has been lacking. Our results provide the first continuous, cell-resolved account of how HR is encoded in a vertebrate vagal sensory ganglion.

Implications for interoception and autonomic disorders

Recent studies have shown that interoceptive signals from the heart influence emotion, cognition [39], and behavior [40], and disruptions in cardiac sensory coding have been implicated in conditions ranging from anxiety to dysautonomia [41, 42]. The discovery of a canonical, phase-shifted kernel code for heart rate provides a mechanistic basis for how the brain may reconstruct cardiac dynamics from distributed afferents. In disorders such as postural orthostatic tachycardia syndrome [43, 44] or long-COVID dysautonomia [45], where feedback control of heart rate is impaired, alterations in this encoding scheme may contribute to pathophysiology. Similarly, aberrant development of sympathetic motifs or premotor inhibitory control could destabilize cardiac regulation. By establishing zebrafish as a model in which all three limbs of the autonomic circuit can be interrogated optically, our work opens the door to dissecting how genetic, developmental, or environmental perturbations disrupt interoceptive coding.

Conclusion

We have shown that the vertebrate nervous system assembles neural circuits for heart–brain regulation through a staged sequence that first establishes descending control, then incorporates sympathetic modulation, and finally engages sensory nerves. Each stage introduces distinct computational motifs: proportional kernels in the motor vagus, proportional-integral-derivative motifs in the sympathetic ganglia, and phase-shifted population codes in vagal sensory neurons. This architecture ensures that cardiac control is robust, adaptable, and ready to support survival as soon as larvae begin to interact with their environment. More broadly, these findings provide a framework for probing the developmental logic of interoceptive circuits and for exploring how their disruption contributes to autonomic disease.

Materials and methods

Key resources table

2. Table with the age and number of animals per experiment

Here we report the number of animals used for results that do not explicitly have the number of animals in the main figure.

Experimental model and subject details

Zebrafish housing and maintenance

For all experiments, we raised small groups of 20–30 larvae in filtered fish water (exchanged daily) in Petri dishes (diameter: 9 cm) at a constant 28 °C. After 5 dpf, larval fish were fed with paramecia once per day, after 9dpf twice per day. We performed experiments with larvae at ages 2-12 dpf for developmental experiments. All animals were kept in accordance with The Harvard Institutional Animal Care and Use Committee (IACUC) [54]. The sex of the fish tested was unknown because it is only distinguishable after 20 dpf [55].

Transgenesis

Transgenic zebrafish larvae were in casper or nacre background [56]. The fish (UAS:GCaMP6f)^jf46 and (UAS:CoChR-eGFP)^jf44 are described in Mu et al.[52]. The experiments that did not involve calcium imaging were conducted with Nacre fish. Anatomical observations were conducted using the progeny of Tg(chata:Gal4VP16, UAS:GFP)^mpn202 and Tg2(phox2bb: eGFP) for Figures 1 and 2 and the respective supplementary figures.

For calcium imaging, we used transgenic zebrafish lines expressing:

• Cytosolic GCaMP6s in glutamatergic neurons, Tg(vglut2-Gal4, UAS:GCaMP6s) (Filippi et al.[49]) (Fig. 7,8)

• Cytosolic GCaMP6s in cholinergic neurons, Tg(chata:Gal4VP16, UAS:GCaMP6s)^mpn202 (Förster et al.[46]) (Fig. 3)

• Cytosolic GCaMP6f in cells that express Tyrosine Hydroxylase, Tg(ZK10205-TH1:Gal4-KI, UAS:GCaMP7f) (Kim et al. [50] and Mu et al. [52]) (Fig. 5)

For optogenetic experiments, we used transgenic zebrafish lines expressing:

• CoChR in Tyrosine Hydroxylase cells Tg(ZK10205-TH1:Gal4-KI, UAS:CoChR-eGFP) (JFRC) (Fig. 5)

• CoChR in cholinergic neurons, Tg(chata:Gal4VP16, UAS:CoChR-tdTomato) (Fig. 3)

• CoChR in cardiomyocites Tg(vglut2-Gal4, UAS:GCaMP6s, myl7-UAS-CoChR-tdTomato) (this study) (Fig. 8)

• For mosaic optogenetic experiments, we used the injected embryos (f0) from vglut2-Gal4, UAS:GCaMP6s, myl7-UAS-CoChR-tdTomato.

Experimental Methods

Immunohistochemistry

For myl7 and acetylated-alpha tubulin staining, animals were fixed in 4% PFA (Sigma-Aldrich) in PBS (Gibco, Grand Island, NY, USA). Primary antibodies were diluted in 3% bovine serum albumin in PBS and incubated at 4°C overnight. Secondary antibodies were diluted in 0.1% Triton in PBS and incubated for 2 h at room temperature. Alexa Fluor 488 was used for acetylated-alpha tubulin or Tyrosine Hydroxylase and Alexa Fluor 546 for myl7.

Drug treatment

Propranolol treatments for Supp. Fig. 5 were systemic, with fish being placed in filtered fish water solutions of propranolol for 2 hours before measuring their responses to optogenetic stimulation of the APG. We used propranolol from Sigma-Aldrich (P0884). The concentration used was 20 µM of Propranolol (Supp. Fig. 5). One hour before the hemodynamic measurements, fish were partially immobilized using 2% agarose gel (Invitrogen, 16520). After 5 minutes of agarose solidifying, the tails of fish were released using a feather micro ophthalmic scalpel (Fisher Scientific link). Partially restrained fish were also in filtered fish water solutions of propranolol after embedding and during the heart rate measurements. The drug was washed out and the fish were tested for the control measurements 3-5 hours after the initial measurement. The concentrations of propranolol was selected based on the reported concentrations used in [57].

Anatomical imaging

For fixated samples

Fixated samples were mounted on bottom glass Petri dishes (Thermo Fisher Scientific 150680) using 1.2% agarose and positioned using plastic tips. Mounted samples were then imaged using a Zeiss LSM 900 inverted scanning confocal microscope. Images were later processed in Zen (Zeiss microscopy software) and ImageJ to reconstruct the volumes.

For live animals

Live zebrafish were embedded in 1.8% low melting point agarose between 1 and 2 hours before imaging. They were positioned laterally using plastic tips on a bottom glass Petri dish. The fish were then imaged using a Zeiss LSM 900 inverted scanning confocal microscope. The step size and speed of recording were not set to the optimal confocal setting as we did for fixated samples. Instead, we selected a speed that allowed visualization of neurites that move with the heart as it beats (Supp. Vids. 1, 2, 3). Before imaging, the most fluorescent fish were selected by screening them using an Olympus MVX10 epifluorescence microscope. For fish that were imaged on both sides, fish were recovered after the first round of imaging using a micro ophthalmic scalpel to release the fish and then embedded on the other side after 1-2 hours.

Serial section electron microscopy

The following protocol was specifically modified from Tapia et al.[58] to enhance extra-cellular space preservation, which improves synapse detection [59]. Unless otherwise noted, all steps were performed at room temperature (RT). A 7 dpf larva was anesthetized and embedded in agarose. The water was replaced with a dissection solution (64 mM NaCl, 2.9 mM KCl, 10 mM HEPES, 10 mM glucose, 164 mM sucrose, 1.2 mM MgCl2, 2.1 mM CaCl2, pH 7.5) supplemented with 0.02% tricaine3. Small slits were cut in the agarose to expose the eyes, and bilateral enucleations were performed to enhance ultrastructural preservation of the extra-cellular space and improve heavy metal staining. A custom-made hook was carefully inserted behind the eyes to minimize brain damage. The larva was then immediately transferred to a cold (4°C) fixation solution composed of 2.5% glutaraldehyde in 0.1 M cacodylate buffer (pH 7.4) supplemented with 4.0% mannitol. The cacodylate buffer was prepared with 0.3 M sodium cacodylate and 6 mM CaCl2, adjusted to pH 7.4. To improve fixation, the tissue was rapidly microwaved in the fixative solution using a microwave system (cat. no. 36700, Ted Pella) equipped with a power controller, steady-temperature water recirculator, and cold spot. The microwaving protocol, based on Tapia et al. (2012)⁵⁷, was performed as follows: at power level 1 (100 W) for 1 min on, 1 min off, 1 min on, followed by power level 3 (300 W) for 20 s on, 20 s off, 20 s on, repeated three times. Fixation was continued overnight at 4°C in the same solution. The following day, the sample was washed three times (30 min each) in 0.5x cacodylate buffer before osmication with 2% OsO4 in 0.5x cacodylate buffer for 90 minutes. After a brief wash in 0.5x cacodylate buffer (<1 min), the sample was incubated in 2.5% potassium ferrocyanide in 0.5x cacodylate buffer for 90 minutes. The sample was then washed with filtered water (three exchanges, 30 min each) before being incubated in 1% (w/v) thiocarbohydrazide (TCH) in filtered water for 45 minutes. Due to poor dissolution of TCH, the solution was preheated at 60°C for ~90 minutes with occasional shaking, then cooled to RT for 5 minutes prior to incubation. After incubation, the sample was washed with filtered water (three exchanges, 30 min each). A second osmication was performed in 2% OsO4 in filtered water for 90 minutes, followed by three additional washes (30 min each). En-bloc staining was then carried out overnight at 4°C using 1% uranyl acetate in filtered water. The uranyl acetate solution was sonicated for 90 minutes and filtered with a 0.22 μm syringe filter before use. All steps involving uranyl acetate were performed in the dark. The next day, the sample was washed with filtered water (three exchanges, 30 min each) and dehydrated through a graded ethanol series (25%, 50%, 75%, 90%, 100%, and 100% ethanol, 10 minutes each step). The sample was then transferred to 100% propylene oxide (PO) (two exchanges, 30 minutes each). Infiltration was performed with LX112 epoxy resin (Ladd, 21212) mixed with PO in a series of graded steps (25% resin/75% PO, 50% resin/50% PO, 75% resin/25% PO, and 100% resin, with each step lasting 4 hours). The sample was mounted in fresh resin in a mouse brain support tissue [60], with the head exposed to facilitate cutting. The mouse tissue, fixed using standard procedures [61], was cut into 2–3 mm cubes, which were pierced with a 0.75 mm puncher (EMS, 57395) to insert the larva. The cubes were stained alongside the fish samples using the protocol described above, except that the uranyl acetate overnight step was performed at RT. The samples were then cured with support tissue for 3 days at 60°C. During all steps, a rotator was used. Aqueous solutions were prepared with water purified through a filtration system (Arium 611VF, Sartorius Stedim Biotech). The entire protocol, including surgery, fixation, staining, and resin embedding, was completed in 5 consecutive days, followed by 3 days of resin curing. The cured block then was trimmed into a diamond shape. Sections (30–35 nm thick) were automatically collected on carbon-coated tape using a custom tape collection device (ATUM) mounted to a commercial ultramicrotome [62]. The tape was cut into strips and deposited on 30 silicon wafers, which were post-stained as previously described [60]. Wafers were then mounted on the 61-beam multiSEM (MultiSEM 505, ZEISS) stage, and the position of each section was determined using a reflected light microscope to guide high-resolution imaging. Imaging was performed at 4×4 nm pixel resolution using secondary electron emission with a dwell time of 400 ns per pixel. The quality of each individual section was assessed using previously described methods [63]. Cells within a 1 µm-thick slab, positioned posterior to and partially encompassing a cardiac neuron, were manually annotated using VAST8. Note that this sample was identical to Boulanger-Weill et al.[64], in preparation.

Custom setup for hemodynamic measurements

Functional imaging for heart rate during looming was conducted using a custom-built optical path based on the Navitar 6000 zoom system. We used a compounded magnification of 8x. To avoid visual artifacts, we used a near-infrared (NIR) ring of LEDs at 750nm and a 750nm band-pass filter (Semrock, 88-013). We acquired 84 frames per second using a Grasshopper 3.0 NIR camera (Teledyne Flir Model: GS3-U3-41C6NIR-C: CMOSIS CMV4000-3E12, NIR). The illumination LEDs were driven using pulse width modulation at 10% duty cycle to avoid heating the sample. The fish was embedded in 2% agarose and then the tail was released using a micro scalpel. This is the same embedding protocol described above under “Drug treatment”. Fish were embedded 1 hour before recordings. Once a fish was embedded, it was placed over a small opening of filter paper that was used to project images from a mini projector. Looming stimuli were designed to be a dot that increased in size and became darker following the derivative of an arc tangent function, which is the equivalent of assuming a predator approaches at constant acceleration. The maximum size of the dot was reached in 6 seconds, and after that, the fish remained in darkness for another 12 seconds before the background light appeared again. To prevent adaptation, fish were exposed to only one stimulus every 5 minutes. Stimuli were repeated 5 times per fish. Similarly, for dark flash experiments background red light was turned off for 32 seconds every five minutes. Tail flicks were tracked using custom Python code previously presented in [65].

Heart rate quantification

To calculate heart rate, we first calculated the “motion energy” (similar to rodent behavior estimation Stringer et al [66]), which is the per-pixel brightness derivative over time. For the image I which consists of the pixels of one cardiac chamber (atrium or ventricle), the motion energy M_E would be . We then calculated the squared root of the change in brightness in M_E to identify the heart rate trace HR_T, which was the subject to an inverse Fourier Transform to calculate the spectrum of the heart rate trace . The most likely frequency band in a small vicinity of a time point t is the heart rate HR(t). When this method failed, we used a custom script to detect the distances between local peaks in the denoised and detrended contraction wave, implemented in MATLAB.

Functional imaging

Similar to the protocol described for embedding fish in the Heart Rate Quantification portion of this Methods section, fish were transferred into the measurement rig after 2 hours of embedding. We used a custom-built two-photon microscope, operated by custom-written Python 3.7-based software [65]. A femtosecond-pulsed MaiTai Ti:Sapphire laser (Spectra Physics) tuned to 920 nm was used to image for GCaMP6s/f. A set of x/ygalvanometers (Cambridge Technology), and a 20× Olympus infrared-optimized objective (XLUMPLFLN) were used to scan over the brain. We collected fluorophore emission using a green photomultiplier tube (PMT), amplified by a current preamplifiers (Stanford SR570). We used frame acquisition rates of around 1.5 Hz in brain imaging experiments. We adjusted laser power to 12.7 mW at the specimen, a low enough value that did not seem to interfere with the physiology of larvae. In our brainstem imaging experiments (Fig. 3), we imaged each plane at a spatial resolution of 0.35 μm per pixel (700 × 700 pixels) for 25 min, while presenting ~5 trials of 12 s of stimuli interleaved by 5 minutes of rest. We then moved the objective 10 μm to the next imaging plane and repeated the procedure. We acquired 4-6 planes, resulting in a total imaging time of around 1.5 h per fish.

In the case of vagal and sympathetic ganglia, we needed the sagittal view to visualize the entire ganglia, and 2-photon microscopy in this position resulted in slow overheating given the pigmentation in that region. Because these ganglia are easily accessible with a single-photon microscope from the sagittal view, we used a Zeiss Multiphoton LSM 980 to image these fish. Because the imaging region needed for these ganglia is smaller, we conducted functional imaging at 16.6 Hz. The scanning laser was restricted to the ganglia and therefore did not result in visual artifcats. In these experiments, UV flash at 410nm was used as a visual stressor.

Optogenetic stimulation

Optogenetic stimulation experiments of the motor limbs of the autonomic system were conducted using galvo mirrors to confine the illumination region of a 488nm laser to small subsets of neurons. The specificity of the activated region depends on (1) the sparseness of the genetic marker, and (2) the targeted region with the laser. Given the precision of the galvo mirrors used to direct the laser the tolerance for spatial illumination is around ~1.2 µm. The stimulation pulse lasted 1.2 seconds at 75.6% laser power with a 488nm laser and was repeated every 5 minutes. The preparation of larvae before the experiment was analogous to the functional imaging preparation. In the cases of the sympathetic system and cardiomyocytes, heart rate was measured before and after stimulation, but not during stimulation, because of the location of the stimulus pulse. Each fish was tested for 5 trials. Frames were acquired every 60 ms, and during the 1.2s of optogenetic stimulation, the stimulated region was scanned 20 times. Fish were screened for tdTomato (which is co-expressed with CoChR in our transgenic fish) at least 2 hours before the experiments.

Computation and Modeling Methods

This section details preprocessing, generalized linear models (GLMs), kernel estimation, cross-validation, dimensionality reduction for kernel shape, circular phase statistics, clustering, and quantification of optogenetic effects used throughout the study. Unless otherwise stated, analyses were performed in MATLAB and Python (NumPy/SciPy), with custom code available upon publication on our GitHub page.

Signals, Notation, and Preprocessing

Let HR(t) denote the heart-rate time series (Hz) sampled at interval Δt; let r_i(t) be the calcium signal (ΔF/Fo) for neuron i. Signals were detrended (second-order polynomial) and z-scored within fish/session.

Causal FIR Design Matrices

For models using discrete time lags, we constructed a lag matrix of HR with causal finite-impulse-response (FIR) supports. For a maximum lag L and step δ, the causal HR design matrix X was constructed with columns:

GLM for Sensory Encoding (HR → Neural Activity)

Neural activity was modeled as a convolution of HR with a causal kernel k_i(τ), plus a bias b_i.

We applied a static nonlinearity for saturation effects:

Kernels were estimated by ridge regression:

With D=I unless otherwise stated, λ was chosen by inner cross-validation.

GLM for Motor Decoding (Neural → HR)

For decoding HR from neural populations:

With g() a sigmoid as in eq. (5).

Cross-Validation and Performance

Blocked 5-fold cross-validation: 80% train, 20% test per fold. Model performance was quantified as crossvalidated R²:

Kernel Normalization and Shape Metrics

Kernels were L²-normalized for visualization in heatmaps:

Kernel PCA and Phase Extraction

Normalized kernels were projected into kernel principal component analysis (PCA) space. Phases were defined from the first two PCs:

Circular Statistics and Permutation Tests

For a set of phases {ϕ_i}, the circular mean vector was:

Permutation tests compared observed R to a null distribution obtained by shuffling or circularly shifting predictors:

Here the “#{·}” notation means “the number of elements satisfying the condition.” In words:

• R_obs is the resultant length of the actual data.

• R_b are resultant lengths from each of B permutations.

• #{b:R_b≥R_obs} counts how many permutation values were greater than or equal to the observed one.

This is the standard way of estimating a permutation p-value (with the +1 in numerator/denominator for finitesample correction).

Determining the number of clusters of kernels

To determine the optimal number of kernel clusters (K) for Fig. 5, we applied three complementary criteria:

1. Silhouette analysis

   • We computed the average silhouette coefficient across neurons for a range of candidate cluster numbers (K = 2–8).

   • The silhouette score measures how similar each neuron’s kernel is to others in its assigned cluster compared to neighboring clusters.

2. Calinski–Harabasz index (Variance Ratio Criterion)

   • This criterion quantifies the ratio of between-cluster variance to within-cluster variance.

   • Higher values indicate that clusters are compact and well separated.

3. Gap statistic

   • The gap statistic compares the within-cluster dispersion of the observed data to that of a null distribution generated from randomized (reference) datasets.

   • The optimal K is the smallest number of clusters such that the observed within-cluster dispersion is substantially smaller than expected under the null.

   • This method guards against overfitting to spurious structure.

All three metrics were computed for candidate values of K, and the chosen solution (K = 6) was supported consistently across methods (Supplementary Fig. 4). We further confirmed stability of the clustering by bootstrap resampling (subsampling neurons and repeating PCA+k-means), which yielded reproducible assignments.

Data availability

All physiological, calcium imaging, and optogenetic data will be deposited to DANDI (Distributed Archives for Neurophysiology Data) and all code will be made available on our lab GitHub page (https://github.com/BrainBodyInteractionsLab/) .

Acknowledgements

We thank members of the Hernandez-Nunez Lab for fruitful discussions and suggestions, in particular Shulin Zhang, Aliya Ablitip, Renzo Mendoza, Ting Ting Yan, and Sofia Melnychuck. We also want to thank members of the Engert and Fishman Labs: Vickie Wang, Yasuko Isoe, Kristian Herrera, Soma Singareddy, and Erin Song for early contributions and for useful discussions and feedback on this project. We would like to thank Richard Schalek and Juan-Carlos Tapia for their help with the serial section electron microscopy sample preparation and imaging. We thank the Harvard Center for Biological Imaging (RRID:SCR_018673) for infrastructure and support.

Luis Hernandez-Nunez received funding from the Harvard Mind Brain and Behavior Fellowship and Young Investigator Award, by the Life Sciences Research Foundation/Additional Ventures Fellowship, by the Warren Alpert Distinguished Scholar Award, by the Branco Weiss Fellowship, by the Burroughs Wellcome CASI grant, and by the National Institutes of Health (R34NS138096). Florian Engert received funding from the National Institutes of Health (U19NS104653 and R01NS124017), and the Simons Foundation (SCGB 542973 and NCGB-CULM-00003241-02). Mark Fishman received funding from the National Institutes of Health (U19NS104653 and R01NS124017).

Additional information

Author Contributions

L.H.-N., F. E., and M. C. F. designed the study, interpreted the results, and wrote the manuscript with feedback from all authors. L.H.-N. developed the instruments and experimental techniques; wrote the software for data analysis; analyzed data; performed anatomical, behavioral, physiological, functional imaging, and optogenetic experiments; and built the encoding and decoding models. J.A. performed physiological, functional imaging, and optogenetic experiments related to the Sensory Vagus System (Fig. 7 and 8) with guidance from L.H.-N. S.S. performed physiological, functional imaging, and optogenetic experiments related to the Motor Vagus System (Fig. 2 and 3) with guidance from L.H.-N. A. M. performed physiological, functional imaging, and optogenetic experiments related to the Sympathetic System (Fig. 4, 5) with guidance from L.H.-N. A.K. conducted anatomical experiments of the Sensory Vagus Nerve that were the initial step for the results of Fig. 6 with guidance from L.H.-N. J. B.-W. contributed the electron microscopy reconstruction of the sinoatrial plexus of 7 dpf fish (Fig. 2G-I). V.R. and M.A. contributed early results on optogenetic activation of the motor vagus nerve, which were the initial steps toward Fig. 3R and S. A. Z.-S. generated the transgenic fish used for Fig. 8.

Nomenclature:

• A:: Atrium

• V:: Ventricle

• BA:: Bulbous Arteriosus

• SAP:: SinoAtrial Plexus

• AVP:: AtrioVentricular Plexus

• ICN:: IntraCardiac Neuron

• ICG:: IntraCardiac Ganglion

Funding

National Institutes of Health (R34NS138096)

National Institutes of Health (U19NS104653)

National Institutes of Health (R01NS124017)

Warren Alpert Foundation (Distinguished Scholar Award)

• Luis Hernandez-Nunez

Burroughs Wellcome Fund (Career Award at the Scientific Interface)

• Luis Hernandez-Nunez

Branco Weiss, Society in Science (Branco Weiss Fellowship)

• Luis Hernandez-Nunez

Life Sciences Research Foundation (LSRF Fellowship)

• Luis Hernandez-Nunez

Additional files

Supplementary video 1. Example of a 5 dpf phox2bb-GFP zebrafish. We identify neurites in the SAP.

Supplementary video 2. Example of a 7 dpf phox2bb-GFP zebrafish. We identify the first ICN.

Supplementary video 3. Example of a 12 dpf phox2bb-GFP zebrafish. At this age, the ICG has appeared in the atrium, and more neurons are in the SAP and AVP.

Supplementary video 4. Example of a 4 dpf ChaTA-Gal4, UAS-GFP larval zebrafish. At this age the cardiac branch of the motor vagus nerve is not yet developed, thus there are no neurites in the recording of the beating heart.

Supplementary video 5. Example of a 5 dpf ChaTA-Gal4, UAS-GFP larval zebrafish. At this age the cardiac branch of the motor vagus nerve has reached the heart and it can be seen moving with the SAP during atrium contractions (yellow arrow in the video).

Supplementary video 6. Example of a 7 dpf ChaTA-Gal4, UAS-GFP zebrafish. At this age the first ICNS neuron has reached the heart, and it moves with atrium contractions.

Supplementary video 7. Example of a 12 dpf ChaTA-Gal4, UAS-GFP zebrafish. At this age the neurites have expanded to reach the AVP (arrow in the video).