Introduction

Social attachments, such as the bonds between parents and children, family and kin, friends, and the long-term relationships between mates and romantic partners form the basis of interpersonal and group dynamics within societies^1–7. Various selective pressures have given rise to diverse patterns of complex social organization involving long-term relationships between members of a species, across mammals and other taxa, one of the most robust and intriguing of which is the enduring bond between mates^2,4,8–13. In mammals, ∼4% of species demonstrate social monogamy, forming persistent pair bonds between mating partners^9,14–16. Comparative studies between closely related species that display distinct patterns of these social behaviors, including social monogamy, reveal correlated species-specific differences in the nonapeptide hormones oxytocin (Oxt) and arginine vasopressin (Avp), or their homologs, suggesting that they serve as critical modulators of affiliative behaviors and pair bonding^{9–11,16–20}. Despite the central role of social attachment in human behavior and a rich history of investigation into the cognitive processes underlying social cognition and affiliative behaviors, it has been difficult to probe the neural mechanisms underlying these fascinating behaviors at the cellular, molecular, and genetic level, as most species that serve as model organisms for such studies in the laboratory do not display long-term social attachment as adults ^6,21–27.

Prairie voles (Microtus ochrogaster) display long-term social attachment such that mating partners show an enduring pair bond and social monogamy^13,28–31. Like many complex innate behaviors, pair bonding consists of distinct modules: the development of increased prosocial behaviors with partners over strangers (partner preference), and increased agonistic behaviors towards novel opposite sex potential mates (stranger rejection)^{13,30,32–37}. As in humans and other species that show affiliative social behaviors, disruption of pair bonds results in increased anxiety-like behaviors and markers of stress, supporting integrated neural and physiologic mechanisms that facilitate the preservation of such attachments^38–40. Pioneering work in prairie voles identified Oxt and Avp as critical mediators of pair bonding^41,42. Interspecific variations in the patterns of expression of the oxytocin receptor (Oxtr) and the vasopressin 1a receptor (V1ar) correlate with the potential for pair bonding within closely related vole species^43–49. Consistently, pharmacologic inhibition of Oxt and Avp signaling via their respective cognate receptors disrupts pair bonding, while exogenous administration of these hormones promotes pair bonding without mating^28,50,51.

Given the central role of Oxt signaling in affiliative behaviors, we tested the genetic requirement for signaling via Oxtr for pair bonding and social behaviors in prairie voles. We generated prairie voles lacking Oxtr and unexpectedly observed that these animals continue to display a strong preference for partners following mating and cohabitation⁵². To gain deeper insight into the modulation of attachment and social behaviors by Oxtr function, we developed a series of behavioral assays to probe distinct aspects of pair bonding and associated social behaviors. To interrogate the molecular mechanisms mediating social attachment, we profiled patterns of gene expression in the nucleus accumbens (NAc), a key site of species differences associated with social monogamy that is dramatically enriched for Oxt binding in prairie voles^45,53,54. We find that Oxtr controls the timing of pair bonding in prairie voles and promiscuity independent of partner preference formation. Oxtr sex-specifically controls reciprocal social behaviors, such that naïve wildtype (WT) males preferentially interact with wildtype females over those lacking Oxtr, while wildtype females show no such bias towards WT males. Oxtr function regulates patterns of gene expression in the NAc, as well as other regions of the vole brain implicated in social and attachment behaviors, in distinct ways in each sex. Finally, we find that Oxtr function regulates the development of Oxt and Avp expressing neurons in the paraventricular nucleus, a principal source of these neuropeptides. Our studies reveal genetically separable components of pair bonding and suggest that Oxtr signaling sex-specifically participates in the development of the circuits underlying attachment behaviors.

Results

Oxtr function facilitates the formation of partner preference

We previously found that Oxtr is not required for the display of partner preference following 1 week of cohabitation⁵². However, acute modulation of Oxtr signaling can alter patterns of social behaviors between individuals^29,50,55. These observations suggest that changes in Oxtr function may influence early social interactions that contribute to the trajectory of pair bond formation and/or specific aspects of attachment behaviors following the formation of partner preference. We first wished to determine whether Oxtr function facilitated pair bond formation by examining the display of partner preference following shorter periods of cohabitation. Consistent with previous studies showing that the cohabitation time required for the display of partner preference differs between males and females, we find that WT females demonstrate partner preference after only 6 hours of cohousing with a WT male partner, while WT males require 5 days following introduction to a potential WT mate (Fig. S1A-C)⁵⁶. We therefore cohoused WT or Oxtr^1-/- females with WT males for 6 hours, and WT or Oxtr^1-/- males with WT females for 5 days following introduction (short cohabitation), or in these combinations for 1 or 7 days (long cohabitation) for females or males, respectively (Fig. 1A, Fig. S1A).

Figure 1:

Both males and females lacking Oxtr show an absence of partner preference following short cohabitation with WT partners (Fig. 1B). Oxtr^1-/- females, but not males, showed a decrease in the fraction of huddling time with their partner and an increase in the fraction of huddling time with a stranger compared to WT controls (Fig. 1C-D). Oxtr^1-/- males present an absence of partner preference (Fig. 1B), displayed due to both decreased time with their cohoused partner and increased time with strangers (Fig. 1C-D). Consistent with our previous studies, we find no difference between WT and Oxtr^1-/- females when the pairs are cohoused for the longer duration (Fig. S1B-F).

To further understand the dynamics of social interactions that contribute to the difference in huddling behavior between WT animals vs their siblings lacking Oxtr following short cohabitation, we quantified the cumulative time subjects spent huddling with their cohoused partners or strangers over the course of the assay. Oxtr^1-/- females show a significant decrease in the amount of time spent huddling with partners in the last 30 mins of the assay (150-180mins) (Fig. 1E). We also observe a significant increase in huddling time of Oxtr^1-/- males with strangers during the middle periods (75 – 130mins) (Fig. 1F), suggesting that they engage in more early prosocial behavior with strangers than WT males. Consistent with our previous work, we find no difference between cumulative huddling times with either partner or stranger between WTs and Oxtr^1-/- animals of both sexes after long cohabitation periods (Fig. 1G-H)⁴⁵. These observations demonstrate that Oxtr function changes the temporal dynamics of social interactions with both familiar and unfamiliar animal in a sex-specific manner. Taken together, our findings suggest that Oxtr function acts early in the process of pair bonding to facilitate the formation of the partner preference.

Oxtr suppresses promiscuous behaviors in a state dependent manner in females

Oxtr influences the dynamics of interactions with partners and strangers differently between the sexes. We therefore wished to examine the detailed patterns of these interactions separately in females and males (Fig. 2A). We examined the distribution of huddles by WT vs Oxtr^1-/- females following short cohabitation as a measure of prosocial behavior displayed by these animals. We find that the top quartile of huddle durations (long huddles) from each subject animal make up an average of ∼90% of the total huddle duration, functionally separating these interactions from incidental side-by-side contact that may not constitute affiliative huddling (Fig. S2A Left). We classified long huddles as partner huddles or stranger huddles and discovered that, consistent with the overall time spent, the median huddle duration carried out by Oxtr^1-/- females is significantly shorter with partners and significantly longer with strangers, compared to WT females (Fig. 2B). The longest huddle displayed with partner is significantly longer for WT females compared to Oxtr^1-/- females though we find no difference in the length of the overall longest huddle. This is due to some Oxtr^1-/- females conducting their longest huddle with the stranger instead (Fig. S2B). Oxtr^1-/- females additionally executed a larger number of huddles only with the stranger, but not the partner, compared to WT controls (Fig. 2C). Taken together, these results suggest that Oxtr function in females both facilitates prosocial behavior with partners and suppresses such behavior with strangers during the early stages of pair bond formation.

Figure 2:

We next wished to determine whether the loss of Oxtr function influences patterns of prosocial behaviors after the display of partner preference. We therefore examined patterns of huddling, in particular the duration and frequency of these interactions, following long cohabitation, a time point at which both WT and Oxtr^1-/- animals display robust preference for partners over strangers⁵² (Fig. S1D). As with short cohabitation, the top quartile of huddles (by duration) for each subject animal following long cohabitation, described an average of 90% of total huddle time (Fig. S2A right). Following long cohabitation, preference and huddling behavior in Oxtr^1-/- females is indistinguishable from their WT siblings (Fig. 2D-E). Comparing changes in the frequency of top quartile huddles displayed by Oxtr^1-/- females after short vs long cohabitation, we find that both WT and Oxtr^1-/- females show decreased huddle frequency with stranger following long cohabitation (Fig. S2C-D). On comparing durations of long huddles, we find that Oxtr^1-/- females significantly increase their median huddle duration with partner post long cohousing (Fig. 2F-G). Thus, loss of Oxtr function in females reduces the duration of early prosocial interactions (short cohabitation) with partners and strangers while Oxtr-independent mechanisms eventually promote WT patterns of affiliative huddling and partner preference following longer periods of cohabitation.

Our results suggest that Oxtr function influences specific aspects of early pair bonding and attachment behaviors. We next wished to dissociate the influence of Oxtr on the formation and display of partner preference from the display of agonistic rejection behaviors towards novel strangers following cohabitation. We therefore developed a battery of assays to examine initial social interactions between mating partners at early and later stages of pair bonding, as well as stranger rejection following cohabitation (Fig.2A). We examined the behavior of subject females from each genotype with WT males when naive animals were first introduced to each other (Introduction), mating behaviors following estrus induction (Timed Mating), prosocial interactions following reunification with their bonded mate (Partner Reunification), and agonistic rejection of novel, opposite sex strangers (Stranger Rejection) (Fig. 2A see Methods). We find no differences in social behaviors during Introduction, Timed Mating and Partner Reunification (Fig. 2H-J). In contrast, during Stranger Rejection, Oxtr^1-/- females show more social behaviors than WT females towards novel naive males in the absence of their partner (Fig. 2K). Taken together, these studies demonstrate that Oxtr function regulates specific components of pair bonding in females in a state-dependent manner, suppressing promiscuous social interactions with stranger males after cohabitation sufficient for the display of partner preference.

Oxtr function influences different prosocial behaviors across the stages of pair bonding in males

Given the Oxtr-dependent difference in huddling behavior between males and females after short cohabitation (Fig. 1C-D), we wished to assess sex differences in the patterns of behaviors displayed by Oxtr^1-/- animals and to further investigate Oxtr function in males (Fig. 3A). Examining long huddles (as above, Fig. 2), we find that such huddles by Oxtr^1-/- males with their partners following short cohabitation are significantly shorter when compared to WT siblings, contributing to diminished partner preference early in pair bonding (Fig. 3B Left). Interestingly, while Oxtr^1-/- males show no difference in the total amount of time spent huddling with strangers when compared to WT siblings (Fig. 1D), they display shorter and more frequent huddles towards them at this time (Fig. 3B, 3C Right). In males, Oxtr thus influences the pattern as well as amount of early affiliative behavior to distinct behavioral stimuli in different ways.

Figure 3:

Oxtr^1-/- males show no difference in time spent with either partners or strangers after long cohabitation (Fig. S1E-F). However, these males show a persistent reduction in the duration of their long huddles with strangers at this time point (Fig. 3D). Following long cohabitation, WT and Oxtr^1-/- males show no difference in the duration of long huddles with partners and a reduction in the duration of long huddles with strangers (Fig. S3A-B). We find no difference in the frequency of long huddles between Oxtr^1-/- and WT males after long cohabitation (Fig. 3E). Consistently, comparing Oxtr^1-/- males’ behavior following short vs long cohabitation reveals that the frequency with which Oxtr^1-/- males huddle with strangers decreases following long cohabitation (Fig. 3F-G). Thus, specifically in males, loss of Oxtr function does not change the overall time spent huddling with strangers after short or long cohabitation but does cause a continuing disruption in the patterns of affiliative behavior that comprises huddling with strangers.

Our observations suggest that the loss of Oxtr disrupts certain characteristics of prosocial behaviors in males regardless of the time spent bonding with their partner. We next wished to dissect whether such loss of function resulted in atypical affiliative or agonistic behaviors across stages of bonding. We examined the behavior of Oxtr^1-/- males and their WT siblings in the same series of paradigms described previously to characterize the formation and maintenance of their pair bonds (Fig. 3A). Oxtr^1-/- males spend similar amounts of time as WT males engaging in prosocial behaviors with partners regardless of bonding state (Fig. 3H-J, S3E). Similar to our observations in Oxtr^1-/- females, Oxtr^1-/- males spend more time engaging in prosocial behaviors with strangers in the absence of their partners (Fig. 3K).

In contrast to time spent conducting mating behaviors with WT female partners (Fig. S3E), Oxtr^1-/- males demonstrate increased frequency of mating attempts when these females are in estrus (Fig. 3I). Examining all instances of mounting, we find that Oxtr^1-/- males carry out shorter mounting bouts more frequently than WT males (Fig. S3F). To determine whether the shorter bouts of mounting in Oxtr^1-/- males influence the receptivity of their partners, we compared the mount-to-intromission ratio for each animal⁵⁷ and find that females are equally receptive to Oxtr^1-/- males as their WT siblings (Fig. S3G). Thus, consistent with our previous observations, Oxtr influences patterns, but not the amount or success, of male mating behaviors, suggesting that such changes in males’ behaviors are equally able to elicit receptivity in females.

Compared to their WT siblings, naïve Oxtr^1-/- males show decreased agonistic behaviors towards an unfamiliar WT female upon first encountering her but display similar levels of rejection of novel (unfamiliar) strangers upon bonding (Fig. S3C-D). Thus, Oxtr function influences distinct components of males’ attachment-related social behaviors towards unfamiliar females at different points across pair-bonding, promoting agonistic displays and lengthening mating bout durations when naïve, and suppressing promiscuous prosocial behaviors upon bonding. Taken together our results suggest that Oxtr influences patterns of social behaviors with future partners during the early phases of pair-bonding differently in females and males, but upon bond formation, functions in both sexes to suppress promiscuous affiliative behaviors towards strangers (Fig. S3H) without affecting agonistic behaviors.

Oxtr sex-specifically influences social interactions between potential mates

Oxtr function influences patterns of social interactions between mates early during the formation of a pair-bond (Short Cohabitation in both sexes, Introduction and Timed Mating in males). We therefore investigated the effect of Oxtr loss on the choice made by a potential mate. We developed a variation of the PPT three-chamber assay to present a naive WT “chooser” with a choice between a naive opposite sex WT or naive opposite sex Oxtr^1-/- animal (Fig. 4A, Fig. S4A, see Methods). Strikingly, we find that Oxtr influences the choice by a WT mate in a sex-specific manner. WT females choose WT or Oxtr^1-/- males with equal frequency (Fig. 4B). In contrast, WT male choosers show significant preference for WT females compared to Oxtr^1-/- females, particularly during the middle of the assay (150-280 mins) (Fig. 4B-C, Fig. S4B). These observations suggest that disruption of Oxtr function specifically in females may alter their behavior such that WT males prefer WT future partners over such mutants.

Figure 4:

While naïve WT females show no difference in overall social choice between naive WT and Oxtr^1-/- males (WT-Oxtr^1-/-, Fig. 4B), they may show distinct behavioral dynamics in the presence of an Oxtr^1-/- male. To assess this, we compared the WT-Oxtr^1-/- condition with a control group in which the naive female chooser is given access to two naive WT males (WT-WT) (Fig. 4D). We designated the tethered animal with whom the subject spent More of Assay Time (“chosen”) as the MOAT. Intriguingly, we found that the presence of an Oxtr^1-/- male (WT-Oxtr^1-/-) influenced the behavior of the WT female chooser, such that they spent significantly more time with WT or Oxtr^1-/- MOAT males compared to the time spent with WT MOAT males from the WT-WT condition (Fig. 4E). Consistent with this finding, naive WT female choosers break contact more frequently in WT-WT conditions when compared to WT-Oxtr^1-/- (Fig. 4F). Similarly, the divergence in the amount of cumulative time spent with the MOAT for WT female choosers occurs after 5 hours have elapsed (320-360 mins) suggesting that complex interactions with males lacking Oxtr influence how females “choose” a male partner regardless of whether they end up choosing a WT or Oxtr^1-/- MOAT (Fig. 4G). In contrast, we find no difference in the contact time or number of such transitions when comparing naive WT male choosers with MOAT females of a WT-WT vs WT-Oxtr^1-/- condition (Fig. 4E-F).

We next examined the pattern of a WT chooser’s interactions with both chosen (MOAT) and unchosen animals (Less of Assay Time, LOAT) over time using a sliding preference index (see Methods). We observed periods of convergence when choosers from a condition (WT-WT or WT-Oxtr^1-/-) displayed significantly concurrent interactions with their MOAT, indicated by the black bar above the trace. WT female choosers show a longer and more consistent preference for their MOAT (regardless of genotype) in WT-Oxtr^1-/- vs WT-WT conditions, whereas WT male choosers show an earlier onset but equivalent length of preference for MOATs in the WT-Oxtr^1-/- condition compared with the WT-WT condition (Fig. 4H). These observations suggest that interactions with males lacking Oxtr influence the behavior of the WT females making a social choice, even when these females interact more with and eventually choose a WT male. Taken together, loss of Oxtr leads to differences in prosocial behaviors during initial interactions with WT naive animals of the opposite sex that may later impact aspects of pair bonding.

Finally, to determine if naive prairie voles lacking Oxtr display generalized changes to their behaviors independent of attachment-associated social behavior, we examined general pro-social, anxiety-related, and locomotor behaviors. With the exception of an increase in locomotion displayed by Oxtr^1-/- females when compared to WT females in an open field paradigm, we found no differences between Oxtr^1-/- and WT voles of either sex (Fig. S4C).

We have examined social interactions during the initial contact between mates, the early stages of pair bonding, and during social attachment behaviors following partner preference formation. Consistent with our findings and previous work, Oxtr function is not required for the eventual display of partner preference. Here, we demonstrate that Oxtr function sex-specifically facilitates early social interactions that contribute to the formation of partner preference, and controls promiscuity.

Oxtr controls the molecular signature of pair bonding in the nucleus accumbens

We previously observed that Oxtr function is not required for the formation of partner preference and here demonstrate that it suppresses prosocial behaviors towards novel strangers following pair bonding. These results suggest that distinct components of pair bonding are differentially influenced by Oxtr function or compensated upon its loss. We next characterized the molecular signature associated with pair bonding to determine whether changes in gene expression following pair bond formation are controlled by Oxtr. Comparative studies demonstrate that Oxt binding in multiple brain regions differs significantly between promiscuous and monogamous rodent species. Oxt binding is particularly enriched in the NAc in prairie voles, while vasopressin binding is largely absent in this area^45,47. Thus, analyzing gene expression in the NAc allows us to identify bonding-related expression that is primarily due to Oxtr function and largely independent of the direct consequences of AVP signaling. Additionally, given the well-established role of the NAc in the association of reward or aversion with specific sensory and environmental cues^58–60 this region is well poised to play a role in the development of partner preference, even without Oxtr signaling, as well as stranger rejection, which appears to be dependent upon Oxtr function.

We performed bulk RNA-sequencing on NAc tissue from WT and Oxtr^1-/- animals of both sexes that were either: group housed with same-sex peers and sexually naïve (pre-pairing); or paired with WT opposite-sex partners for six days, using our timed mating paradigm (post-pairing) (Fig. 5A). PCA showed robust separation of WT and Oxtr^1-/- samples (Fig. S5A). Differential expression (DE) analysis of pairing condition, regardless of genotype or sex, identified 13 genes that show significant changes in expression post-pairing (Fig. S5B), all of which showed decreased expression (Fig. 5B). These genes include Agt, Sgk1, Ptgds, Fosb, Pdk4, and Acsm5 (Fig. 5C) and, while these genes have not previously been linked to shared gene regulatory or functional pathways, the predicted protein-protein interaction network for the 13 DE genes shows more connections than would be expected by chance (p=2.47 e^-05) (Fig. S5C). Intriguingly, even in the limited set of genes with altered expression in the context of pair bonding, changes are driven primarily by expression in WT females (Fig. 5C, S5D). Subgroup analyses by genetic background and sex identified additional genes including Adora2a, Fkbp5, and Spred3 in WT females and collagen-related genes, including Col4a1 and Col4a2, in Oxtr^1-/- females (Fig. 5D, S5D). Only one gene, Sema3b, is significantly different in Oxtr^1-/- males following pairing (Fig. 5D).

Figure 5:

To further examine changes in patterns of gene expression in the NAc associated with pair bonding, we performed weighted gene coexpression network analysis (WGCNA)⁶¹ and calculated three sets of modules (M) from different sample subsets (modules^WT, modules^Mut, and modules^All) (Fig. S5F, S6D). Using the modules^WT subset, we identified 15 coexpression modules (Fig. S5F, Table 2) representing genes that share highly similar expression patterns within the NAc. We identified one module (Mpink^WT) that significantly and specifically correlated with pre-pairing status in WT animals (Pearson R correlation of 0.7, Bonferroni adjusted p-value (p_adj)=0.018) (Fig. 5E, S5E-F). Genes with high module membership for the Mpink^WT are also those with high gene significance for pairing condition supporting the module’s association with behavior condition (Fig. 5SG). Further, genes significantly differentially expressed with bonding are significantly enriched in Mpink^WT (Bonferroni p_adj=7.44 e^-08) and serve as hub genes within this module (Fig. S5E, G). These results suggest that the Mpink^WT module might be driving, or under the influence of pathways driving, the pre-pairing to post-pairing state change.

We asked whether any of the modules^WT identified were regulated by Oxtr and thus lost in Oxtr null animals. To address this, we assessed how well modules^WT are preserved in mutant samples. We calculated a Z_summary measure that combines module density and intramodular connectivity metrics into a composite statistic. We find that of the modules^WT, Mpink^WT had the lowest Z_summary (Z_summary = 2) in Oxtr^1-/- samples (Fig. 5F), indicating that dynamics of gene expression that correlate with bonding states are lost in the absence of Oxtr. Finally, the module eigengene (ME) for Mpink^WT when calculated for all 30 samples and compared between pre- and post-pairing, was significant only for WT females (Wilcoxon rank sum, Benjamini-Hochberg p_adj=0.03) and not for WT males or Oxtr^1-/- animals of either sex (Fig. 5G). These observations suggest there is sex-specific, coordinated regulation of gene expression in the NAc following pairing of mates that may be associated with pair bond associated behaviors.

Loss of Oxtr sex-specifically disrupts patterns of gene expression associated with neurodevelopmental disorders

We next wished to determine how Oxtr signaling contributes to changes in gene expression in the NAc associated with the naive state or the process of pair bonding (Fig. 6A) by comparing the transcriptional profile of WT and Oxtr^1-/- animals of each sex. We observed significant genotype specific differences in gene expression in animals from both sexes, with 466 genes showing differences between WT and Oxtr^1-/- females, 270 between WT and Oxtr^1-/- males, and 1,014 genes when samples from both sexes were combined (using a p_adj<0.05 and absolute (log₂(FC)) > 0.25 as criteria). Of the combined DE gene set, 370 genes are upregulated and 644 are downregulated (Fig. 6B, S6A-B). Genes with the greatest fold change across both sexes include Slc2a3, Kcna1, Rgs8, Cnksr2, Acvr, and Cnot6 (Fig. 6B). Genes specifically downregulated in Oxtr^1-/- females with log₂(FC) less than 1 include Abcb11 (a gene recently associated with treatment resistant schizophrenia⁶²), Ngfr, Pon3, and Shh. While there was no significant difference in Oxtr transcript levels by genotype (p_adj = 0.753)—consistent with minimal nonsense-mediated decay despite a premature stop codon—we have previously shown that no functional protein is produced in Oxtr^1-/- animals⁵². V1aR showed a modest decrease in expression in Oxtr^1-/- animals (p_adj=0.03) (Fig. S6C), though this is likely influenced by variability in WT samples. WGCNA of all samples identified two modules (Mbrown^All and Mturquoise^All) that were significantly positively correlated with the WT genotype compared to mutant (Pearson R = 0.55, p_adj = 0.011 for brown and Pearson R = 0.8, p_adj = 6.6e-7 for turquoise) (Fig. S6D). None of the identified modules showed correlation with pairing condition (Fig. S6D). Our results thus reveal a dramatic change in patterns of gene expression in the NAc with global loss of Oxtr signaling, independent of bonding state.

Figure 6:

Examination of the sets of genes whose expression is altered by the loss of Oxtr, regardless of bonding state or sex, revealed several distinct cellular processes regulated by Oxtr signaling. These processes comprise larger domains of cell function including extracellular channel binding activity, axoneme assembly complex, axon glial ensheathment development, and metal cation transport (Fig. 6C), as well as pathways for G protein signaling and dopaminergic neurogenesis (Fig. S6E). When analyzed by sex, more of the identified processes were enriched in female Oxtr gene sets than in male ones (female: 19/62, male: 6/62). To further link differentially expressed genes to behavioral traits, we asked whether genes altered by loss of Oxtr are selectively enriched for high-confidence associations with autism spectrum disorders (ASD) and other neurodevelopmental diseases (NDD) affecting social attachment in humans as compared to epilepsy, a disorder that does not disrupt social attachment^63–66. Consistently, we find that the genes upregulated with loss of Oxtr (p_adj<0.01, log₂(FC)>0) are specifically enriched for ASD- and NDD-associations but not epilepsy associations (Fig. 6D “Joint”). This enrichment is again primarily driven by genes upregulated in females and not males (Fig. 6D). Thus, similar to the factors that influence gene expression changes associated with pair bonding, many of the changes we identify in the NAc resulting from loss of Oxtr are primarily driven by female-specific effects that may play a conserved role in mediating social attachment in humans.

Oxtr regulates coordinated gene expression across multiple neural populations implicated in social attachment behaviors

Given the expression of Oxtr in regions of the brain known to mediate social behavior in voles and other species^67,68, we wished to determine if genes whose expression is altered in the NAc by the loss of Oxtr were expressed not only in specific populations of cells within this region, but in other regions of the prairie vole brain implicated in social and attachment behaviors. We therefore examined the expression of 15 genes that showed significant changes in expression in WT females when compared to those lacking Oxtr using in situ hybridization (ISH) to both validate our RNA-seq findings as well as examine the expression patterns of these genes. Several of these genes, such as phosphodiesterase 10A (Pde10A) and diacylglycerol kinase beta (Dgkb), show robust expression throughout the striatum (data not shown), while genes such as the calcitonin receptor (Calcr) and delta like non-canonical notch ligand 1 (Dlk-1) are expressed in restricted populations of cells in the NAc (Fig. 6F-H, top panels). The expression of these genes is significantly reduced in Oxtr^1-/- females (Fig. 6E-H), supporting the specificity and resolution of our RNA-seq findings.

In parallel to our examination of gene expression in the NAc, our ISH revealed that loss of Oxtr in females altered the expression of the genes we examined in other regions of the prairie vole brain. We examined the expression of these genes in the lateral septum (LS), a region that has been implicated in agonistic and attachment behaviors in prairie voles and other rodents and shows species dimorphic patterns of Oxt and Avp binding. We observe that both Calcr and Dlk-1 are decreased in the absence of Oxtr (Fig. 6G, H middle panel) in the LS. We also find that expression of Calcr and Dlk-1 is decreased in the paraventricular nucleus of the hypothalamus (PVN), a principal source of Oxt and Avp in the rodent brain which sends projections to the NAc (Fig. 6G-H bottom panels)^69–71. These findings suggest that Oxtr function regulates patterns of gene expression in multiple regions of the brain implicated in social and attachment behaviors in a coordinated and region-specific manner.

Loss of Oxtr results in the decrease in density of Oxt and Avp expressing neurons in the PVN

We find here that loss of Oxtr results in coordinated changes in gene expression in multiple neural populations, including the PVN. We therefore wished to determine if the signatures of gene expression regulated by Oxtr in the NAc in prairie voles (ONP set), regardless of sex and bonding state, were enriched for patterns of gene expression associated with distinct neural populations in other regions of the brain, particularly those implicated in social and attachment behaviors. We compared the ONP gene set to the molecular signatures of 33 hypothalamic subtypes previously identified by single-cell RNA sequencing in mice⁷². Of the 33 clusters examined, the four with highest overlap with the ONP gene set represent subsets of PVN neurons, in particular those that express Oxt and AVP (Fig. 7A, clusters 43 and 26), as well as PVN neurons that express Dlk1 and Calcr (cluster 29), both of which show reduced expression in the NAc and PVN in females lacking Oxtr (Fig. 6G-H bottom panel).

Figure 7:

We next sought to determine if the overlap between Oxtr regulated genes and gene expression patterns of neural populations in the PVN correlated with changes in the expression of neuropeptides within PVN itself. We examined patterns of Oxt and Avp expression using immunohistochemistry in naive adults and compared patterns of expression between WT animals of both sexes and those lacking Oxtr (Fig. 7B, Fig. S7A-B). We find that the density of both Oxt and Avp producing neurons is lower in Oxtr^1-/- animals (Fig. 7C-E). To determine if the loss of Oxt or Avp expression was restricted to specific populations of PVN neurons, we examined the expression of these peptides across the anterior-posterior axis of the PVN in both WT and Oxtr^1-/- males and females. We find that animals of both sexes show a preferential loss of Oxt and Avp expressing neurons in the anterior PVN (Fig. 7F-G). Finally, we examined patterns of Oxt and Avp over the course of development to determine when loss of Oxtr function influences the development of their expression. We find that the density of Oxt neurons at P0 (Fig. S7C-F) and both Oxt and Avp neurons at P21 is lower in Oxtr^1-/- females compared to age matched WT females (Fig. S7G-J). In contrast, we find that Oxt and Avp neuron density loss occurs only in adult Oxtr^1-/- males (Fig. 7D-G, Fig. S7B). These observations suggest the coordinated, region-, and sex-specific influence of Oxtr function on the development of neurons expressing Oxt (and Avp) within the PVN, consistent with a trophic role for this pathway identified in mice^73–76.

Discussion

Robust reproductive behaviors are essential for survival. Accordingly, species have evolved diverse strategies, including social monogamy, which mediate locating, identifying, and attracting mates, facilitate reproduction, and promote the survival of progeny^16,77. Signaling via Oxtr influences a wide range of physiology and behavior associated with social interactions and attachment between conspecifics, including pair bonding between prairie vole mates^41,78. Our recent findings that Oxtr is not genetically required for the formation of partner preference in prairie voles suggests that, like other genes and neuromodulatory pathways across diverse species, Oxtr likely functions in parallel to multiple genetically specified pathways to facilitate pair bonding by influencing specific behavioral modules that contribute to the complex behaviors that comprise social attachment^{33,35,79–85}.

Here we demonstrate that loss of Oxtr in prairie voles influences the dynamics of social interactions between sexually naïve animals and WT potential mates in distinct social contexts in sex-specific ways. When given a choice, naive WT males display a preference for WT females rather than for those lacking Oxtr, supporting the role of Oxtr in controlling behavior during an early, sensitive period in bond formation³⁰. The presence of an Oxtr^1-/- male alters levels and patterns of early prosocial behavior exhibited by WT females, without ultimately affecting their choice of potential partner. Taken together, these observations suggest that Oxtr may modulate social behaviors that act early in a sensitive pair bond formation period, and via these behaviors or yet-to-be-determined sensory cues such as pheromone profiles or patterns of vocal communication, influences the reciprocal behavior of partners.

Animals lacking Oxtr take longer to display partner preference, and despite strong partner preference following longer cohabitation, demonstrate promiscuous prosocial behaviors towards novel potential mates, including huddling with strangers. These observations are consistent with a state-dependent role of Oxtr function, similar to findings in multiple species, including humans^86–88. Oxtr function may thus promote prosocial behavior in the naïve state to facilitate attachment with an unfamiliar potential mate, but following pair bond formation, suppresses further prosocial interactions to new individuals, or even promotes further separation from and rejection of unfamiliar individuals or groups. Across species and outside of the context of pair bonding, oxytocin signaling may play a conserved role in mediating other state-specific aversive behavior, for example, defeat-induced social avoidance learning⁸⁹.

Our findings reveal that patterns of social behaviors contributing to 1) pair bond formation, 2) partner preference, and 3) rejection of novel mates are genetically separable components of pair bonding representing distinct behavioral modules regulated by mechanisms that require Oxtr function, i.e. reciprocal behaviors and stranger rejection, or those that can be displayed in its absence, i.e. partner preference. Further work will allow us to determine the Oxtr-independent mechanisms that mediate the closure of an early sensitive period and facilitate the display of partner preference⁵². Where Oxtr functions to control the switch from prosocial behavior towards a potential mate in naïve animals to the suppression of such behaviors towards, and active rejection of, opposite-sex conspecific strangers following bonding remains to be determined.

Examination of the molecular signatures of pair bonding and the role of Oxtr in regulating these changes in gene expression in the NAc reveals modest, Oxtr-dependent changes in the expression of a specific network of genes in WT animals following pair bonding^45,90. Intriguingly, loss of Oxtr preferentially impacts patterns of gene expression in the NAc in females. The NAc demonstrates little endogenous sex steroid hormone receptor expression difference in gene expression in prairie voles and other species^67,91–93. Sex differences in the effects of Oxtr on gene expression this region may therefore suggest that there are significant sex-differences in the inputs into this region that act through Oxtr during development^94–99. The larger impact of loss of Oxtr in females is consistent with its more general conservation in mammals to regulate aspects of female physiology such as lactation and prosocial behaviors or sensory cues that promote interactions with males¹⁰⁰.

Genes with altered expression in the absence of Oxtr are enriched for genes associated with neuropsychiatric disorders, in particular neurodevelopmental disorders and ASD. Given the sex-biased enrichment of these genes, such developmental differences in the impacts of Oxtr function may contribute to sex differences in the pathophysiology of these disorders^101–103. Accordingly, the molecular signatures of attachment in the NAc, both dependent upon and independent of Oxtr function, are subtle relative to these sex-biased differences in the role of Oxtr and possibly reside in other neural populations that function in attachment.

Validation of changes in gene expression in the NAc resulting from loss of Oxtr suggests coordinated regulation of these genes by Oxtr across multiple brain regions implicated in social and attachment behaviors. Oxtr may be one of multiple factors that organizes specific components of the circuits underlying these displays to control discreet modules of behavior. Consistent with a trophic role for Oxtr function during development, we find that the genes regulated by Oxtr within the NAc are enriched for those that identify cell types in the PVN, a principal source of neuropeptides implicated in pair bonding, which sends oxytocinergic projections to the NAc^{68,70,104–106}. Loss of Oxtr preferentially impacts Oxt and Avp expression in the anterior PVN, suggesting that specific populations of these cells may be more dependent upon such trophic function for their development, survival, or expression of these hormones. Whole brain analysis will likely reveal if such neurons have distinct projections when compared to those that remain in the absence of Oxtr, contributing to our understanding of how these circuits and pathways regulate distinct aspects of social and attachment behaviors and physiology^4,78,107.

Taken together, our behavioral and molecular studies reveal previously unappreciated nuances in the control of social behavior and pair bonding by Oxtr function. They uncover intriguing sex differences in the patterns of behavior influenced by loss of Oxtr function. These observations suggest that circuits with sex-specific patterns of development, or even function in ancestral species, may mechanistically converge in the context of Oxtr function to generate displays of reciprocal interactions, pair bonding and attachment behaviors in WT prairie voles that are similar between males and females^108,109. Manipulation of Oxtr expression at specific developmental timepoints and in specific neural populations will help to elucidate how this pathway influences social and attachment behaviors in the context of normal development to regulate distinct components of pair bonding. Understanding the specific neural populations that express Oxtr and other molecular pathways that contribute to pair bonding will allow us to understand the mechanisms controlling the early experience-dependent, sensitive period of bond formation. In parallel, examination of the activity in these populations following manipulations of Oxtr function will allow us to determine the circuit mechanisms by which similar but distinct sensory cues contribute to the identification of partners or novel conspecifics and consequently, the display of dramatically different patterns of social behaviors. The sum of these mechanisms represents the preference for partners over strangers and reflects the formation and maintenance of pair bonds between mates, affiliation between individuals, and the behavioral consequences of identifying others as unfamiliar or outside of these bonds.

Materials and methods

Star Methods

Key resources table

Resource availability Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Devanand Manoli (devanand.manoli@ucsf.edu).

Materials availability

All unique reagents generated in this study are available from the lead contact upon request.

Data

Transcriptomics data discussed in this publication are accessible through GEO Series accession number GSE279248.

Code

All code used for RNA-seq data processing and differential expression analyses, as well as post-processed expression data, can be found at https://github.com/aseveritt/OXTR_prairie_voles

Experimental model and subject details

Animals

Subjects were laboratory-bred prairie voles (Microtus ochrogaster) which originated through systematic outbreeding of a wild stock captured near Champaign, Illinois. CRISPR/Cas9 editing was used to generate 2 independent null alleles, Oxtr^1-/- and Oxtr^2-/-, both of which were validated using autoradiography to show loss of receptor function⁵² (data not shown). Oxtr^1-/- was backcrossed for >7 generations to mitigate any possible off-target effects and this line was used for all subsequent experiments.

Sexually naïve male and female animals were group weaned at 21 ± 1 days and separated to group housing with same-sex siblings and age-matched same-sex non-siblings. Voles were maintained under a 14:10 h light-dark cycle in clear plastic cages (45 × 25 × 15 cm) with bedding, nesting material (nestlet), and a PVC hiding tube.

Rooms were maintained at approximately 20°C, and food and water were available ad libitum.

Breeding pairs were established between two heterozygotes or a homozygous Oxtr^1-/- mutant male and heterozygous female partner from a breeding line maintained in the lab. Voles were randomly assigned into experimental groups when they reached 7-9 weeks of age at the start of testing. This study was carried out in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals published by the National Research Council. The protocol was approved by the Institutional Animal Care and Use Committee at each respective institution.

Genotyping

The following primers were used to amplify a region where XcmI restriction for Oxtr¹ or Bpu10I restriction for Oxtr² determined the genotype of all animals studied. XcmI cuts the mutant allele in the case of Oxtr¹ and Bpu10I restricts the WT sequence for Oxtr². (Forward: ACTGGAGCTTCGAGTTGGAC, Reverse: ATGCCCACCACTTGCAAGTA)

Method details

Behavior

Triadic naive choice assay

The linear three chamber partner preference apparatus with open top and 10 x 32 inch walls⁵² was modified to capture the interactions of the naive, freely moving WT subject and the naive, tethered, opposite sex potential partner (either 2 WTs or a WT and an Oxtr^1-/-) at either end of the chamber. Behaviors were recorded for 6 hours using a top view camera and the videos generated were analyzed using a contour detection algorithm. Specifically, contours were detected on individual video frames using a threshold-based strategy that took advantage of the contrast between dark-colored vole fur coats and light-colored bedding and chamber apparatus. A standalone vole was detected as a single contour, and two voles in physical contact were likewise detected as a single contour. Consistent vole identities across video frames were determined using the Hungarian algorithm. Since the assay began with one vole in each of the three chambers, and since the subject vole remained freely mobile throughout, it was possible to unambiguously assign each contour to a single vole identity (e.g., subject, right chamber, or left chamber) or to that of a pair of voles (e.g., subject and right chamber voles or subject and left chamber vole), resulting in continuous trajectories for each of the three animals. Time in chamber was determined by the chamber location of the subject vole in each frame, since the chamber location of the right chamber and left chamber animals was fixed throughout the assay by physical tethering and any time the contour of the subject animal merged with that of the tethered animal (right or left), the subject was termed in contact with it. The algorithm also determined the boundaries of the chambers and calculated chamber times. We used the chamber time data from when both ends of the chamber contained tethered WTs (WT-WT condition) to ensure that there was no bias in which side the choosers preferred.

Bond formation assays

Initial introduction of Oxtr^1-/- and WT animals of either sex with an opposite sex WT animal was filmed for 30 minutes. Observed behaviors included anogenital sniffing, investigation (subject animal’s nose contacts any non-genital region), duration of stationary huddling or >50% side-by-side contact between animals, rearing, frequency of aggressive behavior including lunges, receipt of aggression, and tussles, as well as mating behavior including mounting, intromission, and ejaculation. These animals remained in their home cage for 18 hours following pairing, after which, they were separated by a clear, plastic barrier. The barrier was left in place for 24 hours and then removed allowing the subject free access to the partner animal thereafter. A 30-minute period directly following barrier removal was also filmed to assess mating behavior, termed the timed mating assay.

Subjects remained cohoused with the opposite-sex, wildtype partner for durations as indicated in the text. The partner preference test (PPT), as previously published⁵², was used to assess pair bond formation. Subject, partner, and an opposite sex stranger are placed in a three-chamber apparatus with open top and 10 x 32-inch walls. The partner and stranger animals are tethered on either end of the three-chamber arena and the subject allowed free access over a period of 3 hours. Behaviors are recorded from a top view camera capturing the entire apparatus. Videos were scored using BORIS¹¹⁰ post-test by validated scorers blind to condition. Observed behaviors included location (i.e., duration of time in partner, stranger, and neutral chambers), duration of stationary huddling or >50% side-by-side contact between the partner and stranger animals, and frequency of aggressive behavior (i.e., lunges). Of note, analysis of partner preference between wild type siblings from the Oxtr^1-/- background and wild type animals from the breeding colony showed no difference in time spent huddling with partner vs stranger animals.

Bond maintenance assays

Four days after the timed mating assay, female animals underwent a separation-reunification assay. Females remained in their home cage while their male partners were placed in a separate cage in another room for one hour. The males were then reintroduced into the home cage and behavior was video recorded for 30 minutes. Male subject animals underwent the same procedure five days following timed mating. Two days following the separation reunification test, animals underwent a selective aggression assay. Partners were again removed from the home cage and separated to a separate room for an hour. After one hour, a novel sexually naïve animal of the opposite sex was introduced into the home cage with the subject animal and behavior observed for 20 minutes. If significant aggression occurred, defined as >3 tussles during the assay, the animals were separated and the assay terminated. Overall, two WT female assays and one Oxtr^1-/- were terminated for aggression. The same behaviors listed above for the introduction were observed and recorded in the reunification and aggression assays.

RNA Sequencing

Sample collection

Adult 7-9 week old WT and Oxtr^1-/- voles were either cohoused with same sex-siblings or paired with an animal of the opposite sex for four days following the timed mating protocol described above. These animals were euthanized by CO₂ inhalation and decapitated. The brain was quickly dissected out and sectioned into 500μm coronal slices using a brain matrix mold (BrainTree Scientific) chilled on ice. Slices were floated in ice-cold phosphate buffered saline (PBS), the NAc was identified using anatomical landmarks, and dissected using a Zeiss microscope. A total of 45 WT females, 39 WT males, 21 Oxtr^1-/- females, and 18 Oxtr^1-/- males were split between paired and unpaired conditions and tissue from three animals per condition was pooled prior to RNA extraction. Tissue was flash-frozen in liquid nitrogen and stored at -80°C until further processing. Total RNA was extracted using TRIzol according to manufacturer’s instructions and quantified using NanoDrop (Thermo Scientific). Library preparation for next generation sequencing, comprising 35 samples and 4 technical replicates, was done using TruSeq Stranded mRNA Kit according to manufacturer recommendations and sequenced on a NovaSeq 6000 to an average depth of 2.9e07 reads per sample. Sequencing was performed across three runs, but all data was processed together.

Read processing

Adapter and low-quality sequence trimming was performed with Trimmomatic v0.39 using parameters: ‘TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:25’ (Bolger et al., 2014). Trimmed RNA-Seq reads were aligned to a custom MicOch1.0 prairie vole genome, which incorporated a known missing gene V1A, using STAR v2.7.3a in gene annotation mode (Dobin et al., 2013). Alignment, RNA-Seq, and Insert Size quality control metrics were generated using Picard v 2.10.10 (http://broadinstitute.github.io/picard). Homologous genes between vole, mouse, and human references were identified using biomart v2.56.1. When multiple mouse or human genes mapped to the same vole gene, the homologous gene with the highest percent identity to vole was selected.

Differential expression analysis

To remove non- and lowly-expressed genes, only genes with more than 2 counts per million (cpm) in at least 3 samples were retained. Principal component analysis of the cpm revealed that technical replicates all cluster tightly but showed a slight batch effect due to sequencing run (Fig. S5A). To remove this effect, we used RUVSeq v1.34.0 to estimate one factor of unwanted variation, W_1, which we later used to normalize the data prior to differential expression¹¹¹. The factor analysis was performed on the deviance residuals from an initial generalized linear model regression of the upper quartile normalized counts on our covariates of interest (pair bonding status, genotype, sex). Expressed genes were tested for differential expression (DE), after combining technical replicates, with DESeq2 v1.40.2 using a design matrix that included the covariates and W_1¹¹². The adaptive shrinkage estimator ashr v2.2-63 was used in order to be more robust to genes with low cpm values¹¹³. For DE analyses in general, while most changes were driven by female samples, the effect was maintained in the males at a smaller magnitude; thus, analyses were primarily performed for both sexes jointly to maximize accuracy and power.

Bonding effect

To compare pre- to post-pairing individuals, we performed the DE analysis in two ways. First, we compared pairing status regardless of genotype or sex to increase our statistical power and to identify a reproducible gene set of thirteen genes (Fig. 5B, 5C, S5B). While the relationship between the genes is unknown, there is support that the proteins are functionally related. Using the STRING protein database (STRINGdb_2.12.1), the human-homologs of the DE genes share more connections than we would expect by chance (Fig. S5C, p-value= 3.08E-0.7). Next, we compared pairing status using all three covariates, to identify subgroup-specific genes of interest and effect sizes (Fig. 5D, S5D). In both analyses, genes were considered significantly DE with a BH adjusted p-value < 0.05.

Genotype effect

To examine the effect of genotype, we first performed DE analyses between WT and Oxtr^-/- individuals, not including sex as a covariate to identify a reproducible and stringent gene set. We considered genes DE with an adjusted p-value < 0.05 and absolute log2 fold-change > 0.25 (Fig. 6B, S6A, S6B). To examine the functional relationship of the DE genes, we performed a Gene Ontology enrichment analysis using goseq v1.52.0 which controls for gene length bias¹¹⁴. We considered all expressed vole genes with a mouse homolog as our background and all DE genes with a mouse homolog as our foreground. GO categories with an adjusted p-value < 0.05 were considered significantly enriched. To help visualize the results, enriched GO categories were clustered according to their semantic similarity calculated with the R package GOSemSim v2.26.1 using method “Wang”¹¹⁵. Large, general, GO categories with over 1000 genes were excluded in visualizations (Fig. 6C). WikiPathways¹¹⁶ over-representation analysis was performed with R package clusterProfiler v4.8.3¹¹⁷ (Fig. S6D).

For disorder and disease set overrepresentation analysis, Chi-squared tests were used to assess whether Oxtr^1-/- disrupted genes occur in disease-associated gene lists at a frequency higher than we would expect by chance. Genes lists for autism spectrum disorder⁶³, schizophrenia⁶⁴, neurodevelopmental disorders⁶⁶, and epilepsy⁶⁵ were sourced and all alternate or outdated gene symbols were updated manually. For this test, we were interested in sex- and direction-specific enrichments, so we performed a DE analysis separately for males and females. Across all DE analyses, we defined DE genes more stringently (adjusted p-value < 0.01). For enrichment tests, our universe was all expressed vole genes with a human homolog. We split this analysis by fold-change direction to highlight the distinct enrichment in upregulated genes. P-values underwent BH correction across all 24 comparisons.

To compare to mouse neuronal populations, we leveraged single-cell data on mouse hypothalamus development⁷². The released Seurat object GSE132730_subset_cca_to_velo.rds.gz was accessed at GSE132730. Using the authors’ defined clusters (accessed at “prim_walktrap”), we assessed how similar each clusters’ gene expression is to the genes disrupted by the Oxtr^1-/- mutation. This was calculated as the Jaccard similarity coefficient using R package clustifyR v1.12.0¹¹⁸.

Weighted Gene Coexpression Network Analysis

Data

Starting with the RNAseq dataset of n=30 samples (one WT sample and 4 replicate samples were removed) and n=13357 genes, we removed genes with low expression (n=1738 genes with expression level in lowest quantile for all 30 samples) and/or low variability (n=8230 genes with coefficient of variation < 0.05), resulting in a trimmed dataset of n=3389 genes.

Identifying co-expression modules

We used the R package WGCNA (v1.70-3) to conduct WGCNA analyses⁶¹. Briefly, we calculated Pearson correlations for all gene pairs. Next, we constructed an unsigned weighted correlation network to identify co-expression modules comprised of highly correlated genes (positive or negative correlation) with high topological overlap¹¹⁹. Modules were defined as branches of a hierarchical cluster tree using the top-down dynamic tree cut method⁶¹. For each module, the expression pattern was summarized by the module eigengene (ME), which is defined as the right singular vector of the expression patterns. Pairs of modules with high module eigengene correlations (r > 0.75) were merged to maintain a level of decorrelation among MEs.

In more detail, a weighted unsigned network was computed based on a fit to scale-free topology, and the lowest thresholding power that resulted in a scale-free R² fit of 0.9 was selected, and the pairwise topological overlap (TO) between genes was calculated¹¹⁹. We constructed a TO dendrogram using hierarchical clustering (R::hclust function, method= ”average”), and defined modules using the WGCNA::cutreeHybrid function with minimum module size set to 30 genes. We then merged modules with high module eigengene correlations as describe above.

We calculated 3 sets of modules from different sample subsets:

1. Modules^WT: n=15 modules constructed from n=18 WT samples with thresholding power β = 5.

2. Modules^Mut: n= 43 modules constructed from n=12 Oxtr¹-mutant samples with thresholding power β = 6.

3. Modules^All: n=7 modules constructed from all n=30 samples with thresholding power β = 3

Relating modules to sample traits

In order to correlate modules with experimental conditions, we calculated the Pearson correlations between module eigengenes and sample metrics, including species (Mut=0, WT=1), condition (Pre=1, Post=0), sex (F=1, M=0), and quality control (QC) metrics (RIN, yield, sequencing run). We identified three modules that significantly correlated with a sample metric of interest but not with any of the QC metrics.

Significant correlation was defined to be p_adj <= 0.05, where p_adj = BH correction by number of modules. These modules (M) were:

1. Mpink^WT: positive correlation with pre-bonding status (Pearson R = 0.7, p_adj = 0.018)

2. Mbrown^All: positive correlation with WT status (Pearson R = 0.55, p_adj = 0.011)

3. Mturquoise^All: positive correlation with WT status (Pearson R = 0.8, p_adj = 6.6e-7)

We calculated the module eigengenes (ME) for Mpink^WT, Mbrown^All, and Mturquoise^All for all n=30 samples using the WGCNA::moduleEigengenes function. We then conducted Wilcoxon rank sum tests (ggpubr::stat_compare_means function) to assess whether MEs were significantly different between the following groups:

1. MEpink^WT: pre- versus post-bonded animals, separated by sex and species

2. Mbrown^All: WT versus mutant animals, separated by sex and condition

3. Mturquoise^All: WT versus mutant animals, separated by sex and condition

Comparing gene significance with module membership

We assessed whether genes whose expression highly correlate with module eigengenes of modules of interest (Mpink^WT, Mbrown^All, and Mturquoise^All) are themselves highly correlated with traits of interest (bonding status or mutation status). Specifically, for the genes in Mpink^WT, Mbrown^All, and Mturquoise^All, we calculated 1) gene significance (GS), defined as the absolute Pearson’s correlation between expression of a given gene and a trait of interest; and 2) module membership (MM), which is defined as the absolute Pearson’s correlation between the expression profile of a gene with the ME of a module to quantify the relationship between a gene and a given module⁶¹. We then calculated the Pearson correlations between:

1. GS^condition and MMpink^WT for Mpink^WT genes (R = 0.69, p = 5.5e-18)

2. GS^species and MMbrown^All for Mbrown^All genes (R = 0.63, p = 2.2e-32)

3. GS^species and MMturquoise^All for MMturquoise^All genes (R= 0.82, p<1e-200) Hub genes were defined as those genes in the top 25% for both the absolute value of pink module membership and the absolute value of correlation with bonding status.

Assessing whether modules are enriched for DE genes

We were interested in assessing whether Mpink^WT genes are significantly enriched for any the DE pre vs post pairing gene set described above. We assessed for enrichment of DE gene sets using a hypergeometric test.

WT module preservation analysis (Zsummary)

We assessed how well modules^WT are preserved in mutant samples. We calculated a Z_summary statistic for each module^{WT 120} using the WGCNA::modulePreservation function with 200 permutations. The Z_summary measure combines module density and intramodular connectivity metrics to a composite statistic where Z>2 suggests moderate preservation and Z>10 suggests high preservation. We found of all the n=15 modules^WT, Mpink^WT had the lowest Z_summary (Z_summary=2) in mutant samples.

Histology

Immunohistochemistry (IHC)

Sexually naïve voles group housed by sex after weaning (P0, P21, adults:10 - 14 weeks old), were perfused with 25mL ice cold phosphate buffer (PBS) and fixed in 25mL of 4% paraformaldehyde in PBS (PFA). Whole brain was dissected and post fixed overnight at 4°C in 4% PFA in the dark. The tissue was cryoprotected for 24 hours in 30% sucrose in PBS and then embedded in OCT over dry ice to prevent cracking. 60μm serial sections from the entire PVN were collected in PBS and processed as free-floating sections.

Sections were blocked using 10% donkey serum in 0.1% Triton X 100 in PBS and then stained with primary antibodies (Oxytocin Polyclonal Antibody, Invitrogen, Rabbit Polyclonal, Cat #PA5-26701, 1:1000 and Anti-Neurophysin 2/NP-AVP Antibody, clone PS 41, Millipore Sigma, Mouse Monoclonal, Cat#MABN845, 1:10,000) in 1% donkey serum with 0.1% Triton X 100 (staining buffer) in PBS at 4°C with shaking. The sections were rinsed and soaked in secondary antibodies (Cy™3 AffiniPure Donkey Anti-Rabbit IgG, Jackson ImmunoResearch Laboratory, Cat #711-165-152, 1:500 and Alexa Flour 488 Donkey Anti Mouse lgG, Thermo Fisher Scientific, Cat #A21202, 1:500) in staining buffer for 2 hours at room temperature. Sections were DAPI stained and transferred to glass slides in sequence and mounted in Aquamount. Slides were imaged using a confocal microscope (Zeiss LSM 700) and the zstacks were quantified (cell count by hand and area hand selected and measured) in ImageJ. Comparisons of cell density (count/area) were made in age and sex matched samples.

In situ hybridization (ISH)

ISH was performed using RNA probes prepared as described previously¹²¹, to detect expression of Calcr, Dlk1, and Dgkb mRNA in the adult vole brain. We prepared RNA sense and anti-sense probes corresponding to 706bp (CalcR), 773bp (Dgkb), and 796bp (Dlk1) for genes of interest identified in the WT vs Oxtr DEX list from our RNAseq data. For ISH on adult animals, 7-9 week old WT and Oxtr^1-/- voles of both sexes were either group housed with same sex conspecifics or paired for 4 days with a WT animal of the opposite sex following the timed mating protocol described above. The brains were processed as described for IHC, then embedded in OCT and sectioned at 50μm, mounted immediately after sectioning onto SuperFrost Plus glass microscope slides and stored at -80°C until ISH staining. Slides were fixed in 4% PFA for 20 minutes, rinsed in PBS, treated with proteinase K (10μg/mL, Roche) for 20 minutes, rinsed in PBS, and fixed again in 4% PFA for 5 minutes at room temperature. Slides were acylated for 10 mins then rinsed in 1% Triton100X in PBS followed by PBS. The slides were equilibrated in warm hybridization solution for 1-2 hours at 65°C and subsequently incubated with a temporary paraffin cover for 14 - 18 hours at 65°C in fresh hybridization buffer containing RNA probe. After incubation, the slides were dipped in 5X SSC warmed to 72°C to remove the parafilm coverslips and then washed in 0.2X SSC warmed to 72°C. The slides were blocked with 10% heat inactivated sheep serum (HISS) solution then stained for 12 - 18 hours at 4°C in buffer containing 10% HISS and alkaline phosphatase-conjugated sheep anti-digoxigenin antibody (1:2000, Roche).

After extensive washing, the slides were incubated for ∼72 hours at 37°C in staining solution containing nitro blue tetrazolium and 5- bromo-4-chloro-3-indolyl-phosphate (Roche). The slides were finally washed, fixed in 4% PFA, and coverslipped. Slides were imaged using a Keyence fluorescence microscope and the z-stacks were quantified (cell count and area) in ImageJ. For signal quantitation, positive puncta were counted in every third coronal section encompassing the entire anterior to posterior gene expression domain. Cell counts were normalized by total area for a gene expression domain in each section. Positive signal was quantified by calculating total area of signal multiplied by the average intensity of signal (total area x mean pixel intensity) in gene expression domains of a section; every third section was quantified and added. Comparisons of cell counts were made in age and sex matched samples.

Quantification and Behavioral Analysis

For detailed statistical methods please see Table S1. The number of animals used was based on previous studies in the field by our group and others, combined with a power analysis. Assumptions of independence and normality were considered with multiple tests. Preference index (PI) was calculated as (partner huddle duration – stranger huddle duration)/(partner huddle duration + stranger huddle duration). Sliding partner preference was calculated by calculating PI for a single animal in 20-minute intervals where each interval overlapped the next by 15 minutes (e.g. 0-20mins, 5-25mins…). To look at convergence in sliding PI the area under the curve (Area above 0 – Area below 0, AUC) for each hour of the assay was determined for each animal and then pooled for each choosing animal by sex and then condition. Preference was determined as significantly convergent if the AUC from the hour of assay for all animals within the group significantly differed from 0 using a One sample t Test. The level of statistical significance for each test was set at p =< 0.05. Outliers were detected using a combination of z-score (+/- 3) and plots of original and log-transformed data. Analyses were completed in RStudio (version 2024.04.2+764) and GraphPad Prism (version 10.3.1)

Acknowledgements

The authors would like to thank members of the Manoli lab for assistance, advice, helpful discussion and comments on the manuscript. J. Tollkuhn and M. Brainard provided valuable insights to the structure, content and design of the manuscript. Additionally, we would like to thank S. Garg for his assistance in data analysis and help designing the analytic pipeline. The authors received funding from National Institutes of Health grant R01MH123513 (D.S.M.), National Science Foundation grant 1556974 (D.S.M.), Burroughs Wellcome Fund 1015667 (D.S.M.), Whitehall Foundation grant 2018-08-83 (D.S.M.), One Mind Foundation A137726 (D.S.M.), National Institutes of Health grant R25MH060482 (K.M.B.), AP Giannini Foundation Fellowship P0534952 (K.M.B.), Larry L. Hillblom Foundation Fellowship 2020-A-023-FEL (K.M.B.), and National Institute of Mental Health grant R01MH123178 (K.S.P.)

Additional information

Author contributions

Conceptualization, R.S., K.M.B., K.L., and D.S.M.; methodology, A.B., A.J.W., D.S.M., K.L.B., K.M.B., K.S.P., and R.S.; software, K.Q., and S.W.; formal analysis, A.E., R.S., and K.M.B.; investigation, R.S., K.M.B., A.E., B.W., G.W., S.W., K.Q., R.D.L., K.L., N.H., B.A.S., M.C.H., M.S., R.K., A.C., D.G., L.C.C., N.L.G.; visualization, R.S., K.M.B., B.W., and A.E.; funding acquisition, D.S.M., K.S.P., and A.J.W; writing– original draft, R.S., and K.M.B.; writing-review & editing, R.S., K.M.B., D.S.M., N.H., K.S.P., A.E., and A.B.

Additional files

Supplementary Figures

Supplementary Figure Legends

Table 1

Table 2