Single-nucleus and spatial transcriptomics identify brain landscape of gene regulatory networks associated with behavioral maturation in honeybees

Abstract

Animal behavior is linked to the gene regulatory network (GRN) coordinating gene expression in the brain. Eusocial honeybees, with their natural behavioral plasticity, provide an excellent model for exploring the connection between brain activity and behavior. Using single-nucleus RNA sequencing and spatial transcriptomics, we analyze the expression patterns of brain cells associated with the behavioral maturation from nursing to foraging. Integrating spatial and cellular data uncovered cell-type and spatial heterogeneity in GRN organization. Interestingly, the stripe regulon is explicitly activated in foragers' small Keyon cells, which are implicated in spatial learning and navigation. When worker age is controlled in artificial colonies, stripe and its key targets remained highly expressed in the KC regions of bees performing foraging tasks. These findings suggest that specific GRNs coordinate individual brain cell activity during behavioral transitions, shedding light on GRN-driven brain heterogeneity and its role in the division of labor of social life.

Fig. 1. Cell types of the honeybee worker brains identified by single-nucleus transcriptomes.

a Schematic diagram of the experimental design. b Annotated tSNE visualization of the clustering of 121,247 single-nucleus transcriptomes obtained from nurse (70,060) and forager (51,187) bees. KC: Kenyon cells; Glia: glial cells; OLC: optic lobe cells; OPN: olfactory projection neurons; Hem: hemocytes; IPC: insulin-producing cells; Others: unknown cell type. Cluster identifiers from Seurat clustering at a resolution of 0.2 are indicated in parentheses. c Selected marker genes for the cell clusters annotated in nurse and forager bee brains. The expression levels of indicated genes were evaluated by the collapsed pseudo-bulk expression in each cluster as a percentage of unique molecular identifiers (UMIs) of the total cluster. Bars represent the means of 54 biological replicates + SEM. Cluster labels correspond to those shown in Supplementary Fig. S4. Source data are provided as a Source Data file. d Normalized UMIs per cell for the markers of Kenyon (dopR2, mAChR-DM1) and glial (repo, glnS) cells Heatmap plotted over global tSNE.

Fig. 2. Spatial transcriptomic atlas of the brain sections of nurse and forager bees revealed subclusters in anatomical regions.

a A schematic view of the honeybee brain with the main structures, including the mushroom body (MB), ocelli (OC), vertical lobe (VL), antennal lobe (AL), lobula (LO), medulla (ME), lamina (LA), and retina (RE). b Hematoxylin-eosin staining of one of the representative brain sections of the honeybee. c Nucleic-acid dye staining of the brain section shown in (b) attached on the Stereo-seq chip. Unsupervised clustering of the brain sections from nurse (d) and forager (e) bees based on the Stereo-seq data at bin-50 (50 × 50 DNB) resolution. Bins are colored based on the annotation of the clustering. Projection of our snRNA-seq data onto the brain sections of nurse (f) and forager (g) bees using Tangram. The color represents the relative probability of each cell type being mapped to the specific locations within the sections. Scale bars, 0.3 mm.

Fig. 3. Subclustering of Kenyon cells from honeybee worker brains.

a tSNE visualization for the subclustering of Kenyon cells. lKC: Class I large-type Kenyon cells; mKC: Class I middle-type Kenyon cells; sKC: Class I small-type Kenyon cells; II-KC: Class II Kenyon cells. Cluster identifiers from Seurat clustering at a resolution of 0.2 are indicated in parentheses. b Discernible marker genes for the subclusters annotated in Kenyon cells. Bars represent the means of 54 biological replicates +SEM. c Heatmap showing the average expression scale for marker genes in subclusters of Kenyon cells. d Representative GO terms enriched by the DEGs in different subclusters of KCs. Fold change is calculated as the proportion of detected genes within each GO term. Dot size represents the normalized fold change. Significance is assessed using the two-sided hypergeometric test. e Correspondence of KC clusters between Apis mellifera, Monomorium pharaonis, and Drosophila melanogaster determined by MetaNeighbor. The lines with different thickness link KC clusters from A. mellifera with those from M. pharaonis and D. melanogaster according to the AUROC scores. Only transcription similarities with AUROC > 0.6 are shown (Supplementary Fig. S13). Source data are provided as a Source Data file.

Fig. 4. Subclustering of optic lobe cells from honeybee worker brains.

a tSNE visualization for the subclustering of OLCs from the worker brains, which are grouped into eight subtypes. Cluster identifiers from Seurat clustering at a resolution of 0.2 are indicated in parentheses. b Correspondence of OLC clusters between Apis mellifera, Monomorium pharaonis, and Drosophila melanogaster. The lines with different thickness link OLC clusters from A. mellifera with those from M. pharaonic and D. melanogaster according to the AUROC scores. Source data are provided as a Source Data file. c A schematic diagram of the structural characteristics of neurons from the optic ganglia of the honeybees (adapted from Ribi and Scheel). svf, short visual fibers; lvf, long visual fibers; L, L-fibre; Tnf, tangential fiber with narrow-field ending; Twf, tangential fiber with wide-field ending; Tm, transmedullary neuron. Spatial visualization of the distribution of the OLC4 subcluster in the optic lobe regions of the brain sections of nurse (d) and forager (e) bees inferred by Tangram. The color of each merged bin is scaled according to the Tangram scores based on the projection of the snRNA-seq dataset. Scale bars, 0.3 mm.

Fig. 5. Worker castes in honeybees induce changes in brain cell populations and gene expression profiles.

Bar plots display DEGs identified in different cell types of (a) foragers (n = 1265) and (b) nurses (n = 1016). Significance assessed by two-sided quasi-likelihood F-tests in edgeR, adjusted by Benjamini–Hochberg. DEGs were identified as genes with an absolute log₂FC > 1 and an adjusted P value < 0.05. Bars indicate whether a gene is a DEG in a given cell type. Overlap of DEGs between single-cell clusters and bulk brain^,, in forager (c) and nurse (d) bees. The color of each square represents the proportion of overlapping DEGs, normalized within each bulk dataset. Significance is assessed using the two-sided Fisher’s exact test. e Heatmap showing upregulated DEGs from different subclusters of cells in forager (upper) or nurse (down). Colors represent the normalized log₂FC between nurses and foragers across different cell types. Source data are provided as a Source Data file.

Fig. 6. Gene regulatory networks underlie different cell types.

a Normalized enrichment score (NES) and the size of all 184 regulons identified from honeybee brains. Nineteen regulons for transcription factors (TFs) previously associated with behavioral maturation of honeybees^,, are indicated in red, with their highest-scoring motif shown. The inset shows the size distribution for target genes of each regulon. b SCENIC-based tSNE showing 14 clusters based on regulon activity (AUCell score) for each cell. Proportions of different (c) cell types and (d) subtypes in each SCENIC cluster. Cluster labels correspond to those shown in (b). SCENIC-based tSNEs showing the cluster of cell types (e) and subtypes (f) with similar regulatory states. Source data are provided as a Source Data file.

Fig. 7. Spatial heterogeneity of the gene regulatory networks in the honeybee brain associated with behavioral maturation.

a Differential regulon activity between forager and nurse bee brains across cell subtypes. Nineteen regulons previously associated with bee maturation are shown here^,,. The size of the square represents a -log₁₀(adj. P) of the two-sided t-test, and P value is adjusted using the Benjamini–Hochberg method. The color indicates the normalized log₂FC across different cell types. Log₂FC > 0 indicates upregulation in foragers. b Regulatory network for the co-regulatory factors CrebB, stripe, Hr38. Regulons for each TF are represented by different edge colors. The target genes are shown in circles; the TFs and the corresponding motifs are shown in hexagon nodes. c Volcano plot showing the regulons with differential activity scores in functional subregions of the mushroom body and the optic lobe. The significance is determined by two-sided t-test, adjusted by Benjamini–Hochberg correction. d Spatial visualization of the normalized activity score of the regulon stripe at mushroom body regions of nurse and forager brains. Scale bars, 0.3 mm. Fluorescent in situ hybridization of stripe (e), Hr38 (f) and caveolin-3 (g) in the bee brains. Simultaneous detection of target gene mRNA (red) and (4’, 6-diamidino-2-phenylindole) DAPI-stained nuclei (white), using confocal microscopy. Brain sections from 22-day forager, 10-day forager and 10-day nurse bees, respectively. Scale bars = 20 um. Each experiment was repeated independently three times. h The bar plot displays the variation in fluorescence intensity within the mushroom body for stripe (top) and its target genes Hr38 (middle) and caveolin-3 (bottom). The comparative analysis was conducted using the two-sided LSD test, with a significance level of 0.05, adjusted by Benjamini–Hochberg correction. Each point represents a randomly selected fluorescence region (n = 5). Bars represent the average intensity + SEM. Source data are provided as a Source Data file.