Cellular atlas of the anemonefish forebrainClusters accurately reflect major cell classes

We placed male A. ocellaris anemonefish together in pairs (n = 12 individuals, 6 pairs). After one member of each pair had completely transformed into a female and the pair had produced viable offspring numerous times (see Methods), we dissected the forebrain (containing the entire telencephalon and POA) for snRNA-seq. In total, we sequenced and analyzed 21,929 nuclei isolated from these 12 individuals (Fig. 1AB; see Supplementary Excel File 5 for fish body size and sequencing pool). On average, 2,026 genes and 4,322 transcripts were detected per nucleus (see Supplementary Fig. 1 for complete summary). Nuclei were first clustered into 19 “parent” and, nested within those parent clusters, 48 “child” clusters based on gene expression profiles (Fig. 1B). Clusters were then assigned to one of nine major cell classes (Fig. 1C) based on expression of canonical cell type-specific marker genes (Supplementary Fig. 2).

Fig. 1

Cellular atlas of the anemonefish forebrain contains all expected cell types and identifies key cell types involved in reproduction. (A) Experimental design and workflow. (B) Single nuclei plotted as points in UMAP coordinate space and colored according to cluster assignment. Clusters were assigned based on transcriptomic similarity/dissimilarity with other nuclei. Parent clusters were identified with numbers (1–19), and if a parent cluster contained multiple child clusters then each child cluster was assigned a letter (a, b, c, etc.). Colored fields surround parent clusters to aid in identifying their boundaries. (C) Clusters were assigned to cell types based on expression of canonical marker genes among the cluster’s marker-DEGs. For more information on marker genes used to assign cell types, see Supplementary Excel File 1. D-E) Dotplots depict the expression of select genes within each cluster, chosen a priori due to their documented involvement in the differential regulation of reproduction between sexes. Dot size indicates the percent of cells within the cluster that expressed the given gene. Dot color indicates the average normalized expression value within the cluster. These genes fall into two general categories: those involved in steroid signaling (D), and those that mark neuronal populations that release particular neuropeptides (E). Genes associated with small molecule neurotransmitters are shown in Supplementary Fig. 4

Neurons were broadly classified as glutamatergic (slc17a6a + and/or slc17a6b+) or GABAergic (slc32a1+, gad2 + and/or gad1b+), and glutamatergic and GABAergic neurons segregated strongly across clusters (Supplementary Fig. 2). Radial glia (cyp19a1b + and/or gfap+) [61, 62] were split into four child clusters, one of which (14d) strongly and selectively expressed crocc2 (an essential cilia protein component) and was thus identified as a putative ependymal or ciliated cell population [63]. Radial glial child clusters 14a-c expressed several genes involved in neurotransmitter reuptake (slc6a11b, slc1a2b, glula) that are expressed in mammalian astrocytes [64], supporting the hypothesis that radial glial cells in fish perform functions characteristic of mammalian astrocytes [61]. Other less abundant cell classes included oligodendrocytes (cluster 15a; mbpa+), oligodendrocyte precursor cells (cluster 16; cspg4+), microglia (cluster 17; p2ry12+), leukocytes (non-microglial) (cluster 18; ptprc+/p2ry12-), and endothelial cells (cluster 19; vegfd+) (Supplementary Fig. 2). All major cell types were detected in both sexes and in all individuals in roughly similar proportions (Supplementary Fig. 3). These clusters align strongly with major known cell types in the teleost forebrain.

Cluster transcriptomes identify distinct neuromodulatory signaling systems

To further characterize cluster biology, we investigated expression patterns of genes associated with specific neuromodulatory systems, i.e. major small molecule neurotransmitters, neuropeptides, and steroid hormone receptors (a subset is shown in Fig. 1D-E; complete results in Supplementary Figs. 4–6). The analysis revealed two cholinergic (chat + and/or slc18a3b+) neuronal clusters (7a, 13e) as well as clusters containing multiple nuclei that co-expressed the dopamine transporter and synthesis enzyme genes (th + and slc6a3+; most frequently in cluster 9f, less frequently in 13d, 13b, 13a, and 9d). Only a small number of nuclei co-expressed transporter and synthesis enzyme genes for norepinephrine (dbh + and slc6a2+; clusters 2a, 6a, and 9b) and serotonin (tph1 + and slc6a4b+; cluster 13d; Supplementary Fig. 4). We speculate that the low prevalence of nuclei co-expressing both transporter and synthesis enzymes reflects the high degree of sparsity of snRNA-seq data generally, and the low recovery of dbh, tph1, and slc6a4b transcripts in particular. Several clusters were enriched for neuropeptide genes implicated in reproduction in fish, including gnrh1, kiss1, galn, npy, oxt, avp, pdyn, penka, tac1, tac3a, ccka and cckb, scg2b, and adcyap1a and adcyap1b (Fig. 1E; Supplementary Fig. 5). Notably, cluster 13d selectively expressed a subset of neuropeptide-related genes (avp, oxt, galn) known to be selectively expressed in the fish POA, and multiple transcription factor markers known from mouse POA (hmx3a, hmx2, hmx3), supporting the conclusion that nuclei in this cluster are POA-derived [65]. Other POA neuropeptide genes (e.g. crhb, trh, crh) were expressed in cluster 13d as well as other clusters (e.g. clusters 1a, 1b, 11, 13a, 13b, 13c). Interestingly, gnrh1, a neuropeptide that is selectively expressed in the fish POA, was strongly and preferentially expressed in glutamatergic cluster 1d. One plausible explanation is that gnrh1 + nuclei in cluster 1d are anatomically derived from the POA but are developmentally derived from a distinct cell type that does not express other POA markers hmx3a, hmx2, or hmx3. Ultimately, these data are consistent with recent work in zebrafish showing that the POA is populated by multiple neuronal subpopulations from distinct neurodevelopmental cell lineages [66].

In addition to monoamine and neuropeptide systems, we also investigated steroid systems associated with sex and social dominance. Sex steroids (e.g. estrogens, androgens) are key mediators of sexual differentiation across vertebrate life [2], and glucocorticoid signaling is a key mediator of dominance status between the sexes [67, 68]. To identify cell populations that may respond to steroid hormone signals, we investigated cluster-specific expression patterns of sex steroid and glucocorticoid receptor genes (Fig. 1D; Supplementary Fig. 6). Putatively POA-derived cluster 13d was enriched for all sex steroid receptors (ar1, ar2, esr2a, esr2b, and pgr), consistent with known steroid hormone sensitivity of the POA in teleosts and other vertebrates species [2, 69, 70]. Notably, all radial glial child clusters (including putative ependymal cells) were enriched for esr2b and ar2, consistent with steroid modulation of neurogenesis [71]. Glucocorticoid receptor (nr3c1) was widely expressed and was enriched in parent clusters 1, 2, and 9 (specifically, glutamatergic 1a-d, f and g; glutamatergic 2a, b, d, e, h and i; and GABAergic 9a-f). Together these data support neuromodulatory signaling architecture as a major axis distinguishing clusters in this study and reveal a suite of candidate steroid hormone-sensitive neuronal and glial populations, including putatively POA-derived cluster 13d.

Gene expression supports strong anatomical delineation of clusters

Previous work in other species has shown that snRNA-seq clusters can segregate strongly along neuroanatomical dimensions [43, 72]. To investigate the potential neuroanatomical origins of nuclei in our dataset, we first examined expression patterns of previously published teleost neuroanatomical marker genes [43] across clusters. This revealed a strong pattern whereby glutamatergic neuronal clusters tended to express genetic markers of the dorsal pallium and GABAergic neuronal clusters tended to express genetic markers of the subpallium, mirroring findings in other teleosts [73]. Next, we employed SAMap to align anemonefish snRNA-seq clusters with previously published cichlid snRNA-seq clusters, which had been correlated to forebrain regions through integration with spatial transcriptomics. Briefly, SAMap estimates the cross-species transcriptional similarity of snRNA-seq cell types, accounting for the degree of protein sequence similarity between species [58]. We observed remarkably strong correspondence between anemonefish and cichlid clusters, with most clusters being mutual best matches, a pattern that spanned radial glia, oligodendrocytes, oligodendrocyte precursor cells, microglia, and many neuronal subtypes (Supplementary Fig. 7). These results support strong conservation of forebrain cell populations across teleost species. Many glutamatergic child clusters in anemonefish were transcriptionally similar to glutamatergic clusters in cichlids that spatially mapped to specific forebrain regions [72]. For example, anemonefish cluster 2 aligned strongly with cichlid cluster 10_Glut, which spatially mapped to the dorsal medial region of the telencephalon (Dm). In contrast, anemonefish cluster 3 aligned strongly with cichlid cluster 12_Glut, which spatially mapped to the granule subregion within the dorsal lateral pallium (Dl-g), a subregion that was recently shown to bear strong transcriptional similarities with the mammalian retrosplenial cortex [72]. These results support a level of neuroanatomical dimensionality to the anemonefish clusters, identifying specific brain regions from which specific clusters were likely derived. This provides a plausible anatomical basis to enrich our understanding of the sex differences revealed in subsequent analyses (see below).

Predicted cell-cell communication between clusters identifies predominant sending and receiving populations

In the brain, heterogeneous cell populations are interconnected to molecularly communicate in functional circuits. To investigate possible patterns of molecular communication among cell populations identified here, we used NeuronChat to perform cell-cell communication analysis [74] Briefly, this analysis uses gene expression associated with known ligand-receptor interactions to estimate the molecular potential for directional communication (“connection weight”) between all possible pairs of putative “sender” and “receiver” cell populations. We analyzed predicted communication potential among parent clusters, child clusters, and gene-defined subpopulations of the putatively POA-derived cluster 13d (Supplementary Excel File 6). The connections with highest weights were dominated by GABAergic and putatively interneuronal clusters 10 (pvalb6+) and 11 (sst+) as sending populations, followed by clusters 13a, 13b, 13f, 9c, and 13d galn + and 13d esr2a + neurons. Interestingly, on average, oligodendrocyte precursor cells had the strongest receiving weights, followed by glutamatergic neuronal cluster 3. Within putatively POA-derived cluster 13d, the strongest receiver subpopulations were ghrhrb+, kissr+, and galn + nuclei. Overall, predicted communication was strongly driven by expression patterns of Neurexin–Neuroligin gene pairs, followed by genes associated with GABAergic signaling (Supplementary Fig. 8). These results confirm that clusters in this study are heterogeneous in predicted communication with one another. Furthermore, they provide a plausible means of inferring not just sex differences across individual clusters but across full putative neural circuits in subsequent analyses (see below).

Sex change leads to extensive reorganization of the forebrainSex differences in cell type-specific gene expression are pervasive and pronounced

Previous work in anemonefish and other sex-changing fishes has shown sex differences in brain gene expression, but only for a small number of genes [35,36,37, 75]. However, these studies were limited by analysis of whole brain or whole brain region homogenate rather than analysis of individual cells or nuclei. Because the forebrain contains conserved brain regions and heterogeneous cell populations that regulate gonadal physiology, reproductive behavior, and social decision-making, we hypothesized that a subset of cell populations would display robust and widespread sex differences in gene expression. In support of this hypothesis, a large suite of neuronal and glial clusters exhibited strong sex differences in gene expression (Fig. 2AB; Supplementary Excel Files 7 and 8). Among parent clusters, a disproportionately large number of genes exhibited sex-associated expression (sex-DEGs) across a suite of neuronal populations (clusters 1, 2, 3, 9, 13) as well as microglia (cluster 17) (Supplementary Excel File 7). These patterns were explained in part by a set of constituent child clusters exhibiting a disproportionately large number of sex-DEGs (1a, 1b, 2a, 9a, 9c, 9d, 13a, 13b) as well as radial glial cluster 14b and mature oligodendrocyte cluster 15a (Fig. 2A, B; Supplementary Excel File 8). In gene ontology (GO) enrichment testing for biological processes, sex-DEGs were enriched for distinct functions in specific clusters (see Supplementary Results for more information). For example, sex-DEGs in cluster glutamatergic cluster 2a were uniquely enriched for “GABA-A receptor complex” and “chloride channel complex”, possibly reflecting transcriptional effects of sex differences in inhibitory input.

Several clusters showed strong directional biases in the proportions of sex-DEGs that were downregulated versus upregulated in females versus males (Fig. 2B; Supplementary Excel Files 7 and 8). More sex-DEGs were upregulated in females in glutamatergic neuronal parent clusters 3 and 13 (and a similar trend was observed in radial glial cluster 14); and more sex-DEGs were downregulated in females in glutamatergic neuronal parent clusters 1, 2, and 15. These patterns were mirrored in a set of constituent child clusters (2a, 13a, 13b, 13c), and were also observed in GABAergic cluster 9c, which had more sex-DEGs upregulated in females versus males. Enrichment testing revealed that sex-DEGs that were upregulated in females versus males (and vice versa) were related to specific biological functions (Supplementary Excel Files 7 and 8). For example, in parent clusters 1, 2, and 15, sex-DEGs that were upregulated in females were differentially enriched for GO categories related to apoptosis (clusters 2 and 15) or autophagy (cluster 1), whereas sex-DEGs that were downregulated in females were differentially enriched for multiple categories related to neuronal differentiation (Supplementary Results).

Fig. 2

Widespread sex differences in cell type-specific transcriptomes and proportions of cells in the anemonefish forebrain. (A) UMAP plot indicating the number of sex-DEGs by child cluster. Only clusters that were significantly enriched for sex-DEGs are shown. (B) Bar plot showing the number of sex-DEGs per cluster. Asterisks indicate clusters that were significantly enriched for sex-DEGs, also shown in panel A. Red filled portions of the bars show the number of sex-DEGs that were upregulated in females compared to males. Blue filled portions show the number of sex-DEGs that were upregulated in males compared to females. (C) UMAP plot showing sex differences in cluster proportions. Sex differences in proportion are presented as percent change from male to female. Only clusters that showed a significant proportion difference (q < 0.05) are colored. Eleven of 48 clusters (representing 27% of all cells) differed in proportion between sexes. D-E) The percentage of nuclei in parent cluster 2 and child cluster 2a in individual males and females are shown, females displayed higher percentages than males. F) In cluster 2a, average CytoTRACE scores (across all nuclei within individuals) were greater (indicating relatively less mature) in males than females. G) Within 2a, CytoTRACE scores indicated 2a.1 as the more differentiated/mature subcluster and 2a.2 as the less differentiated and more immature subcluster within 2a. H) The relatively mature subcluster 2a.1 was approximately twice as abundant in females than males. I) The relatively immature subcluster 2a.2 was similar in abundance between males and females. * above the bars indicates significance at q < 0.05

Sex differences in cellular composition of the forebrain are widespread and prominent

Previous work in teleosts (including anemonefish) and mammals identify sex differences in the relative proportions of specific cell populations in the forebrain, including within the POA [2, 28, 29, 39, 47, 76]. Based on this work, we tested for sex differences in the relative proportions of clusters and putatively POA-derived cell populations. A number of both neuronal and glial populations showed sex differences in relative proportions (Fig. 2C, complete results for parent clusters, child clusters, and POA-derived populations shown in Supplementary Excel Files 9–11, respectively). Among parent clusters, sex differences were observed in glutamatergic neuronal cluster 2 (more abundant in females) and radial glial cluster 14 (more abundant in males) (Supplementary Excel File 9). Among child clusters, the sexes differed in the relative proportions of multiple glutamatergic neuronal clusters (1f, 1 g, 2a, 2b, 2c, 2e, 13c and 15b), GABAergic cluster 13a, and radial glial clusters 14a and 14b (Fig. 2C; Supplementary Excel File 10). Taken together, these results support the idea that sex change leads to widespread differences in the neuronal and glial composition of the anemonefish forebrain, with 11 of 48 child clusters (comprising about 27% of all cells) being sexually differentiated in cell abundance.

Cluster 2 stands out as strongly differentiated in gene expression and cell abundance

Glutamatergic neuronal parent cluster 2 (Fig. 2D), and several embedded child clusters (2a, b, c and e) were 35–95% more abundant in females versus males (Supplementary Excel Files 9 and 10). Cluster 2 also had a disproportionately large number of sex-DEGs, a disproportionately large number of female-downregulated sex-DEGs, and sex-DEGs were enriched for biological processes related to neuronal differentiation and apoptosis (Fig. 2B, see “Cell type-specific gene expression” above). Among its constituent child clusters, sex differences in proportions were the most pronounced in 2a (Fig. 2E), which preferentially expressed cckb, the gene encoding cholecystokinin, and galr2a, the gene encoding galanin receptor 2. To explore cluster 2a in more detail, we investigated sex differences in neuronal maturity using CytoTRACE, a bioinformatics tool that estimates cellular maturity based on the gene expression. Cluster 2a nuclei were predicted to be relatively more mature in females compared to males (Fig. 2F; cluster 2a was the only child cluster with a sex difference in CytoTRACE score; Supplementary Excel File 12; 0 indicates high relative maturity, 1 indicates low relative maturity). To better understand this pattern, we re-clustered 2a independently from other nuclei and identified two subpopulations of 2a that differed in relative maturity: cluster 2a.1 was predicted to be relatively more mature and 2a.2 was predicted to be relatively less mature (Fig. 2G). Females displayed approximately double the proportion of the relatively mature cells in subcluster 2a.1 compared to males (Fig. 2H), but similar proportions of the relatively immature cells in subcluster 2a.2 (Fig. 2I). The relatively mature subpopulation expressed both cckb and galr2a. These data support a model in which protandrous sex change involves maturation and increased abundance of a population of glutamatergic, galanin-sensitive, CCK-secreting neurons in the Dm, while maintaining consistent proportions of the relatively immature subpopulation. CCK signaling is known to play a role in the neural control of reproduction in fish and a recent study in zebrafish shows that CCK signaling is required for the development of female gonads [69, 77].

Sex change involves substantial reorganization of the steroid-sensitive and neuropeptidergic subpopulations of the POA

Given their established roles in regulating sexually-differentiated physiology and behavior in anemonefish [78,79,80] we hypothesized that steroid-sensitive and neuropeptide-synthesizing neuronal populations would also exhibit sexually dimorphic proportions and/or transcriptome. To investigate this, we tested for differences in the proportions of nuclei expressing specific steroid hormone receptors, neuropeptides or neuropeptide receptors known to be involved in the neural control of reproduction, specifically within putatively POA-derived cluster 13d (Supplementary Excel File 11). In support of our hypothesis, we found males displayed a significantly greater proportion of nuclei expressing adcyap1a, esr2a (Fig. 3A), npy7r, and tacr3, whereas females displayed a greater proportion of nuclei expressing esr2b (Fig. 3B), cckbr2 (Fig. 3C), and pgr (Fig. 3D). Notably, esr2b and pgr were strongly co-expressed across nuclei (see Supplementary Fig. 9), suggesting that increased proportions of pgr + and esr2b + in females may be driven by the same underlying cell population either increasing in abundance or, alternatively, upregulating esr2b and pgr. Taken together, these results support the idea that transformation from male to female is associated with substantial molecular and cellular reorganization within the preoptic area, including differentiation of signaling systems known to regulate reproductive physiology in other teleosts and plausibly playing a similar role in anemonefish (e.g. estrogen, cholecystokin). Such changes may support the transition from male to female and/or maintenance of sex-specific physiology and behavior.

Fig. 3

The putative POA-derived cluster contains sexually-differentiated subpopulations marked by receptors for estrogen, progesterone, and CCK. Boxplots show sex differences in proportions of sex steroid receptor expressing neurons in cluster 13d that map to the POA. A) Males display a greater proportion of cells in 13d that are esr2a + than females. UMAP plot of cells in cluster 13d, cells colored blue if they express esr2a. B-D) Females display a greater proportion of cells in 13d that are esr2b+, cckbr2+, and pgr + as compared to males. * above the bars indicates significance p < 0.05

Sex differences in radial glia are correlated with specific neuronal sex differences

In adult teleosts, radial glial cells supply new neurons throughout the brain and are critically involved in the differentiation and maturation of progenitor neurons [12, 43, 61, 62, 71, 81]. In our earlier analyses, radial glial clusters 14a and 14b were more abundant in males, and cluster 14b contained a disproportionately large number of sex-DEGs (see above). Earlier analyses also identified sex differences in neuronal proportions and expression of genes related to neurogenesis. Given these findings together, we hypothesized that radial glia play a role in sex change and/or maintenance of sex-specific phenotypes. To test this hypothesis, we re-clustered radial glia independently of all other nuclei and identified nine subclusters. These subclusters differed in transcriptional signatures of cell maturity, radial glial functional states, and expression of notable marker genes including aromatase (strongly expressed in subclusters 2, 4, 5 and 6; Supplementary Fig. 10, Supplementary Excel Files 12 and 13). Further analysis suggested sex differences in multiple dimensions of radial glial biology. Radial glial subclusters 0 and 3 displayed a disproportionately large number of sex-DEGs (Supplementary Excel File 14; Supplementary Fig. 11). The composition of specific subclusters also differed between the sexes, with females exhibiting greater proportions of the relatively immature subcluster 8 and males exhibiting greater proportions of the relatively mature subcluster 5 (Supplementary Excel File 15). Females showed stronger expression of aromatase compared to males in subcluster 4 (Supplementary Excel File 14), mirroring whole-brain patterns in A. ocellaris and other species [35, 75, 82, 83].

We next hypothesized that radial glia may play a role in producing some of the neuronal sex differences identified above. To test one aspect of this, we asked whether radial glial proportions or aromatase expression were correlated with proportions of particular neuronal populations. Indeed, the relative proportions of radial glial subclusters 0 and 5 were strongly and inversely correlated with the proportions of a subset of neuronal populations that differed in abundance between the sexes (relatively mature 2a, 2c, 2e, 13c, 15b, 1 g), and multiple putatively POA-derived subpopulations (13d tacr3+, 13d adcyap1+, 13d pgr