Structural and genetic determinants of zebrafish functional brain networks

Abstract

Network science has revealed universal brain connectivity principles across species. However, several macroscopic network features established in human neuroimaging studies remain underexplored at cellular scales in small animal models. Here, we use whole-brain calcium imaging in larval zebrafish to investigate the structural and genetic basis of functional brain networks. Mesoscopic functional connectivity (FC) robustly captures the individuality of larvae and reflects structural connectivity (SC) derived from single-neuron reconstructions. Several connectome properties, including diffusion mechanisms and indirect pathways, predict interregional correlations. SC and FC share a hierarchical modular architecture, with structural modules shaping spontaneous and stimulus-driven activity patterns. Visual stimuli and tail monitoring reveal a functional gradient that coincides with sensorimotor functions. Last, regional expression levels of specific genes predict interregional FC. Our findings reproduce key mammalian brain network features, demonstrating larval zebrafish as a powerful model for studying large-scale network phenomena in a small and optically accessible vertebrate brain.

Fig. 1.. Brain-wide imaging of functional networks in zebrafish larvae.

(A) Experimental configuration for whole-brain two-photon imaging, light stimulus delivery, and behavioral monitoring in Tg(elavl3:H2B-GCaMP6s) zebrafish larvae. (B) Five planes highlighted from 21 functional imaging planes at different depths (left); side projection in one example fish (maximum intensity, right); white arrow indicates the dorsal side; refer to fig. S1 for image scales. (C) Centroids of automatically identified neurons from one larva, registered to mapZebrain and mapped into 70 brain regions (pseudocolors; retina/eye is not displayed). (D) Regional fluorescence time series from one representative larva, ordered from anterior to posterior regions; a pink rectangle highlights a period of dark-flash visual stimuli. (E) Group-averaged FC matrix ( $n=22$ larvae), ordered from anterior to posterior regions. (F) Network visualization of group-averaged FC, using the same color map as the previous panel; nodes represent brain region centroids, and bottom quartile FC edges are not displayed; network edges are mirrored across both hemispheres for visualization. (G) Functional network similarity across larvae, ordered from the most globally similar to the most globally dissimilar individual; five example individual networks are plotted on the right. (H) Functional network similarity of seven larvae imaged on consecutive days; similarity scores are averaged across both temporal directions, from dataset 1 to dataset 2 and vice versa; white asterisks denote maximal similarity values on each row. (I) Network similarity scores are significantly higher when comparing individuals to themselves [diagonal values from matrix in (H)] rather than different individuals [off-diagonal values from matrix in (H)] ( $P=3\times {10}^{-5}$ , t test). A, anterior; P, posterior; V, ventral; D, dorsal; L, left; R, right.

Fig. 2.. Structure-function coupling of zebrafish brain networks.

(A) 3D rendering of 1912 reconstructed neurons generated on mapZebrain.org. (B) Mesoscopic SC overlaid on brain anatomy; edge weights are color mapped with the next panel, and bottom quartile edges are not displayed. (C) Undirected SC matrix compared with FC; matrices are arbitrarily scaled for visualization. (D) Distributions of individually sampled FC values from edges that have either bidirectional, unidirectional, or no underlying structural connections (one-way ANOVA, $F=5048$ , $P=0$ ). (E) Functional and structural node degrees with arbitrary rescaled values for visualization; structural degrees are displayed on the left hemisphere, and functional degrees are displayed on the right. (F) Linear regression between structural and functional degrees (Pearson’s $r=0.693$ ). (G) Depiction of the FC modeling approach used in following panels; properties of the structural network, such as the shortest path between two nodes, are used to predict FC. (H) Four example matrices from an array of 40 graph properties used as predictors of FC in this study. (I) FC variance explained by each of the predictors, subtracted by the variance explained by SC; predictor families are indicated on the right, in the same order as the bar chart labels. n.s., not significant. (J) Relative (red) and cumulative (black) variance explained by each PC derived from the set of predictors. (K) Structural PC1 compared to FC; matrices are arbitrarily scaled for visualization. (L) Linear regression between upper triangle values of structural PC1 and FC (Pearson’s $r=0.684$ ). (M) Null distributions of structural PC1 correlations derived from null SC matrices (1000 matrices per distribution, empirical $P<0.001$).

Fig. 3.. Structural and functional modules are spatially overlapped at multiple hierarchical levels.

(A and D) Module coassignment probability matrices, hierarchically clustered using Ward’s linkage (for SC and FC, respectively); four modules are chosen arbitrarily as an intermediary hierarchical level for visualization. (B and E) Connectivity matrices, reordered by modules according to the previous panels (SC and FC, respectively); colors are mapped to squared edge weights for better contrast. (C and F) Brain region centroids colored by module membership (SC and FC, respectively); scale bars, 100 μm. (G) Structural and functional modularity at different hierarchical levels, compared against various null models. SC modularity curves are plotted above (purple), and FC curves are below (red); 99th percentiles of null distributions are plotted for null models. SCCM, spatially constrained configuration model; CM, configuration model; PR, phase-randomized time series (see Materials and Methods). (H) Module identity vectors plotted at different hierarchical levels, for both empirical data and null models; $N_m$ denotes the number of modules; representative labels are shown for each null model. (I) Adjusted Rand index computed on pairs of module labels from the previous panel at different hierarchical levels; 99th percentiles of null distributions are plotted. (J) Left: Visual examples of spatially compact (top) versus distributed (bottom) modules. Right: Boxplots of pairwise internal distances between regions belonging to the same modules, compiled across all hierarchical levels; all distributions are different with the exception of two comparisons, as indicated below the figure (one-way ANOVA, $P<0.001$ , followed by Tukey’s post hoc test for pairwise comparisons).

Fig. 4.. Mesoscopic coactivation patterns coincide with the modular organization of SC.

(A) Raster plot of z-scored regional activity across 22 larvae. Left and right hemispheres are plotted in the upper and lower halves, respectively. (B) Root-mean-square (RMS) of regional coactivation values for each of the corresponding frames in (A); the red line indicates a statistical threshold for high coactivation events ( $P<0.05$ , false discovery rate fixed at $\mathrm{\alpha }=0.05$ using the Benjamini-Hochberg procedure). (C) Correlation matrix of 243 detected events, separated into three clusters that are illustrated in the following panel. (D) t-SNE projection of 243 high-amplitude events, with three clusters identified using density-based clustering. (E and F) Raw data projections of the two main clusters, coregistered and averaged across larvae ( $n=22$ ); top; high-amplitude motor events; bottom: high-amplitude visual responses to dark stimuli; contrast is adjusted after averaging (see Materials and Methods). a.u., arbitrary units. (G and H) Average coactivation matrices of the two main clusters, reordered according to four structural modules; statistical significance of module coactivations is indicated on right matrix borders. (I and J) Null distributions of maximal module coactivation values in SA-preserving and SA-breaking surrogates; dashed lines indicate empirical coactivation values. (K) Two independent analyses recover modular forebrain activity that was weakly detected as cluster 3 in the previous clustering analysis. (L) Averaged raw fluorescence from frames corresponding to the 95th percentile of forebrain IC activity. (M) Average coactivation matrix of frames belonging to the 95th percentile of forebrain IC activity. (N) Null distributions of forebrain modularity, similar to (I) and (J).

Fig. 5.. Identification of sensorimotor neuronal populations.

(A) Example camera frame from high-speed monitoring of a head-restrained larva with tail tracking points. (B) Angular time series of 10 tail segments displaying two successive swimming events; colors correspond to tracking points on the previous panel. (C) Cumulative tail angle of one larva over a full experiment, with detected swimming events below. (D) Overlapped densities of motor-positive, motor-negative, dark-responsive, and light-responsive cells, projected on the mapZebrain template brain ( $n=18$ larvae, spatial)  $P<0.025$ . (E) Top; Motor and stimulus event vectors. Bottom: 50 example cells per functional category (color traces) with average population activity (black traces). (F) Individual cell densities from (D); the color map represents arbitrary spatial density units, relative to the maximum density observed in dark-responsive cells; the number of cells aggregated across animals and contributing to each distribution is indicated below each panel; spatial $P<0.025$ . (G) Number of cells per functional category per larva; error bars indicate SD; statistical comparisons are avoided due to slight methodological differences between motor and visual cell identification (see Materials and Methods). (H) Left: Regional measure of functional diversity using adjusted Shannon entropy; nodes reflect region centroids, and node sizes are redundant with the color map. Right: Schematized examples of low and high diversities in a brain region. (I) Number of overlapping cells for each nonorthogonal pair of functional categories per larva; error bars indicate SD. (J) Spatial density of polyfunctional cells detected across all animals, spatial $P<0.025$ . A, anterior; P, posterior; M+, motor positive; M−, motor negative; D, dark; L, light. Scale bars, 100 μm.

Fig. 6.. Regional sensorimotor functions coincide with the main functional gradient.

(A) Top: SC matrix and its first two diffusion gradients, with nodes denoting brain region centroids. Bottom: FC matrix and its first two diffusion gradients. (B) Eigenvalues of the first 50 structural diffusion gradients. (C) Absolute correlation $\mid r\mid$ between the first 50 SC and FC gradients; functional gradients are reordered to optimally match the structural gradients. (D) Eigenvalues of the first 50 functional gradients. (E) Top: Dark and motor-positive neuron centroids from Fig. 5, in red and blue, respectively. Bottom: Sensorimotor index of each brain region. (F) Top: FC matrix under spontaneous conditions. Middle: First FC gradient of the corresponding matrix above. Bottom: Correlation between sensorimotor index (SI) and first FC gradient. $r_s$ denotes the Spearman correlation coefficient, and P values are obtained through SA-preserving permutations (1000 permutations). (G) Similar to the previous panel, except with FC computed while including visual stimuli. (H) Similar to the previous panel, except with high-resolution FC computed from spontaneous activity. (I) Similar to the previous panel, except with high-resolution FC computed while including visual stimuli. Scale bars, 100 μm.

Fig. 7.. Gene coexpression predicts FC.

(A) Five example fluorescence in situ hybridization markers from mapZebrain (left: top view, 90th percentile projections), with corresponding regional annotations (right: in black outlines); intensities are scaled arbitrarily per gene for visualization. (B) Regional gene expression matrix; each row (brain region) is z-scored independently. (C) Left: CGE matrix across all 290 genes, compared with the FC matrix; $r_s$ denotes the Spearman coefficient. Right: CGE matrix of 12 optimized genes obtained through simulated annealing, compared with FC. (D) Average (full lines) and 95th percentile (dotted lines) of optimized NMI between CGE and FC for varying numbers of genes used in simulated annealing runs; pink curves correspond to empirical gene sets, whereas black curves correspond to spatially shuffled genes; a vertical line indicates the elbow of the average empirical curve; 1000 optimization runs per number of genes. (E) NMI distributions at the predetermined elbow point for empirical and null optimization results using 12 genes ( $P<0.001$ , t test). (F) Selection frequency of each gene across 1000 simulated annealing runs; pink bars correspond to empirical genes, whereas black bars correspond to spatially shuffled gene datasets; a horizontal bar denotes a selection threshold, with significant gene names indicated above. (G) 90th percentile intensity projections of 22 significant gene markers; gene names are indicated in the next panel, ordered by selection probability (left to right, top to bottom). (H) Overlap of 22 pseudocolored gene markers; pixelwise color contrast is used to highlight the most prominent genes at each location; colors are used to accentuate the global patterns, rather than to precisely distinguish individual genes (which are displayed separately in the previous panel). A, anterior; P, posterior; SA, spatial autocorrelation. Scale bars, 100 μm.