Introduction

Large, heritable differences in complex behavior are observed even across closely-related species and across sexes within the same species. What are the molecular substrates, encoded by genes, that facilitate this variation in innate animal behavior? Two molecular traits known to govern behavior are: 1) abundances of specific neuronal cell types, which may be controlled by developmental genes involved in neurogenesis and apoptotic pathways, and 2) abundances of proteins within neuronal cells, which may be controlled by genes involved in transcriptional and translational regulation.

Variations in both types of traits have been documented to change behavior across species and sex. In Drosophila, loss of song-patterning neurons underlies a species-specific loss of a type of courtship song¹ while the expansion of an olfactory neuronal population contributes to Drosophila sechellia’s species-specific specialization for feeding on the noni fruit². In mammals, neuronal number differences mediating behavioral traits have largely been described across sexes rather than species. Neuronal population sizes, particularly of cell types within the preoptic area and bed nucleus of the stria terminalis, are thought to be regulated by sex-specific levels of neuroestrogens that control apoptosis, thereby leading to sex-biased cell abundances and sex-specific control of behavior^3,4.

Between cell populations, gene expression level differences of neuronally-important proteins such as neurotransmitter receptors, smell and taste receptors, and protein kinases are also known to cause behavioral differences across species (see [5] for review). For example, the expression level of the vasopressin 1a receptor gene in the ventral pallidum contributes to the evolution of partner preference formation in monogamous voles^6,7. Across sexes, differential expression of sex-steroid receptors in the mammalian brain is well-known⁸, and the resulting consequences for sex-specific gene expression and control of sex-specific behavior is an area of heavy investigation^9–12.

Advances in single-cell and single-nucleus RNA-sequencing have enabled high-throughput characterization of both cell type abundance evolution and cell type-specific expression evolution, and recent single-cell profiling of the brain of behaviorally-divergent species have revealed intriguing insights about both modes of evolution^13–15. However, the majority of comparative single-cell RNA-sequencing in mammals has been performed across fairly diverged species with many phenotypic differences, such as human and mouse (∼90 million years diverged¹⁶) or across primates^17–21 such as humans and marmosets (∼40 million years diverged¹⁶), lending limited insight into which of the specific behavioral difference may be affected by which of the many observed molecular changes.

Closely-related species of Peromyscus deer mice that exhibit variation in heritable behaviors gives us an opportunity to more precisely dissect the molecular neuronal mechanisms that can alter innate mammalian behaviors^22,23. In particular, P. maniculatus and P. polionotus diverged only ∼2 million years ago²⁴, but lie on opposite ends of the monogamy-promiscuity spectrum^25,26 and exhibit large, heritable differences in their parental care behavior^25,26. Both males and females of the monogamous P. polionotus are more parental than the promiscuous P. maniculatus^25,26. Further, this difference is much more distinct across males: while promiscuous P. maniculatus fathers exhibit little to no parental care in laboratory assays, monogamous P. polionotus fathers provide care to nearly the same level as P. polionotus mothers^25,26. These differences provide an opportunity to identify the molecular products that evolve to modulate social behavior, and how they do so in a sexually dimorphic manner.

Evolution of the hypothalamus has been implicated in social behavior differences across many species^27–29, including across these two species of Peromyscus²⁶. Here, we used single-nucleus RNA-sequencing (snRNA-seq) to profile the medial preoptic area (MPOA), a region of the hypothalamus known for its many cell types involved in mating and parenting behavior³⁰. We profiled males and females of both P. maniculatus and P. polionotus species and, by comparing to Mus data, identified MPOA cell types conserved across rodents. We identified cell type abundance and gene expression differences across species and sex, and hypothesized how these differences may contribute to species- and sex-specific behaviors. Finally, we characterized how sex-bias is evolving across these two Peromyscus species and propose potential cell type abundance and gene expression changes that may be co-evolving with monogamy.

Results

A transcriptional cell atlas of Peromyscus MPOA

We performed droplet-based snRNA-seq on the MPOA of 6 males and 6 females from each of the P. maniculatus and P. polionotus species (Methods) (Fig 1A). An atlas of the MPOA in Mus musculus was previously developed that combined single-cell RNA-seq data with imaging-based spatial information and data on specific neuronal populations activated by social behaviors³¹. In order to create a comparable dataset, we used a dissection strategy identical to that used for creating the Mus MPOA atlas to dissect a region containing the MPOA and surrounding nuclei (∼2.5 mm x 2.5 mm x .75 mm equivalent to Mus Bregma +0.5 to −0.6) (Methods, Fig 1B). To reduce batch effect, thereby increasing statistical power for detecting species and sex differences, we took advantage of the genetic diversity of our outbred study system and pooled nuclei from a male and female from each of the P. maniculatus and P. polionotus species (i.e. 4 samples) into each snRNA-seq run (Fig 1C). We collected a total of 105,647 nuclei across 6 replicates (i.e. 24 samples) (Methods, Supp Table 1). To computationally demultiplex the pooled data, we used bulk RNA-seq or whole genome sequencing data to call single nucleotide polymorphisms (SNPs) for each animal. We then used SNP data to assign sample information to each nucleus in the snRNA-seq dataset (Methods, Supp Figs 1-3, Supp Table 2). In total, we were able to confidently assign a single sample to 95,057 nuclei. The remaining 10,590 nuclei appeared to be multiplets containing two or more samples and were removed from subsequent analyses.

Figure 1.

Figure 2.

Figure 3.

Table 1.

We combined data across all samples and clustered the data using Seurat v4.1.3³² (Methods). We first performed a low-granularity clustering using the first 30 principal components (PCs) of our data and a low resolution parameter of 0.2 to obtain a small number of clusters. This procedure identified 21 major cell classes which were assigned as neurons, microglia, astrocytes, immature and mature oligodendrocytes, endothelial, and ependymal cells using known marker genes from Mus hypothalamus³³ (Supp Table 3, Supp Fig 4). We then focused only on neuronal nuclei (52,121 nuclei), and reclustered the data at a higher granularity using the first 100 PCs of our data (Methods) and a resolution parameter of 1.8. This resulted in 53 neuronal cell clusters (Fig 1D). We did not find major biases in sex or species across our clusters, as shown by even distribution of each sex and species across the UMAP (Supp Fig 5). Our cell clusters reflect known cell types of the MPOA. For example, arginine vasopressin (Avp) and oxytocin (Oxt) are exclusively and highly expressed in clusters 34 and 40 (Supp Fig 6), respectively, and reflect the well-known populations of AVP and OXT neurons that play important roles in social and reproductive behaviors. Additionally, vasoactive intestinal protein (Vip) is highly expressed in clusters 41 and 46 (Supp Fig 6), which together represent the VIP neurons of the suprachiasmatic nucleus that play a major role in regulating circadian rhythms. Further inspection of these two clusters finds specific expression of neuromedin S (Nms) only in cluster 41 and gastric releasing peptide (Grp) only in cluster 46 (Supp Fig 6), suggesting that these two clusters actually represent two subpopulations of VIP neurons, Vip/Nms and Vip/Grp, that have recently been reported as genetically and functionally distinct cell types^34,35.

Fig 4.

Fig 5.

Figure 6.

While some neuropeptidergic cell types were obviously identifiable in our data, they were the minority of our cell clusters. We sought to provide more biological context for the rest of our data by identifying homology between our cell clusters and cell types defined from the Mus MPOA atlas³¹. We used Seurat’s integration workflow³⁶ to identify conserved gene correlation structures across the Mus and Peromyscus datasets and transfer Mus cell type labels onto each Peromyscus nucleus (Methods, Fig 1E). Because the Mus MPOA atlas consists of separate inhibitory and excitatory datasets, we first defined each of our nuclei as inhibitory or excitatory based on whether it expressed higher levels of inhibitory markers Gad1 and Gad2 or excitatory marker Slc17a6 (Methods, Supp Fig 7). We then transferred inhibitory and excitatory Mus cell type labels only onto inhibitory and excitatory nuclei, respectively.

For each Peromyscus cluster, we calculated the proportion of nuclei within the cluster assigned to each Mus cell type and labeled clusters if at least 15% of the cluster was predicted to belong to the same Mus cell type label (Methods). This procedure identified mappings for 34 of our 53 neuronal cell clusters (Fig 1E, 1F, Supp Fig 8). (Mus cell type labels were denoted by [31] as inhibitory (i) or excitatory (e), and named after marker genes identified from Mus MPOA atlas.) The cell clusters for which we were able to assign Mus homology include four cell types activated by relevant social behaviors in Mus: cluster 5 (i18: Gal/Tac2), which is activated by mating behaviors; cluster 17 (i23: Crh/Nts) which is activated by parenting behaviors; and clusters 20 (i16: Gal/Th) and 9 (i20:Gal/Moxd1), which are activated by both mating and parenting behaviors³¹.

The Peromyscus cell clusters for which we were unable to assign Mus homology likely fall into two categories: First, dissection and/or anatomical differences between Peromyscus and Mus may influence the presence and absence of certain cell types. For example, cluster 2 forms a clear “island” on our UMAP and is transcriptionally specified by high expression of Pvalb, Gad1, and Gad2 (Fig 1D, Supp Figs 6-7). This cluster very likely consists of inhibitory interneurons that have been reported to reside primarily in the thalamus in Mus³⁷ but may have been inadvertently captured in our dataset due to differences in dissection or differences in spatial distribution of these interneurons across species. Second, some cell types may not be transcriptionally well-defined and homology would not be assignable based on gene expression data alone. For example, cell clusters that lie in the middle of our UMAP (e.g. clusters 0, 13, 25, 8, etc.) likely do not have strong transcriptional signatures captured by our snRNA-seq and were also less likely to be mapped to a homologous Mus cell type. Still, we were able to assign homology for the majority of our nuclei enabling us to better interpret cellular and molecular changes in the context of their biological functions.

Differential abundance of cell clusters across species and sex

With our Peromyscus MPOA cell atlas, we first asked if there were differences in cell abundance across species or sex (Fig 2A). To do so, we treated cell abundance as count data and used edgeR³⁸ to fit a generalized linear model that included replicate, sex, and species as covariates (Fig 2A, Methods). Importantly, edgeR controls for differential numbers of neurons sequenced across samples (Methods). We used edgeR to test for significant coefficients (i.e. differential abundance across species or sex) and accounted for multiple hypothesis testing with a false discovery rate (FDR) correction³⁹. Because we would later validate our results experimentally, we used a lenient FDR cutoff of 0.3 to initially identify candidate clusters showing differential abundance.

Across our 53 neuronal clusters (Supp Table 4), we found four clusters that were called differentially abundant across species and two that were called differentially abundant across sexes (Fig 2B, Supp Table 5). Arginine vasopressin (AVP) and oxytocin (OXT) neurons (clusters 34 and 40, respectively) were both ∼1.6x more abundant (AVP FDR = 0.19; OXT FDR = 0.09) in the promiscuous P. maniculatus than the P. polionotus (Fig 2C). The AVP neuron number difference corroborates a previous finding that P. maniculatus had significantly more AVP neurons than P. polionotus (∼1.7x) in the medial paraventricular nucleus⁴⁰. Furthermore, we previously reported based on bulk RNA-seq data that the expression level of Avp was nearly 3x higher in P. maniculatus hypothalamus compared to P. polionotus, and that this expression difference mediates the different amounts of parental nesting behavior performed by each species²⁶. Our data suggests this difference in Avp expression and parental nesting behavior may be driven, at least in part, by a neuronal number difference across species.

Two cell clusters were found to be more abundant in the monogamous P. polionotus compared to the promiscuous P. maniculatus: cluster 43 (i26:Tac1/Prok2) (3.4x higher; FDR = 0.08) and cluster 20 (i16:Gal/Th) (1.6x higher; FDR = 0.02) (Fig 2C). We were intrigued by the latter cluster because of its known roles in parental care. In Mus, i16:Gal/Th is activated by pup exposure and successful mating in both sexes³¹. This cell type likely represents a subset of the galanin-expressing neurons in the MPOA that were previously reported to elicit pup grooming in Mus fathers when optogenetically activated⁴¹. In Peromyscus, we also observed robust and specific expression of galanin in this cell cluster (Supp Figure 9). Previously, we have shown that P. polionotus of both sexes perform more parental care than their P. maniculatus counterparts, and that this behavioral difference is heritable²⁶. Thus, we found a positive association between innate parental care behavior and number of a galanin-expressing neuronal cell cluster across Peromyscus species.

Finally, we identified two cell clusters to be more abundant in males than females in both Peromyscus species: cluster 9 (i20:Gal/Moxd1) (1.7x higher; FDR = 0.007) and cluster 5 (i18:Gal/Tac2) (1.5x higher; FDR = 0.27) (Fig 2D). Both cell types are activated by mating behaviors in male Mus³¹. i20:Gal/Moxd1 neurons are found in both the sexually dimorphic nucleus of the preoptic area (SDN-POA) as well as the bed nucleus of the stria terminalis (BNST), and cannot be transcriptionally differentiated using snRNA-seq data³¹. Both the SDN-POA and BNST have been described to be sexually dimorphic and larger in males of many mammals including rats⁴², Mus⁴³, and humans^44,45. Here, we again found a male-bias of i20:Gal/Moxd1 in both a promiscuous and monogamous species.

Several of our differentially abundant cell types have known, highly-conserved gene markers: AVP neurons (marked by AVP), OXT neurons (marked by OXT), and i20: Gal/Moxd1 (marked by calbindin1 [CALB1])^4,46. We took advantage of this to experimentally validate our single-nucleus data, as well as to examine interspecies differences in the spatial distribution of differentially abundant cell types. We performed immunoreactive staining across the entire hypothalamus, counted cells within hypothalamic subregions, and tested for differential abundance in the direction identified from our single-nucleus data (Methods).

AVP and OXT neurons are formed from a common developmental lineage^47,48 and both primarily reside in the paraventricular nucleus (PVN) and the supraoptic nucleus (SON). (Our MPOA dissections for snRNA-seq would have captured the entire PVN and the anterior portion of the SON.) Our counts of AVP-immunoreactive (ir) neurons found 1.5x more cells in the PVN (p = 7.0e-4, one-sided Mann Whitney test) and 2.1x more cells in SON (p = 6.0e-6) (Fig 3A, 3B) in P. maniculatus compared to P. polionotus. Our counts of OXT-ir neurons found 3.1x more cells in the PVN (p = 2.8e-6) and 4.7x more cells in SON (p = 6.2e-5) (Fig 3C, 3D). Therefore both regional populations of AVP and OXT neurons appear to be contributing to their differential abundance across species.

In contrast, when we examined CALB1-ir cells (marking i20:Gal/Moxd1), we found regional differences in differential abundance. As expected, CALB1-ir cells were found in the SDN-POA and BNST. (Our MPOA dissections would have captured the entire SDN-POA and the ventral portion of the BNST.) Because the cell density of CALB1-ir cells is so high in the BNST, we calculated the area of immunoreactivity in the BNST rather than counting cells individually. The BNST was ∼2x larger in males in both P. maniculatus (p = 0.004, one-sided Mann Whitney test) and P. polionotus species (p = 0.03) (Fig 3E, 3F). However, in the SDN-POA, we found that the number of CALB1-ir cells was 1.8x more abundant in the promiscuous P. maniculatus males (p = 0.03, one-sided Mann Whitney test), but not significantly different across sexes in the monogamous P. polionotus (p = 0.2) (Fig 3E, 3F). The number of CALB1-ir cells in the SDN-POA was very low and suggested that the conserved differential abundance observed in the scRNA-seq data was likely driven primarily by large numbers of i20:Gal/Moxd1 cells of the BNST. However, our immunostaining suggests that the two regional populations of i20:Gal/Moxd1 are evolving differently across species. While sex bias in the BNST has remained conserved, sex bias in the SDN-POA has become more reduced in the monogamous P. polionotus.

In sum, we computationally identified six neuronal cell types of the MPOA that have changed in abundance across species or sex, and experimentally confirmed three cell types with abundance changes. Our findings align with previous work and identify a difference in AVP neurons, which had been reported to be differentially abundant across our focal species, and i20:Gal/Moxd1, which is known to be male-biased across mammalian species. The four neuronal cell types with differential abundance across species are candidates for underlying species differences in behavior. Differential abundance of AVP neurons may be partially responsible for a differential expression level of AVP across species that contributes to a difference in parental nesting behavior²⁶. Additionally, we observed the monogamous species has increased abundance of a galanin-expressing cell type (cluster 20, homologous to i16: Gal/Th) known to govern parental care behavior. Finally, our immunostaining found that regional populations of transcriptionally identical cell types may evolve differently across species, and suggests that regional-specificity of sex bias in i20:Gal/Moxd1 cells may be important to the evolution of sex-specific behaviors.

Differential neuronal gene expression across species and sex

We next asked what gene expression differences there were across species and sex. To do this, we pseudobulked our single-cell expression data and summed up gene counts across cells within the same cell cluster and sample (Fig 4A). We again used edgeR to fit a generalized linear model that included replicate, sex, and species as covariates (Methods). Because edgeR normalizes for library size across samples, cell type abundance is controlled for in our differential expression analysis.

Across species, we found a range of 18 to 3,401 differentially expressed (DE) genes in each cell cluster (FDR < 0.05), with no bias towards either species (Fig 4B, Supp Table 6). The number of DE genes in each cell cluster was highly correlated with the abundance of that cell cluster (Fig 4C). This is expected given that larger abundances is analogous to deeper sequencing of a cell population and results in more genes detected and more accurate estimates of gene counts. In total, we found 8,301 unique genes to be DE in one or more cell cluster. When compared to a background of expression-matched non-DE genes (Methods), DE genes were highly enriched in neuropeptide receptors (2.4x enrichment, FDR = 7e-04), and modestly enriched in neuropeptides (1.5x, FDR = 0.005), neurotransmitter receptors (1.3x, FDR = 0.02), and voltage-gated channels (1.3x, FDR = 0.03) (Fig 4D). We did not find a significant enrichment of genes involved in neurotransmitter synthesis or transport, and only a small, though marginally significant, enrichment of transcription factors (1.1x, FDR = 0.04) (Fig 4D).

Across sexes, we performed differential expression analysis separately for each species in order to identify changes in sex-biased expression across species (Supp Table 7, 8). Remarkably, we observed far fewer sex-biased genes in the monogamous species compared to the promiscuous species. We found a range of 0 to 40 DE genes (FDR < 0.05) in P. maniculatus cell clusters and 2 to 27 DE genes (FDR < 0.05) across P. polionotus cell clusters (Fig 5A). In total, we found 204 female- and 194 male-biased genes in P. maniculatus; in contrast, only 170 (83%) female- and 70 (36%) male-biased genes were detected in P. polionotus (Fig 5B). Our data suggest a substantial reduction of sex-biased genes in the monogamous species, driven primarily by a reduction of genes with male-biased expression.

Despite the short evolutionary timescales between our species, we found very few sex-biased genes shared across species (Fig 5B). Only eight genes were found to be female-biased in the same cell types in both species, including three known female-specific genes associated with X-inactivation (Xist, Tsix, and Jpx) (Supp Table 9). Eleven genes were found to be conserved in male-bias including Nrip1 (cluster 5, homologous to i18:Gal/Tac2), a nuclear protein that mediates estrogen signaling⁴⁹, and Ecel1 (cluster 9, homologous to i20:Gal/Moxd1), a neuronal protease⁵⁰ also known to be male-biased in the Mus musculus MPOA⁵¹ (Supp Table 9).

In both P. maniculatus and P. polionotus, two cell types had the most number of sex-biased genes: cluster 5 (i18:Gal/Tac2) and cluster 9 (i20:Gal/Moxd1) (Fig 5A,5B). When comparing the number of sex-biased genes detected to cell cluster abundance, clusters 5 and 9 were obvious outliers in the number of sex-biased genes detected (Fig 5C). The enrichment of sex-biased genes was apparent in both species, even though the two cell types displayed a large reduction in the number of male-biased genes in the monogamous P. polionotus.

Notably, these were the same two cell types we found to be differentially abundant across sexes (Fig 2D). However, we found the enrichment of sex-biased genes to be true even when we re-performed our analysis on data that was downsampled so that cell abundance was equal across sexes (Supp Fig 10A). Additionally, both clusters exhibited high expression of all four major sex-steroid hormones receptors – estrogen receptor α (Esr1), androgen receptor (Ar), prolactin receptor (Prlr), and progesterone receptor (Pgr) (Supp Fig 10B). These receptors may be playing roles regulating the observed sex differences in both cell abundance and gene expression⁸. Together, our data highlights two cell types that are enriched for sex-biased gene expression in both a promiscuous and monogamous species, suggesting that these cell types function in highly conserved, sex-specific pathways.

Finally, we performed a gene enrichment analysis on all sex-biased genes, compared to a background of expression-matched non-sex-biased genes (Methods). Despite the rapid turnover of sex-biased genes across species, we found sex-biased genes to be overwhelmingly enriched in neuropeptides in both species (29.5x enrichment, FDR < 2e-4 in P. maniculatus; 36.8x enrichment, FDR < 2e-4 in P. polionotus) (Fig 6A, 6B). Out of 69 annotated neuropeptides in our dataset, 21 were found to be sex-biased in at least one cell cluster in either species (Fig 6C). We found neuropeptides were more likely to be female-biased in both species (19 female-biased vs. 9 male-biased), although this did not reach statistical significance (Fisher’s exact test p = 0.06) (Fig 6C). We found examples of both species-specific and conserved sex-bias (Fig 6C) including 6 instances of neuropeptides sharing the same sex-bias in the same cell type across species (Table 1). Therefore, our data suggest that neuropeptides may be playing key roles in sex differences, but that their expression levels appear to be evolutionarily labile, potentially underlying species-specific variation in sex-specific behavior.

Our differential expression analysis identifies striking sex differences in gene expression in the MPOA. Overall, we found rapid turnover of sex-biased genes across closely-related species, including a significant reduction in the number of male-biased genes in the monogamous species. Conserved across species, however, we found an enrichment of sex-biased genes expressed in two cell types: cluster 9 (i20:Gal/Moxd1) and cluster 5 (i18:Gal/Tac2). These same cell types were also found to be differentially abundant across sexes in both species, further underscoring their potential significance in sex-specific biology. Finally, our data suggest that neuropeptides may be a central class of genes mediating sex differences in the brain.

Discussion

How the evolution of neuronal traits, such as cell type abundance and neuronal gene expression, contributes to variation in innate social behaviors is poorly understood. Here, we explored this question by performing single-nucleus RNA-sequencing in the MPOA of males and females of a monogamous, biparental species of Peromyscus and a species that is promiscuous and exhibits maternal-only care. We focused on the MPOA as it is a region known to be important for mating and parenting behaviors, and took advantage of an existing Mus MPOA single-cell atlas to assign homology between Peromyscus cell clusters and functionally characterized cell types.

Differential abundance of parenting-related cell types

Of the four cell types differentially abundant across species, three have obvious ties to parenting behavior: AVP (cluster 34), OXT (cluster 40), and i16:Gal/Th (cluster 20) neurons. The differential abundance we observed in AVP neurons both supports our prior finding that variation in Avp levels contributes to a parental nesting behavior difference across our focal species²⁶, and adds to our understanding of the evolutionary mechanisms involved. Bendesky et al.²⁶ previously used Quantitative Trait Locus (QTL) analysis to link parental nest building to two genomic loci, one of which encompasses the Avp locus. With bulk RNA-sequencing, Bendesky et al. found Avp was 3x more highly expressed in the promiscuous, less parental species (P. maniculatus), and further experimentally showed that this increase in Avp inhibits parental nesting behavior. The genetic linkage between the Avp locus and parental nesting was partially attributed to a cis-effect that leads to differential expression across species, though this cis-effect changes expression by ∼1.1x and does not explain the full difference in expression across species.

In our study, we found AVP neuron count was ∼2x higher in the promiscuous P. maniculatus, suggesting that a difference in neuron number may additionally contribute to Avp expression and behavioral differences. Additionally, we observed 1.8x higher expression of Avp within Avp neurons of P. maniculatus compared to P. polionotus, corroborating the cis-effect found by Bendesky, et al²⁶. Together, the combination of expression difference and abundance difference better matches the full 3x difference observed in bulk RNA-sequencing. Our data proposes that evolution may be modulating behavior by changing expression levels through both local, cis-regulatory pathways and upstream pathways in neurogenesis and neurodevelopment.

We also found OXT neuron number to be ∼3x higher in the promiscuous species. The direction of differential abundance is opposite to what may be expected given oxytocin’s known roles in pair bonding and monogamy⁵². The species difference in oxytocin neurons may simply be due to the shared developmental lineage between AVP and OXT neurons, which are formed from the same precursor cells and only differentiate from each other later in embryogenesis^47,48. Our data suggests that genetic changes upstream of AVP/OXT differentiation have evolved across our focal species. Still, more work needs to be done to understand if the OXT abundance difference leads to behavioral changes, or if, perhaps, they have been compensated for through other changes such as in expression and/or spatial distribution of the oxytocin receptor.

Finally, we find a subtype of galanin neurons, i16:Gal/Th (cluster 20), to be more abundant in monogamous species. Galanin neurons of the MPOA were first reported to be specifically activated by parental care behavior in Mus⁴¹, a finding that has since been replicated in cichlids⁵³ and frogs⁵⁴. These findings suggest galanin has an evolutionarily conserved role in parenting behavior. In a study of biparental poison frogs (Ranitomeya imitator), frogs that performed parental behavior in a laboratory assay had significantly more galanin neurons than non-parental frogs, independent of sex⁵⁴. This study and ours find a correlation between galanin neuronal number and increased propensity for parental care behavior, both within and across species. However, further investigation is required to understand whether these observations are causally linked.

Two cell types play important roles in sex-specific biology

Across sexes, we found two cell types, i20:Gal/Moxd1 (cluster 9) and i18:Gal/Tac2 (cluster 5), that were sex-biased in cell abundance and had an enrichment of sex-biased genes across both Peromyscus species. i20:Gal/Moxd1 (cluster 9) is known to be more abundant in males across many mammalian species in two hypothalamic regions: BNST and SDN-POA (see Campi et al. 2013⁵⁵ for an overview). However, our immunohistochemistry experiments found the male-bias in abundance to be significantly reduced in the SDN-POA of the monogamous species, though conserved across species in the BNST.

Our findings mirror a previous study comparing sister species of polygamous and monogamous voles, which found the volume of the SDN-POA to be sexually dimorphic in the polygamous species, but not in the monogamous species⁵⁶. Another study in P. californicus, an independent lineage of monogamous Peromyscus, found male-bias in the volume of SDN-POA, but noted that the volume difference was half of what had been reported in other, promiscuous rodents⁵⁵. In contrast, the male-bias found in the volume of the BNST was consistent with fold-differences reported in other rodents⁵⁵. Together, these findings suggest that 1) sex differences in the neuronal number of the BNST and SDN-POA can be independently evolving, even within transcriptionally similar cell types, and 2) sex differences in the neuronal number of the SDN-POA, but not BNST, may be co-evolving with monogamy.

Much less is known about the second sexually dimorphic cell type, i18:Gal/Tac2 (cluster 5). In Mus, it is activated by mating behaviors³¹ and also reported to be enriched in estrogen receptor (Esr1) expression and sex-biased genes⁵⁷. Thus, this cell type appears to play conserved sex-specific roles, possibly in mating and copulatory behaviors.

Sexual dimorphism in gene expression is reduced in monogamous species

Finally, we found significantly fewer sex-biased genes in the monogamous species compared to the promiscuous species. A prevailing hypothesis, first proposed by Darwin⁵⁸, suggests that males of promiscuous species have larger variance in reproductive success compared to those of monogamous species, leading to greater intrasexual competition and more extreme sexual dimorphisms. Correlations between monogamy and reduced sexual dimorphism have been repeatedly observed in morphological traits such as body size and animal ornamentation (e.g. coloration or plumage) across diverse taxa including primates, birds, fish, and insects^59–61. However, much less is known about whether the reduction of sexual dimorphism extends to molecular phenotypes such as gene expression. Work by us and others previously failed to find consistent genome-wide patterns in bulk RNA-sequencing data^62,63. The striking reduction we see in our snRNA-seq data highlights the importance of quantifying gene expression at the level of specific cell types to better understand evolution of gene expression and its relationship with organismal-level phenotypes.

Across our study species, there is also a reduction of sexual dimorphism in behavior, with the emergence of biparental care in the monogamous species (in contrast to maternal-only care exhibited by the promiscuous species). Whether the reduced sexual dimorphism in neuronal gene expression may be actually mediating the reduced sexual dimorphism in social behavior is an interesting avenue for future research.

Conclusion

This study begins to characterize molecular differences in a region of the brain responsible for mating and parenting behaviors across two species with divergent mating systems. Overall, we find both cell abundance and gene expression changes across sex and species, both of which contribute to differential levels of neuronally-important proteins. Though not explored here, changes in cell abundance may additionally accompany neuronal circuitry changes that affect behavior. Finally, we find intriguing differences in neuropeptidergic cell types as well as cell types expressing high levels of sex steroid hormone receptors, implicating two key classes of signaling molecules – neuropeptides and sex hormones – as important evolutionary dials for modulating species and sex-specific behavior. Future work, including experimental characterization of the effects of our observed molecular differences, is required to reveal the precise nature of how genes can encode for the substrates needed for generating variation in innate behavior.

Methods

Animal husbandry

P. maniculatus bairdii and P. polionotus subgriseus animals were originally acquired from the Peromyscus Genetic Stock Center (Columbia, SC, USA). We weaned animals at 23 days of age into single-sex, same-species groups. We maintained animals on a 16 h:8 h light:dark cycle at 22°C, housed them in standard mouse cages with corncob bedding and provided them with food and water ad libitum. Animal husbandry and experimental procedures were approved by the Harvard University Faculty of Arts and Sciences Institutional Animal Care and Use Committee (protocol 27-15-3).

MPOA dissection

Virgin animals between ages 60 and 103 days were euthanized by CO₂ and brains were rapidly dissected, flash-frozen in Tissue-Tek® Optimal Cutting Temperature Compound, and stored at - 70°C. The MPOA was dissected on a cryostat (Leica CM3050 S) at −20°C. Brains were sectioned coronally until reaching the widening of the anterior commissure. The dorsal and lateral regions of the brain were removed so that only a ∼2.5 mm x 2.5 mm region that contained the MPOA remained (Fig 1B). Five 150 μm thick sections were then made to capture the MPOA. Sections were placed in a 2 mL tube which was kept at −70°C until the day of nuclei extraction.

Nuclei extraction and 10x library construction

Nuclei were isolated as described in [64]. Briefly, dissected frozen tissue was placed in a single well of a 6 well plate, with 2 ml of extraction buffer. Mechanical dissociation was performed by trituration using a P1000 pipette, pipetting 1 ml of solution slowly up and down, 20 times. Trituration was repeated 5 times, with 2 minutes of rest between each time. The entire volume of the well was then passed twice through a 26 gauge needle into the same well. The tissue solution was transferred into a 50 ml Falcon tube filled with 30 ml wash buffer and spun in a swinging-bucket centrifuge for 10 min at 600g and 4°C. Following spinning, the majority of supernatant was discarded, leaving 500 μl of tissue solution. DAPI was added (1:1,000) and debris was removed from nuclei by FACS before nuclei were loaded into a 10X Genomics 3’ V3 Chip. A detailed protocol can be found at dx.doi.org/10.17504/protocols.io.bi62khge.

Single-nucleus RNA-sequencing

Our snRNA-seq data consists of two experiments: a pilot experiment that aimed to obtain a single replicate of 10,000 nuclei which was sequenced on a NovaSeq S1 lane (rep01), and a full experiment that aimed to obtain 5 replicates of 20,000 nuclei which were sequenced on an NovaSeq S4 flowcell (rep02 - rep06). Each replicate contained pooled nuclei from a male and female from each of the P. maniculatus and P. polionotus species (Fig 1C). Our final dataset captured 105,647 nuclei across 6 replicates sequenced to a depth of 1,623 - 2,320 median UMI counts / cell (Supp Table 1). The 10x library construction and sequencing were performed by the Harvard Bauer Core Facility.

Bulk RNA-sequencing and genotyping

For our pilot study (rep01), we dissected out a portion of the hypothalamus posterior to the MPOA from each animal for bulk RNA-sequencing. Samples were collected in Trizol (Invitrogen), and total RNA was extracted using the RNeasy mini kit (QIAGEN). Poly(A) enrichment, directional mRNA library preparation, and paired-end sequencing of 150 base pair reads on a NovaSeq 6000 were performed by Novogene. Samples were sequenced to a depth of ∼35 million read pairs. Reads were mapped to their respective genomes, Pman_2.1.3 (GCA_003704035.1) or Ppol_1.3.3 (GCA_003704135.2), using STAR with default parameters⁶⁵. Variants were then called using the joint genotyping procedure recommended by the Genome Analysis Toolkit (GATK)⁶⁶.

Whole genome sequencing (WGS) and genotyping

For our full experiment (rep02 - rep06), we performed WGS on animals for genotyping. We extracted genomic DNA from tail tissue using proteinase K digestion followed by the Maxwell RSC (Promega) DNA extraction. We prepared DNA libraries using the Nextera DNA Flex kit with Illumina adapters. All samples were pooled into a single library and sequenced on two NovaSeq S4 lanes with paired-end sequencing of 100 base pair reads. Reads were mapped to their respective genomes using bwa-mem with default parameters⁶⁷. The median read depth for these samples was 267,987,590 mapped reads per sample, which corresponds to ∼20X sequencing coverage across the genome. Variants were called using the GATK procedure described above.

Sample demultiplexing

To demultiplex pooled scRNA-seq runs, we first used 10x Genomics CellRanger v5.0.0 with default parameters to map reads to a combined P. maniculatus and P. polionotus transcriptome. Our in-house transcriptome annotations were generated using Comparative Annotational Toolkit⁶⁸ and GENCODE v15 annotations of Mus musculus (GRCm38 mm10) as reference. We used the standard CellRanger pipeline to filter out droplet barcodes associated with background noise⁶⁹. For each droplet that passed the background filter, we preliminarily assigned each droplet to the species for which the majority of reads mapped to (which we refer to as the “primary” species). Droplets clearly fell into two clusters: one cluster that represents singlets in which a small percentage of reads incorrectly maps to the secondary species and a second, sparser cluster of droplets that represent multiplets in which two species were truly captured (Supp Fig 1). To delineate between the two clusters, we examined the number of reads mapping to the primary genome and binned droplets in increments of 500 reads. Within each bin, we then plotted the distribution of reads mapping to the secondary genome. We fit this distribution with a mixture model of two gaussians using normalmixEM⁷⁰ and identified the point at which the two distributions cross (Supp Fig 1). We then found the line of best fit through all intersection points across the bins and used this line as the threshold between singlets and multiplets (Supp Fig 1). Within each species, we next used demuxlet⁷¹ to infer sample identity based on the known sample genotype data (Supp Fig 2). Demuxlet additionally identifies multiplets when mixtures of sample genotypes are detected.

For each 10x run, we expected our detectable multiplet rate to be ∼9%⁷² and our assigned detectable rate ranged from 8-13% across our replicates (Supp Table 2A). We found that our assigned singlets were evenly distributed across the four samples (Supp Fig 3A, Supp Table 2B) and the distribution of the number of genes detected per droplet was similar across samples and replicates (Supp Fig 3B). Finally, when we checked the expression level of the female-specific gene, Xist, we found that female-assigned cells had significantly higher Xist read counts than male-assigned cells and the median read counts of Xist across male-assigned cells was 0 (Supp Fig 3C). After sample demultiplexing, we used CellRanger v5.0.0 to map reads from each nucleus to the transcriptome of its assigned species and used these read counts for all subsequent analyses.

Cell type clustering

All normalization, data integration, and clustering of scRNA-seq data were performed with Seurat v4.1.3³². Within species, we normalized our scRNA-seq data using sctransform^73,74 and then integrated data across species using Seurat’s integration workflow with default parameters. Briefly, Seurat learns a conserved gene correlation structure between two data sets using canonical correlation analysis and identifies pairwise correspondences between single cells across datasets in order to transform data sets into a shared space³⁶. We created UMAPs using the runUMAP function, and clustered the integrated data by using the FindNeighbors function.

Homology mapping of inhibitory and excitatory cell types

We defined each cell cluster and all nuclei within that cluster as inhibitory or excitatory based on higher expression levels of Gad1 and Gad2 (inhibitory) or Slc17a6 (excitatory) (Supp Fig 7). We then mapped inhibitory and excitatory cell type labels from the Mus MPOA atlas to inhibitory and excitatory nuclei, respectively, using Seurat’s label transfer workflow with default parameters³⁶. Briefly, Seurat first identifies pairwise correspondences between single cells across datasets (“anchors”) and then “transfers” cell type labels by weighting the cell type identities of each anchor with their distance to the nuclei at hand. For each nuclei, this procedure calculates the probability that it belongs to each cell type label and assigns the cell type label which received the highest probability.

To find cluster-to-cell-type mappings, we calculated the proportion of each cluster that was assigned to each Mus cell type label (Supp Fig 8). We removed Mus cell type labels that lacked any notable marker genes and that indiscriminately mapped onto many of our clusters (e1:Glut, i1:Gaba, i7:Gaba, i9:Gaba, i11:Gaba, i13:Gaba, i12:Gaba, i15:Gaba, i39:Gaba) (Supp Fig 8). We then retained only cell clusters and Mus cell type labels if 15% of a cluster or more were predicted to belong to the same Mus cell type label. This resulted in 17 excitatory clusters mapping onto 11 Mus excitatory cell type labels and 17 inhibitory clusters mapping onto 18 Mus inhibitory cell type labels (Figure 1E,1F).

Differential abundance of cell clusters

Differential abundance analysis was conducted with edgeR³⁸. After clustering our neuronal snRNA-seq data, we created a count matrix of the number of nuclei within each cluster that belonged to each sample. Cluster sizes were normalized using TMM normalization, which relies on the assumption that there are minimal differences in cluster sizes across samples⁷⁵. We used the estimateDisp function with default parameters to estimate the dispersions across cell abundance counts and calcNormFactors function with default parameters to control for the differential numbers of total neurons profiled per sample. We then fit a generalized linear model with glmQLFit function and a design that included replicate, species, and sex as covariates. Finally, we used glmQLFTest to test for significant coefficients and corrected p-values with an FDR correction³⁹.

Immunohistochemistry and cell counting

We transcardially perfused mice with ice-cold 1x phosphate-buffered saline (PBS) with heparin (100 U/mL) and then with 4% paraformaldehyde in PBS. Brains were dissected out, postfixed for 24 hours at 4°C, cryopreserved in 30% sucrose, and stored at −70°C until subsequent use. To stain for protein, we sectioned brains at 30 μm. We aimed to capture the entire hypothalamus and collected sections beginning at the appearance of the corpus callosum until the appearance of the thalamus, keeping every second section (48 sections spanning ∼3 mm).

We blocked free floating sections for 1 hour in 10% normal goat serum (NGS) and 0.5% Triton-X in PBS. We incubated sections overnight with either rabbit anti-Avp (1:4000, ImmunoStar 20069), mouse anti-Oxt (1:2000, EMD Millipore MAB5296), or mouse anti-Calb1 antibody (1:1000, Millipore Sigma C9848) and 10% NGS in PBS. We used either donkey anti-rabbit Alexa 647 antibody (1:1000, Thermo Fisher Scientific A31573) or donkey anti-mouse Alexa 546 (1:1000, Thermo Fisher Scientific A11003) for secondary detection and mounted tissues with DAPI Fluoromount-G (SouthernBiotech, 0100-20). Slides were imaged on an AxioScan.Z1 slide scanner (Zeiss).

Following imaging, we exported images to .tif format and counted cells using a custom Fiji/ImageJ macro. Briefly, we manually outlined a region of interest (ROI) around the MPOA. Within the ROI, we binarized the image using the “Moments” method⁷⁶ for calculating the optimal threshold. We then used the watershed segmentation algorithm to separate touching particles. Finally, we used the Analyze Particles function to count the number of particles inside the ROI.

Differential expression analysis

To obtain pseudobulked data, gene counts were summed across all cells belonging to the same cell cluster and sample resulting in a matrix of 30,704 genes by 24 samples by 53 cell clusters. Differential expression analysis was then performed with edgeR³⁸, using the same procedure described above for differential abundance analysis. For differential expression analysis across species, we fit a generalized linear model with a design that included replicate, sex, and species as covariates. For differential expression analysis across sex, we performed edgeR analysis separately for each species and fit a generalized linear model that included only replicate and sex as covariates. We accounted for multiple hypothesis testing with a FDR correction³⁹ and used a FDR cutoff of 0.05 to define significant differential expression.

Gene enrichment analysis

Gene lists were manually curated by using the Gene Ontology Resource⁷⁷ and are available in Supp Table 10. To test for enrichments of each gene category, we created backgrounds with the following procedure: We binned all genes by their expression level in bins of 0.5 log₁₀(counts per million). For each gene in our foreground (i.e. a significantly differentially expressed gene), we randomly chose a gene from our entire dataset from the same expression bin. If the same gene was DE in multiple cell types, we matched the highest expression level. To obtain an empirical p-value, we repeated this procedure 5,000 times and calculated the number of times a background list contained equal or higher numbers of foreground genes within a gene category, and we accounted for multiple hypothesis testing of several gene categories using a FDR correction. We calculated the enrichment score by dividing the number of foreground genes by the mean number of background genes (across 5,000 random samples) within each gene category.

Data availability

All sequencing data has been deposited in GEO under accession number GSE272719.

Author contributions

JC, SRE, and HEH conceived of the study. JC performed all experiments with contribution from CK and PRR. JC analyzed data with contributions from SRE. JC wrote the manuscript with contributions from SRE and HEH. All authors approved of the final manuscript.

Acknowledgements

We thank Dhananjay Bambah-Mukuu and Harris Kaplan for advice on MPOA dissection; Charles Vanderburg and Naeem Nadaf for advice on single-nucleus extraction; and Doug Richardson for advice on cell counting. JC was supported by the Harvard Data Science Initiative and the National Institutes of Health (K99 GM146243-01). SRE and HEH were supported by the Howard Hughes Medical Institute.