Developmental tracking of the intracardiac nervous system reveals a minimal structure for neural control of the heart.

(A) Schematic representation of a custom-built experimental setup for imaging the visceral cavity of zebrafish larvae, including the heart. Visual stimuli can be delivered with a projector or an LED. (B) HR extraction. Top: example frames during ventricular dilation and atrial contraction; bottom: NIR reflectance traces from A (green) and V (cyan) used to compute beat timing (1 s scale bar). (C) HR responses to UV flash across development (4, 5, 7, 12 dpf). Thick black line, mean normalized HR; gray shading, ± SE; vertical violet bars, stimulus onset; n=14–20 fish per age. (D) HR responses to dark flash across the same ages. Plotting as in (C) with gray stimulus bars. (E) Cardiac innervation timeline in fixed hearts stained for acetylated α-tubulin (neurites, cyan) and Myl7 (myocardium, magenta). Innervation appears at SAP by 5 dpf, extends through the atrium by 7 dpf, and reaches AVP and ventricle by 12 dpf. Scale bars, 50 µm. (F) In vivo phox2bb:GFP imaging of the heart region. 4–5 dpf: axons approach the heart but no ICNs are present; 7 dpf: first ICN appears; 12 dpf: multiple ICNs (magenta arrowheads) are evident. Scale bars, 50 µm. (G) Immunohistochemistry for Phox2bb (cell bodies, green) and Myl7 (magenta). 5 dpf: no ICN somata in the heart; 12 dpf: ICNs near SAP (≈6), near AVP (≈6), and a dorsal-posterior atrial cluster (ICG; ≈16). Scale bars, 50 µm. (H) Schematic summary of the developmental sequence: 5 dpf—axon entry at SAP; 7 dpf—first ICN; 12 dpf— ICNs at SAP/AVP and ICG. BA, bulbus arteriosus. (I) Population summary of innervation phenotypes by age (proportion of fish; 4 dpf n=12, 5 dpf n=15, 7 dpf n=20, 12 dpf n=12). Categories: no innervation; innervation only; ≥1 ICN; multiple ICNs. Conventions. HR traces are normalized per fish. SE shading indicates standard error (between-fish variability). All scale bars, 50 µm. Abbreviations: NIR, near-infrared; SAP, sinoatrial plexus; AVP, atrioventricular plexus; ICN, intracardiac neuron; ICG, intracardiac ganglion; BA, bulbus arteriosus.

Formation of a direct cholinergic brain-to-heart motor circuit.

(A–C) Immunostaining for ChAT (green; cholinergic fibers/neurons) and Myl7 (magenta; myocardium) at 5, 7, 12 dpf. Cholinergic fibers reach the sinoatrial region by 5 dpf; the first ChAT-positive ICN appears at 7 dpf; multiple ICNs cluster near the SAP by 12 dpf. Scale bars, 50 µm. (D) In vivo volumetric imaging in ChaTA-Gal4>UAS-GFP larvae with transmitted-light myocardium. Left→right: 4 dpf (no innervation), 5 dpf (emerging fibers), 7 dpf (first ICN), 12 dpf (extension toward AVP). SAP/AVP and the atrium (A) and ventricle (V) are outlined; arrowheads mark neurites/ICN. Scale bars, 100 µm. (E) Schematic progression of cholinergic cardiac innervation across ages: fibers first contact the SAP (5 dpf), a single ICN appears (7 dpf), then 2–6 ICNs are present (12 dpf). BA, bulbus arteriosus. (F) Population summary of cholinergic phenotypes by age: no innervation (gray), innervation only (light teal), one ICN (teal), 2–6 ICNs (dark teal). 4 dpf n=15n=15n=15; 5 dpf n=14n=14n=14; 7 dpf n=17n=17n=17; 12 dpf n=14n=14n=14. (G) ssEM reconstruction of the SAP at 7 dpf with segmented cell types (legend). Neurites approaching the SAP are indicated (arrows). Scale bar, 2 µm. (H) EM micrograph showing a vesicle-rich axonal terminal contacting the ICN (arrowheads). (I) EM micrograph showing a synaptic contact from the ICN onto a cardiomyocyte (arrowheads). Notes. Cholinergic identity of the first ICN was corroborated by HCR-FISH for ChAT and overlap with ChaTAGal4>GFP and phox2bb:GFP (Supp. Fig. 2). No ChAT-positive ventricular innervation was detected through 12 dpf. Acronyms: A, atrium; V, ventricle; SAP, sinoatrial plexus; AVP, atrioventricular plexus; ICN, intracardiac neuron; BA, bulbus arteriosus.

Motor-vagus decoding of heart rate and emergence of a vagal premotor nucleus.

(A) Experimental schematic. Calcium imaging of hindbrain cholinergic neurons (ChaTA-Gal4), simultaneous optical heart-rate (HR) readout, and dark-flash stimuli. (B–E) Example recordings across development (4, 5, 7, 12 dpf). Top: HR (Hz). Bottom: representative neuronal ΔF/F0 heatmaps (rows, neurons; columns, time). (F) Decoding kernels for all predictive neurons (exceeding Otsu threshold on cross-validated performance) plotted as a heatmap and ordered by shape similarity. Kernels are causal finite-impulse responses spanning 0–4 s in 0.25 s steps. (G) Canonical kernel motifs obtained by clustering all predictive kernels: C1 (positive, red) and C2 (negative, blue). Curves show means with shaded standard error (SE); n values indicate the number of neurons assigned to each motif. (H) Decoding model. Each neuron contributes a causal linear kernel; outputs are summed and passed through a static sigmoid nonlinearity to capture HR saturation. (I) Model performance by age using blocked 5-fold cross-validation (80% train, 20% test per fold). Points, perfish CV R2; violins, distributions; horizontal bars, medians. Annotated p values are Wilcoxon rank-sum with BHFDR correction (Methods). Across 31 fish (5, 7, 12 dpf), only 2 at 5 dpf and 2 at 7 dpf fell below CV R2 = 0.20; all others achieved CV R2 between 0.25 and 0.60. (J) Example predictions for two fish near the cohort median (CV R2 = 0.35 and 0.43). Predicted HR closely tracks measured HR and stimulus-locked responses. (K) Left: fraction of non-encoding and encoding neurons per age (encoders <5% at each age; prevalence increases at 12 dpf relative to 5 dpf). Right: within encoders, fraction assigned to C1 (positive) or C2 (negative). The negative motif expands from ~40% at 5 dpf to ~90% at 12 dpf, consistent with the cardioinhibitory role of motor vagus. (L–N) Anatomy maps showing the spatial distribution of encoders at 5, 7, and 12 dpf. Red, C1 (positive); blue, C2 (negative). Encoders cluster medially across ages, with positive encoders more anterior and dispersed, and negative encoders concentrated near the midline. (O) Summary schematic of encoder locations across fish: reproducible medial cluster of negative encoders and broader, less consistent distribution of positive encoders. (P) Optogenetic testing strategy. ChaTA-Gal4 drives CoChR in cholinergic hindbrain neurons; patterned illumination targets the negative-kernel locus; HR recorded simultaneously. (Q) Targeted stimulation region for the premotor locus (example plane). (R) At 4 dpf, before anatomical heart innervation is detectable, photoactivation does not measurably change HR. (S) At 7 dpf, after innervation appears, the same stimulation produces a robust bradycardia. Traces show mean ± SE; vertical bars mark stimulus onset.

Development of the sympathetic cardiac circuit.

(A) Paravertebral sympathetic neurons in vivo. Arrowheads mark sympathetic neurons along the trunk; APG indicates the anterior paravertebral ganglia that project toward the heart (schematic, right). Scale bar, 100 µm. (B) APG morphology over development (5, 7, 12 dpf). TH immunolabeling reveals growth and the characteristic C-shaped configuration. Scale bars, 50 µm. (C) Sympathetic cardiac innervation (TH, green) with myocardium (Myl7, magenta). 5 dpf: no cardiac TH fibers. 7 dpf: TH-positive projections reach the SAP (arrowheads). 12 dpf: SAP innervation strengthened; no AVP or ventricular TH fibers detected. Scale bars, 50 µm. (D) Schematic summary of innervation patterns at 5, 7, and 12 dpf. BA, bulbus arteriosus; A, atrium; V, ventricle; SAP, sinoatrial plexus; AVP, atrioventricular plexus. (E) Population summary of sympathetic cardiac innervation by age (with/without TH-positive fibers at the SAP). 5 dpf n=9: 0% innervated; 7 dpf n=13: 46% innervated; 12 dpf n=10: 100% innervated. Notes. TH marks catecholaminergic sympathetic fibers; all panels use TH (green) and Myl7 (magenta) unless indicated. Innervation remained SAP-restricted up to 12 dpf in this cohort.

Emergence of control motifs and function in sympathetic cardiac circuits.

(A) Assay schematic. TH1-Gal4>UAS-GCaMP6s labels sympathetic neurons in the anterior paravertebral ganglia (APG). Calcium activity and heart rate (HR) are recorded simultaneously during brief UV flashes. (B–D) Representative recordings at 5, 7, and 12 dpf. Top: HR (Hz). Bottom: APG population activity (ΔF/F heatmaps; rows = neurons; columns = time). Vertical purple bars indicate UV flashes. (E) Linear kernels for all predictive neurons, pooled across fish and ordered by shape similarity. Lags: 0–4 s in 0.1-s steps; color denotes normalized kernel amplitude. (F) Canonical kernel motifs (cluster means ± SE). Unsupervised clustering reproducibly yielded six motifs: two biphasic/derivative-like (C1, C4), two monophasic/proportional-like (C2, C5), and two monotonic ramp/integrator-like (C3, C6). n per motif indicated. (G) Decoding GLM diagram. Each neuron contributes a causal kernel; outputs are summed and passed through a static sigmoid nonlinearity to predict HR. (H) Cross-validated performance by age (blocked 5-fold; 80% train/20% test per fold). Points, fish; violins, distributions; bars, medians. P values (above) are Wilcoxon rank-sum with BH-FDR correction. (I) Example held-out predictions from the GLM for two fish (CV R2 = 0.45 and 0.56). Blue, measured z(HR); orange dashed, model prediction. (J) Cluster composition across development. Stacked bars show the proportion of neurons assigned to each motif (C1–C6) per age. Derivative-like motifs are rare at 5 dpf and emerge by 7–12 dpf. (K–M) Spatial organization of APG neurons colored by cluster identity in left and right ganglia at 5, 7, 12 dpf. Dorsal–ventral and anterior–posterior axes indicated; scale bars, 50 µm. Neuron positions vary between fish and between sides. (O) Optogenetic test of function. TH1-Gal4>UAS-CoChR-GFP neurons in the APG are activated with patterned light (galvo-steered) while HR is recorded. (P–R) Trial-averaged HR responses to APG activation at 5, 7, 12 dpf (gray, individual trials; red, mean; blue bar, stimulus). No effect at 5 dpf; increasing tachycardia at 7–12 dpf. Response is β-adrenergic-dependent (abolished by propranolol; see Supp. Fig. 5). Analysis details. Kernels estimated on the masked time base with causal fractional lags (0–4 s, 0.1 s step). Predictive neurons defined by Otsu threshold on per-neuron CV performance (with a 0.2 floor). GLM performance reported as CV R2 (blocked folds). Clusters derived by k-means clustering in kernel PCA space (see Methods and Supp. Material, Clustering Kernels); motif means plotted with standard error (SE).

Developmental establishment of the heart-to-brain interoceptive circuit.

A) Vagal sensory anatomy. Vglut->GFP labels cranial sensory neurons; the nodose ganglion (posterior) and epibranchial ganglia (anterior) are indicated (left). Schematic (right) shows ganglia positions and sensory projections toward the heart. Scale bar, 50 μm. B) In vivo imaging of vagal sensory projections to the heart at 5, 7, and 12 dpf. Dashed outlines mark atrium (A) and ventricle (V). 5 dpf: ganglia present but no cardiac sensory fibers. 7 dpf: sparse fibers in a subset of animals. 12 dpf: robust sensory projections reach the sinoatrial plexus (SAP) and atrioventricular plexus AVP) (arrowheads). Scale bars, 50 μm. C) Schematic summary by age: absence dpf, variable onset at 7 dpf, and established SAP/AVP innervation by 12 dpf. BA, bulbus arteriosus. D) Population summary of sensory vagus cardiac innervation (presence/absence of heart-projecting fibers) across development. 5 dpf n=8: 0% innervated; 7 dpf n=10: 20% innervated; 12 dpf n=11: 100% innervated. Acronyms: A, atrium; V, ventricle; SAP, sinoatrial plexus; AVP, atrioventricular plexus; BA, bulbus arteriosus.

Cardiac encoding in vagal sensory neurons.

(A) Experimental schematic. Calcium imaging of vagal sensory neurons (vglut2-Gal4>UAS-GCaMP6s) was performed simultaneously with optical heart-rate (HR) readout during UV-flash stimuli. (B–D) Representative recordings of HR (top, Hz) and vagal ganglion population activity (bottom, ΔF/F0 heatmaps; rows = neurons; columns = time) at 5, 7, and 12 dpf. Purple bars indicate stimuli. (E) Proportion of cardiac-encoding neurons across development. Encoding was defined as neurons with CV R2 > 0.2 in GLM fits of HR→neural activity. Points are individual fish, violins show distributions, and p values are Wilcoxon rank-sum with BH-FDR correction. (F) Heatmap of normalized linear kernels from all encoding neurons, sorted by shape similarity. A single biphasic motif is expressed with systematic phase shifts. (G) Principal component analysis (PCA) of kernel shapes. Top: distribution of kernel phases in PC1–PC2 space lies along a circular trajectory. Bottom: canonical kernels at four representative phase bins (0, π/2, π, 3π/2). (H) Schematic representation of phase-shifted kernel ensembles. Neurons tile the canonical waveform across four quadrature-like groups. (I) Model performance across development. Violin plots show CV R2 values (blocked 5-fold cross-validation) for left and right vagus at 5, 7, and 12 dpf. Horizontal bars denote medians; annotated p values are Wilcoxon ranksum with BH-FDR correction. (J) Example neurons illustrating predictive power of the GLM. Black, measured ΔF/F0; green, predicted ΔF/F0. Top: neuron with phase π kernel decreases activity when HR rises. Bottom: neuron with phase 0 kernel increases activity when HR rises. (K) Radial histograms of kernel phase distributions at 5, 7, and 12 dpf. Early neurons are biased toward 0 and π; by 12 dpf distributions become bimodal, with encoders at opposite phases. Inset: permutation p-value heatmap showing that 5 dpf distributions differ significantly from 7 and 12 dpf. (L–N) Spatial distribution of encoding neurons colored by kernel phase in left and right nodose ganglia at 5 (L), 7 (M), and 12 (N) dpf. Neurons of different phases are intermixed within each ganglion and across both sides.

Stimulus-invariant encoding of heart rate in vagal sensory neurons.

(A–C) Experimental paradigms for three sources of heart-rate modulation. (A) Spontaneous fluctuations in the absence of external perturbations. (B) Mosaic optogenetic activation of small subsets of cardiomyocytes (blue arrow, galvo-delivered stimulation). (C) Global optogenetic activation of the entire myocardium (blue arrow, galvo-delivered stimulation). Calcium imaging of vagal sensory neurons was performed simultaneously with optical readout of heart rate. (D–F) Representative recordings of heart rate (top, Hz) and vagal ganglion calcium activity (bottom, ΔF/F0 heatmaps; rows = neurons; columns = time) during (D) spontaneous fluctuations, (E) mosaic optogenetic perturbations, and (F) global optogenetic tachycardia. Blue bars mark stimulus onset in optogenetic conditions. (G–I) Heatmaps of causal kernels from GLMs trained to predict neural activity from HR dynamics under each condition. Kernels are aligned by phase, revealing the same biphasic motif expressed with systematic phase shifts, as observed under UV stimulation (Fig. 7F). (J–L) Principal component analysis of kernel shapes. Kernels distribute along a circular trajectory in PC1–PC2 space, consistent with phase-shifted versions of a single canonical filter. (M) Cross-validated model performance (CV R2) across conditions. Violin plots show distributions for left (L) and right (R) vagus at 12 dpf. Median encoding strength was comparable across conditions (0.4–0.55). (N) Canonical kernel shapes recovered across spontaneous, mosaic, and global optogenetic perturbations converge to the same biphasic temporal motif. UV stress condition from Fig. 7 is replotted for comparison. (O–Q) Anatomical maps of encoding neurons colored by kernel phase in left and right nodose ganglia under (O) spontaneous fluctuations, (P) mosaic optogenetic perturbations, and (Q) global heart optogenetic activation. Phase groups are intermixed within each ganglion and across both sides, with a bias toward ventral positions but no strict spatial segregation.

Emergence of Functional Heart-Brain Circuits.

Graphical summary of anatomical connections indicated by arrows and neural computations, indicated by kernel shapes across developmental stages. VPN: Vagal Premotor Nucleus. APG: Anterior Paravertebral Ganglia. VSNs: Vagal Sensory Neurons.

Developmental progression of heart rate responses to visual stimuli.

(A) Schematic of the experimental setup. Larval zebrafish were embedded with the tail released, and visual looming stimuli were projected from below, while near-infrared imaging was used to extract heart rate. (B) Baseline heart rate distributions across development (4, 5, 7, and 12 dpf). Each violin shows the distribution across fish, with median and interquartile ranges indicated (n = 15–20 fish per age). * indicate p<0.0001 and ** indicate p< 0.001. p-values from Steel-Dwass test. (C–F) Average heart rate responses to looming stimuli at (C) 4 dpf, (D) 5 dpf, (E) 7 dpf, and (F) 12 dpf. Thin gray traces represent individual animals, the thick black line is the mean across fish, and shaded regions denote the Standard Error (SE) of the mean. Vertical gray bars mark the stimulus presentation. (G) Representative simultaneous recordings of heart rate (red) and tail angle (cyan) in response to looming at 7 dpf (left) and 4 dpf (right). At both ages, tail flicks indicate behavioral responses, but only the 7 dpf fish show robust heart rate modulation. (H) Quantification of heart rate responses to UV flash, dark flash, and looming across development. Violin plots show normalized heart rate in baseline, peak, and adapted epochs for each fish (5 trials averaged per animal). Each dot represents one fish; horizontal bars denote medians and vertical bars denote the 90% interquartile range. p-values from the Steel–Dwass test are reported for each pairwise comparison.

Fluorescent in situ hybridization confirms cholinergic identity of first intracardiac neuron.

(A) Confocal volumetric imaging of Tg(chata:Gal4, UAS:GFP) larvae. GFP expression (green, 488 nm) labels cholinergic neurons, while hybridization chain reaction (HCR) probes against ChaTA (magenta, 546 nm) confirm transcript expression. Insets show higher-magnification views of a representative neuron (boxed). Merged images demonstrate overlap of GFP and ChaTA HCR signals within intracardiac neurons. (B) Confocal volumetric imaging of Tg(phox2bb:EGFP) larvae. GFP expression (green, 488 nm) labels visceral motor and autonomic neurons, while HCR probes against ChaTA (magenta, 546 nm) confirm cholinergic identity. Insets show higher-magnification views of a representative intracardiac neuron (boxed). Merged images again demonstrate colocalization of GFP and ChaTA HCR signals.

Control optogenetic stimulation of the anterior spinal cord does not affect heart rate.

(A) Schematic of the experimental setup. Larval zebrafish were embedded for simultaneous heart imaging and targeted optogenetic stimulation. Blue light stimulation was directed to the anterior spinal cord region as a control, outside of the vagal premotor population. (B) Confocal image showing the stimulation site (blue shading) relative to the spinal cord and motor vagus nerve. The stimulation region was restricted to the dorsal spinal cord and excluded premotor vagal neurons. Scale bar, 50 μm. (C) Heart rate responses at 7 dpf during control spinal cord stimulation. Thin gray traces represent individual trials on 6 fish, the thick black trace shows the mean across fish, and shading denotes the Standard Error of the Mean. No significant change in heart rate was observed, confirming that bradycardic effects in Fig. 3 arise specifically from vagal premotor activation.

Determining the number of kernel clusters.

(A) Silhouette analysis for candidate values of K. The silhouette score quantifies the similarity of each kernel to its assigned cluster relative to neighboring clusters. The mean silhouette across neurons is plotted; higher values indicate more coherent clusters. The maximum score was observed at K = 6. (B) Gap statistic as a function of K, computed with 30 bootstrap reference datasets (B = 30). The gap statistic measures the separation between observed within-cluster dispersion and that expected under a null distribution. The inflection point at K = 6 supports the chosen cluster number. Error bars denote ± s.e.m. across bootstrap samples. (C) Elbow criterion applied to explained variance by clustering. The ratio of between-cluster to total variance increases with K but shows diminishing returns beyond K = 6, indicating this as the optimal choice. Across all three metrics—the silhouette coefficient, gap statistic, and Calinski–Harabasz index—K = 6 was consistently supported as the optimal number of clusters. Bootstrap resampling further confirmed stability of the cluster assignments.

Propranolol blocks sympathetic optogenetic effects on heart rate.

Optogenetic stimulation of the most anterior sympathetic ganglion (APG; blue bar) produces developmentalstage–specific effects on heart rate. (Left) At 5 dpf, stimulation does not alter heart rate (green, control). Treatment with 20 µM propranolol (red) also shows no effect. (Right) At 9 dpf, stimulation induces a robust increase in heart rate under control conditions (green), which is abolished by propranolol treatment (red), consistent with β-adrenergic signaling mediating the sympathetic drive. Thin lines represent individual trials, thick lines the mean across fish, and shaded regions the s.e.m. Heart rate was normalized by subtracting and dividing by the resting baseline.

Dimensionality and statistical validation of kernel encoding.

(A) Eigenvalue spectrum of the kernel covariance matrix. The first two principal components (PCs) capture nearly all of the variance, with PC1 explaining ~85% and PC2 ~15%. (B) Cumulative variance explained as a function of the number of PCs. A 99% variance target (dashed line) is reached by the first two PCs. (C) Participation ratio (PR) analysis comparing bootstrap resampling (purple) with a shift-null distribution (gray). The observed effective dimensionality (red line, PR = 1.36) is far smaller than expected under the null, confirming low-dimensional structure. (D) Variance explained by the first three PCs. PC1 accounts for 84.5%, PC2 for 14.9%, and PC3 for only 0.3% of the variance, with PC2 ~44 times larger than PC3, indicating that kernels are well described by two dimensions. (E) Empirical cumulative distribution function (ECDF) of pooled cross-validated prediction performance (CV R2) across neurons. Real data (black) deviate significantly from the surrogate null (gray), Kolmogorov–Smirnov test, D = 0.234, p = 1.7 × 10-19. (F) Distribution of permutation p values pooled across neurons. The enrichment of very small p values indicates significant encoding across the population. (G) Fraction of neurons above threshold as a function of CV R2. Real data (black) exceed the null mean (light gray) and 95% confidence interval (dark gray shading) across thresholds, demonstrating more neurons with predictive power than expected by chance. (H) Sorted CV R2 values across all neurons. The heavy-tailed distribution of real neurons (black) lies above the null distribution (gray), showing that a subset of neurons carries disproportionately strong encoding signals.