Sexual coordination in a whole-brain map of prairie vole pair bonding

Sexual bonds are central to the social lives of many species, including humans, and monogamous prairie voles have become the predominant model for investigating such attachments. We developed an automated whole-brain mapping pipeline to identify brain circuits underlying pair-bonding behavior. We identified bonding-related c-Fos induction in 68 brain regions clustered in seven major brain-wide neuronal circuits. These circuits include known regulators of bonding, such as the bed nucleus of the stria terminalis, paraventricular hypothalamus, ventral pallidum, and prefrontal cortex. They also include brain regions previously unknown to shape bonding, such as ventromedial hypothalamus, medial preoptic area and the medial amygdala, but that play essential roles in bonding-relevant processes, such as sexual behavior, social reward and territorial aggression. Contrary to some hypotheses, we found that circuits active during mating and bonding were largely sexually monomorphic. Moreover, c-Fos induction across regions was strikingly consistent between members of a pair, with activity best predicted by rates of ejaculation. A novel cluster of regions centered in the amygdala remained coordinated after bonds had formed, suggesting novel substrates for bond maintenance. Our tools and results provide an unprecedented resource for elucidating the networks that translate sexual experience into an enduring bond.


Introduction
Bonds are essential to the social lives of many species, enabling individuals to coordinate 29 parental care, territorial defense, or other shared activities (Emlen & Oring, 1977;Marlowe, 30 2000; Trivers, 1972). Among humans, friendships and a happy marriage protect against a wide  The neurobiology of bonding is most studied in the prairie vole, a socially monogamous 38 rodent in which males and females form bonds, share a nest, and raise young together 39 (Lieberwirth & Wang, 2014; L. J. Young & Wang, 2004). In prairie voles, as in many pair-40 bonding species, bonds form in response to courtship and repeated mating (Carter et al., 1995; Reference Atlas (ARA; Figure 1A, Figure Supplement 1a, and Video 1) (Dong, 2008). This 100 alignment revealed that the prairie-vole brain is ~30% larger than the mouse brain and distinct in 101 shape ( Figure 1A), but that relative volumes of brain regions were consistent across species and 102 showed no evidence of sexual dimorphism (Figure Supplement 1a, Supplemental Files 1 and 2). 103 To validate and refine the ARA anatomical borders in the prairie-vole brain, we 104 performed whole-brain Nissl staining as well as iDISCO immuno-labelling targeting the cell-   To coordinate the timing of mating, all subjects were isolated for 4-5 days, females were 140 brought into estrus, and both sexes were screened for sexual receptivity (see Figure 2A). Subjects 141 assigned to the bonding condition were paired with a novel opposite-sexed individual; to control 142 for non-sexual social affiliation, remaining subjects were re-paired with a familiar same-sex 143 sibling. Members of each pair were isolated on either side of a divider for 2h. Following this 144 acclimation, the divider was removed, and the pair could interact freely. (Figure 2A). Behavior 145 sessions were terminated following a fixed timeline: animals were euthanized just before barrier 146 removal (0h), following initial mating (2.5h), after initial bond formation (6h), or after bond 147 stabilization (22h). Automated behavioral measures, such as proximity, vocalization, and relative 148 movement were scored throughout the behavioral sessions. We performed detailed manual 149 scoring for 1h focal intervals beginning 2h before euthanasia -a timeframe that reveals the

152
Mate and sibling pairs did not differ in locomotor activity or time spent in side-by-side 153 contact ("huddling"; t-tests, FDR-corrected q=0.17 and 0.37; Figure 2D). Mate pairs progressed 154 through known stages of mating behavior, showing elevated rates of anogenital investigation and 155 vocalization (t-tests, q=0.0003 and q=0.0003); males moved more often toward females, while 156 females moved away from their partners (paired t-test, q=0.01), a behavior consistent with 157 female copulatory pacing (Pfaus et al., 2001). Consummatory aspects of mating -mounting, 158 intromission and ejaculation -were common in both the 2.5h and 6h focal windows; by 22h 159 mating was rare, and levels of male-female huddling resembled those of same-sex siblings   Next we measured brain-wide c-Fos immunostaining, a common proxy for neuronal activity and 180 plasticity (Sheng & Greenberg, 1990). To identify brain areas that differed between mate-paired 181 and sibling controls, we used a generalized linear model (GLM) to compare two alternative 182 models. A null model included terms for sex, time, and block; a full model included these ROIs that differed significantly (permutation test, FDR-corrected α < 0.1; Figure 3B), 68 regions 188 were anatomically distinct.

189
A hierarchical cluster analysis assigned these 68 ROIs into 8 groups ( Figure 3C-D, Figure   190 Supplement 1d-e). Each of these groups, or 'clusters,' included multiple brain regions but often 191 centered on a specific structure or group of structures. The purple "BST" cluster, for example,  The orange "Amyg" cluster contained the lateral hypothalamus and various olfactory related 198 regions, but the majority of the cluster consisted of amygdalar nuclei such as cortical and medial 199 amygdala and the anterior amygdala area. The light orange "AUD" cluster contained some 200 regions within auditory cortex, and the red "Thal" cluster involved a variety of thalamic regions.

201
The light purple "AOB" cluster contained minimally correlated regions such as the accessory 202 olfactory bulb, and so we do not consider this group to be one of the major clusters. Comparing     The second canonical correlate (CC2) captured differences between animals who were isolated 236 or paired and loaded particularly highly on the limbic cortical cluster (PFC); we interpret this 237 dimension as capturing appetitive aspects of behaviors ( Figure S2f-h).

238
The first two CC dimensions showed no evidence of sexual dimorphism, while the third 239 dimension suggested subtle but non-significant differences between male and female mates 240 ( Figure S2f-g). A formal two-model comparison isolating the effect of the sex*pairing 241 interaction revealed no ROIs that differed following FDR correction. If we limit the comparison 242 to just those 68 areas identified as responding to pairing, we find 29 unique ROIs that exhibit 243 evidence of sex*pairing interaction, although these effects are much weaker than those identified 244 in our above analysis of pairing (Video 2).

245
Although differences between males and females were modest, the CCA suggested 246 substantial individual differences among mated pairs. We used these data to test whether pairing 247 involved a coordinated change in activities across brain regions. We find that during mating and 248 bond formation (2.5-6h), activity is strongly correlated across ROIs ( Figure 4E-F and Figure   249 S2i). Similar to these patterns distributed across putative pair-bonding regions, CC1 brain scores 250 show a strong correlation between male and female mates (r=0.93, p<0.0001, Figure 4G).

251
Following bonding (22h), within-pair correlations are confined to the coupling of activity 252 between the hypothalamus and amygdala, including between MPOA and medial amygdala. In 253 siblings, who were housed together since birth, hypothalamus-amygdala correlations are 254 relatively uniform across time points, while correlations outside of these regions are largely 255 absent ( Figure S2i).

256
The principal nucleus of the BST (prBST), a subregion of the posterior BST, loads highly 257 on CC1 ( Figure 4B), and exhibits a strong response to mating in the 2.5 and 6h time points 258 (permutation test, q=0.003, Figure 4C); we found that this brain region also shows strong 259 correlations in heterosexual pairs (r=0.851, p<0.0001, Figure 4G). Remarkably, the number of  Here, we present a brain-wide network of IEG+ circuits active as mating experience elicits a 284 pairbond. We find no evidence of anatomical sexual dimorphism, and only modest evidence of 285 dimorphic function. Our whole-brain mapping of pair-bond formation implicates 68 unique brain 286 regions, including 18 regions that are within the primary brain network proposed for prairie-vole 287 bond formation (Walum & Young, 2018). The 68 identified regions are more strongly 288 anatomically connected to one another than predicted by chance, and so can be interpreted as 289 circuits active during pairing.  Ideas and Speculation 343 This brain-wide map of IEG induction provides a uniquely comprehensive perspective on the 344 circuits that enable sexual behavior to become a bond. To help organize this enormous dataset, 345 we looked to existing literature to compare the implicated circuits to those related to sociosexual 346 behaviors in other species (Figure 5A, B). Overlaying the time course of behaviors and circuit 347 activity suggest distinct stages of bonding and attendant neural function. We divide those stages 348 into mating, bonding, and ongoing bond maintenance.

349
The first canonical correlate captured the major dimensions of both brain and behavioral 350 variation, and showed the largest group differences surrounding the 2.5 hr timepoint.  suggesting that repeated copulation is a means of partner assessment (Dewsbury, 1988; 392 Eberhard, 1996); in prairie voles, we suggest that both males and females are assessing the 393 ability of a male to monopolize a female, a trait that would predict male paternity and ability to  After bond formation and stabilization, we find that that mated pairs are similar to co-411 housed siblings in terms of brain-wide IEG induction patterns ( Figure 5C). More remarkable, is Overall, our data survey the brain to identify circuits active as mating produces a pairbond. Brains were registered to a standardized reference brain as previously described. Initial 496 3D affine transformation was calculated using 6 resolution levels followed by a 3D B-spline was not stained for c-Fos due to major issues with the perfusion. The latter two samples were 554 identified as outliers by a Rosner's test (EnvStats R package (Millard, 2013)); their whole-brain 555 c-Fos counts were higher than the rest of the samples (R = 4.61 and 5.504, P < 0.05).  associations between behavioral measures and ROI c-Fos+ cell counts ( Figure S2e,h,i).

642
To identify brain regions that are sensitive to pairing status, we used a model comparison  This method excluded larger brain regions in lieu of smaller and more localized ROIs.

27
Significant "chosen" ROIs were assigned to groups by using hierarchical clustering with the 674 ward D2 method with a Euclidean distance matrix extracted from c-Fos cell counts (Murtagh & 675 Legendre, 2014; R Development Core Team, 2016). The resulting tree was cut at such that ROIs 676 were grouped into anatomically similar clusters (Figure 3, Figure Supplement 1f,g). We used 677 multi-dimensional scaling (MDS) to further interpret the degree of similarity between chosen 678 ROIs and their clustering groups based on the Euclidean distance matrix (Figure Supplement 1e).

679
The MDS method is a form of non-linear dimensionality reduction to visualize similarity in 680 Cartesian space (Cox & Cox, 2008).

681
To confirm whether chosen ROI clusters represented structural circuits, we compared 682 cluster assignments to published data on structural connectivity in the mouse brain (Knox et al., cases, we used data from the next inclusive ROI that was available (e.g., BST data used for 687 BSTpr). Then, for each cluster, we found the mean normalized connection density between the 688 regions. We excluded data from the matrix diagonals (e.g., BST to BST) to emphasize 689 connections between, rather than within, regions. We took an average of these cluster densities 690 values to capture the overall connection density based on our cluster assignments. A permutation 691 test was used to assess whether this connection density was higher than expected by chance. The 692 rows of the connectivity matrix (rows = origin ROIs, columns = target ROIs) were randomly 693 shuffled prior to computing the average cluster density, and this was done 10,000 times to 694 construct a null distribution. This null distribution was then compared to the observed density to 695 estimate its probability. 696 We used canonical-correlation analysis (CCA) was used to investigate the relationships Atlas.git). Raw behavioral and neural data will be made available upon reasonable request.