Brain functional networks associated with social bonding in monogamous voles

  1. M Fernanda López-Gutiérrez
  2. Zeus Gracia-Tabuenca
  3. Juan J Ortiz
  4. Francisco J Camacho
  5. Larry J Young
  6. Raúl G Paredes
  7. Néstor F Díaz  Is a corresponding author
  8. Wendy Portillo  Is a corresponding author
  9. Sarael Alcauter  Is a corresponding author
  1. Instituto de Neurobiología, Universidad Nacional Autónoma de México, Mexico
  2. Silvio O Conte Center for Oxytocin and Social Cognition, Center for Translational Social Neuroscience, Yerkes National Primate Research Center, Department of Psychiatry and Behavioral Sciences, Emory University, United States
  3. Escuela Nacional de Estudios Superiores, Unidad Juriquilla, Universidad Nacional Autónoma de México, Mexico
  4. Instituto Nacional de Perinatología Isidro Espinosa de los Reyes, Mexico

Abstract

Previous studies have related pair-bonding in Microtus ochrogaster, the prairie vole, with plastic changes in several brain regions. However, the interactions between these socially relevant regions have yet to be described. In this study, we used resting-state magnetic resonance imaging to explore bonding behaviors and functional connectivity of brain regions previously associated with pair-bonding. Thirty-two male and female prairie voles were scanned at baseline, 24 hr, and 2 weeks after the onset of cohabitation. By using network-based statistics, we identified that the functional connectivity of a corticostriatal network predicted the onset of affiliative behavior, while another predicted the amount of social interaction during a partner preference test. Furthermore, a network with significant changes in time was revealed, also showing associations with the level of partner preference. Overall, our findings revealed the association between network-level functional connectivity changes and social bonding.

Introduction

The prairie vole (Microtus ochrogaster) is a rodent native of North America whose natural behavior involves pair-bonding, which can be defined as a long-lasting, strong social relationship between individuals in a breeding pair in monogamous species (Walum and Young, 2018). Pair-bonded voles will usually display selective aggression towards unfamiliar conspecifics; biparental care, including paternal behavior and alloparenting; incest avoidance; and reproductive suppression of adult individuals within a family group (Carter et al., 1995). These behaviors make the prairie vole a valuable model to investigate behaviors associated with a socially monogamous reproductive strategy (Young and Wang, 2004), social bond disruption, social isolation, and social buffering (Lieberwirth and Wang, 2016). Being comparable to human-like social interactions, studying the neurobiology of social behavior in the prairie vole may allow further understanding of human social bonding, its alterations in psychological disorders, and overall impact on health (Kiecolt-Glaser et al., 2010).

Pair-bonding-related behaviors in the prairie vole depend on hormonal mechanisms and the activation of emotional, reward, and sensory brain circuits (Walum and Young, 2018), integrating large functional networks (Johnson and Young, 2017). Among these, the mesolimbic reward system and the social decision-making network (SDMN) (Johnson et al., 2017) are the proposed networks to be involved in pair-bonding, being modulated by steroid hormones, dopamine (DA), oxytocin (OXT), arginine vasopressin (AVP), γ-aminobutyric acid, glutamate, and corticotropin-releasing factor, among others (Walum and Young, 2018). A current hypothesis suggests that pair-bonding consists of two different plastic processes: the formation of a neural representation of the partner, and a selective preference for the partner, that is the maintenance of the pair-bond (Young and Wang, 2004; Walum and Young, 2018). In this process, an association has to be created between the reinforcing properties of sex (mating) and the olfactory signature from the partner (Ulloa et al., 2018; Young and Wang, 2004). In broad terms, both OXT and AVP are necessary and sufficient for the formation of the pair-bond (Young and Wang, 2004; Lieberwirth and Wang, 2016), and their release during social and sexual interactions are the likely triggers and critical modulators of the mentioned network, since most regions included in the SDMN express their corresponding receptor binding (Johnson and Young, 2017). DA would also be released in concert with OXT and AVP in specific regions to modulate the adequate display of behavior and formation of the pair-bond (Young and Wang, 2004).

While prairie voles of both sexes already display changes in behavior by 24 hr of cohabitation with mating as a result of pair-bonding (Wang and Aragona, 2004; Williams et al., 1992), the nucleus accumbens (NAcc), for example, has substantially more D1-like receptor binding 2 weeks after female exposure, relative to non-pair-bonded males (Aragona et al., 2006), a phenomenon that is also been observed in female voles (Resendez et al., 2016). This evidence suggests that long-term plasticity is relevant for pair-bond maintenance. Furthermore, male prairie voles that were pair-bonded for 2 weeks displayed selective aggression toward conspecific male and female strangers (Gobrogge et al., 2007), and the 2 week time frame has also been found relevant in bond disruption (Insel and Hulihan, 1995) and partner loss (Tabbaa et al., 2017). While there is a possibility that these changes in behavior involve the interplay of broad networks, it has not been directly explored in vivo.

Recently, novel electrophysiologic and optogenetic techniques have been employed to demonstrate that the functional connectivity between the NAcc and the medial prefrontal cortex (mPFC) during initial cohabitation in female prairie voles modulates the affiliative behavior with their potential partner (Amadei et al., 2017), providing exciting data of the relevance of such corticostriatal interactions for social bonding. However, this approach does not allow the study of the interaction of multiple brain regions, that is networks, and their relevance of such interactions in the process of pair-bonding. Neuroimaging methods may provide the alternative to explore such networks in a longitudinal fashion, since few studies have made use of positron-emission tomography to explore limited aspects of such longitudinal changes (Bales et al., 2007), providing the first longitudinal evidence of neurophysiological changes associated with pair-bonding. Potentially, functional magnetic resonance imaging (fMRI) may be the ideal tool to explore the longitudinal changes in functional brain networks (Damoiseaux et al., 2006), providing high spatial resolution, wide brain coverage, and being minimally invasive to explore longitudinal changes. In particular, resting-state functional magnetic resonance imaging (rsfMRI) explores the low-frequency fluctuations (<0.1 Hz) of the blood-oxygen-level-dependent (BOLD) signal, which has proven to remain highly synchronized within the sensory, motor, and associative networks in the brain of humans (Damoiseaux et al., 2006), non-human primates (Rilling et al., 2007), and rodents (Gozzi and Schwarz, 2016), including the prairie vole (Ortiz et al., 2018), turning this technique into a promising tool for translational research. Indeed, functional connectivity explored by means of the correlation of rsfMRI signals of anatomically separated brain regions (Friston et al., 1993) has demonstrated to be correlated with neuronal activity (Mateo et al., 2017), and it is subject to change in plastic processes such as learning and memory in both humans (Jolles et al., 2013) and rodents (Nasrallah et al., 2016). A recent study found differences in diffusion-weighted imaging and resting-state functional connectivity between two distinct populations of male prairie voles (Ortiz et al., 2020). Here, we make use of this non-invasive technique to explore behavioral correlates with brain functional connectivity and its longitudinal changes associated with pair-bonding in male and female prairie voles.

Results

Baseline functional connectivity predicts the display of affiliative behavior

Thirty-two 3-month-old sexually naïve female (N = 16) and male (N = 16) prairie voles (M. ochrogaster) were used in the study. Prairie voles underwent three magnetic resonance imaging (MRI) acquisition sessions: a baseline scan before cohabitation, a second scan 24 hr after the onset of cohabitation, and a third scan 2 weeks after the onset of cohabitation (Figure 1A). The final imaging sample consisted of 90 datasets, with only six subjects missing one session (see Materials and methods). The day after the baseline scanning session, female and male voles unrelated to each other were randomly assigned as couples and placed for cohabitation in a new home to promote ad libitum mating and social interaction. Four days before cohabitation, silastic capsules (Dow Corning Silastic Laboratory Tubing; ThermoFisher Scientific, Pittsburg, PA) containing estradiol benzoate (E2B; Sigma–Aldrich, St. Louis, MO) were implanted in previously ovariectomized female voles to enable sexual receptivity and promote mating (see Materials and methods). The first 6 hr of cohabitation was video recorded for analysis of social and mating behavior. Mount (M ± SEM: 65.4 ± 31.7 min), intromission (116 ± 35.3 min), and ejaculation (125 ± 34.4 min) latencies were obtained for male voles (N = 16). Lordosis latency (22.3 ± 13.3 min) was also measured on females (N = 16), and huddling latencies (69.5 ± 15.8 min) were obtained for each male and female pair. Three of the 16 couples did not mate during the recorded period, but all voles displayed huddling and licking/grooming behavior with their sexual partner. Once joined, voles remained housed in couples for the rest of the experiment and were only separated for MRI scanning sessions and behavioral tests.

Figure 1 with 5 supplements see all
Experimental design and brain regions of interest.

(A) Sequence of experiments during a 30 day period: Female voles were bilaterally ovariectomized before MR and behavioral protocols. After being allowed to recover from surgery for 10 days, silastic capsules containing E2B (estradiol benzoate) were implanted via s.c. 4 days before cohabitation for sexual receptivity induction. Once couples went under cohabitation, they were housed together for the rest of the experiment and were only separated for PPT and MR scanning sessions. OVX: ovariectomy surgery. MR: magnetic resonance imaging scanning session. PPT: partner preference test. (B) Regions of interest (ROIs) for network functional connectivity analyses. Antero-posterior coronal slices of the prairie vole template overlayed with ROI masks with the resolution used in the analysis. Each color represents a different ROI. ACC: anterior cingulate cortex. AON: anterior olfactory nucleus. BLA: basolateral amygdala. BNST: bed nucleus of the stria terminalis. DG: dentate gyrus. dHIP: dorsal hippocampus. MeA: medial amygdala. MOB: main olfactory bulb. LS: lateral septum. mPFC: medial prefrontal cortex. NAcc: nucleus accumbens. PVN: paraventricular nucleus. RSC: retrosplenial cortex; VP: ventral pallidum. vHIP: ventral hippocampus. VTA: ventral tegmental area. (C) 3D views of ROI masks embedded within the prairie vole template.

To explore relationships between socio-sexual behavior and functional connectivity, 16 regions of interest (ROIs) were defined according to their previously reported relevance in the process of pair-bond formation and maintenance (Johnson and Young, 2017; Lieberwirth and Wang, 2016; Walum and Young, 2018), which were the following: anterior cingulate cortex (ACC), anterior olfactory nucleus (AON), basolateral amygdala (BLA), bed nucleus of the stria terminalis (BNST), lateral septum (LS), medial amygdala (MeA), main olfactory bulb (MOB), medial prefrontal cortex (mPFC), nucleus accumbens (NAcc), retrosplenial cortex (RSC), paraventricular nucleus of the hypothalamus (PVN), ventral pallidum (VP), ventral tegmental area (VTA), dentate gyrus (DG), dorsal hippocampus (dHIP), and ventral hippocampus (vHIP) (Figure 1B,C). For each subject and session, average time series for each region were extracted from the pre-processed fMRI datasets (Figure 1—figure supplement 1). The latter were used to obtain connectivity matrices, based on the partial-correlation estimates for all possible pairs of ROIs (Figure 1—figure supplement 2).

The amount and proneness to affiliative behavior has shown reliability in describing the level of sociability and mating systems in Microtus voles. Through the measurement of huddling, Salo et al., 1993 were able to distinguish social behavior in four different species of voles by scoring huddling when the pair was either sitting or lying in bodily contact, with the prairie vole having the highest huddling accumulation. Other studies in voles have also measured huddling latency to assess levels of affiliative behavior, since it may influence pair-bond induction and formation (Shapiro and Dewsbury, 1990; Amadei et al., 2017).

In order to identify potential relationships between behavior and functional connectivity with a network perspective, huddling latencies during cohabitation in voles of both sexes (N = 28) were tested as linear covariates of baseline functional connectivity data, that is before cohabitation, with Network-Based Statistic (NBS) framework, using the NBS Toolbox (Zalesky et al., 2010).

NBS analysis found significant negative linear relationship (p≤0.001) with a large network including VP, MeA, LS, VTA, RSC, BLA, NAcc, ACC, and the DG (Figure 2A). Our results show that the higher the connectivity between these regions before cohabitation, the shorter the huddling latencies during cohabitation in voles of both sexes. Additionally, a posteriori Pearson’s correlations confirmed the correlation strength between each connection and huddling latencies: VP–MeA (r(26) = −0.468, p=0.011), MeA–LS (r(26) = −0.372, p=0.051), LS–VTA (r(26) = −0.502, p=0.006), VTA–RSC (r(26) = −0.586, p=0.001), RSC–PVN (r(26) = −0.382, p=0.044), RSC–BLA (r(26) = −0.362, p=0.058), BLA-DG (r(26) = −0.420, p=0.025), BLA–NAcc (r(26) = −0.543, p=0.002), and NAcc–ACC (r(26) = −0.378, p=0.047), (Figure 2B–J). These results show that functional connectivity between these regions reflect the predisposition to display affiliative behavior in both female and male prairie voles.

Relationships between baseline functional connectivity and affiliative behavior (huddling) during cohabitation with mating in male and female prairie voles.

(A) Representation of a prairie vole brain with regions (nodes) that constitute the network with a significant negative association with huddling latency. Scatter-plot graphs (B–J) of the connections in a with best line fit between baseline functional connectivity (Fisher z-transformed partial-correlation values) and huddling latencies (minutes) during cohabitation. The higher the connectivity between these regions before cohabitation, the shorter the huddling latencies during cohabitation in voles of both sexes. ACC: anterior cingulate cortex. BLA: basolateral amygdala. DG: dentate gyrus. LS: lateral septum. MeA: medial amygdala. NAcc: nucleus accumbens. PVN: paraventricular nucleus. RSC: retrosplenial cortex. VP: ventral pallidum. VTA: ventral tegmental area.

Figure 2—source data 1

Functional connectivity values of edges that correlate significantly with huddling latencies in prairie voles.

https://cdn.elifesciences.org/articles/55081/elife-55081-fig2-data1-v2.xlsx

Prairie voles of both sexes show partner preference after cohabitation

Between 48 and 72 hr of cohabitation, partner preference was evaluated on each subject (N = 32) to assess pair-bonding behavior. This protocol was based on a previously described test (Williams et al., 1992Figure 3B, see Materials and methods). The partner preference index revealed a significant difference between the proportion of time spent on the incentive area related to the partner (median = 0.63) with the area related to the stranger vole (median = 0.37) for all subjects (U = 378, p=0.0365, effect size r = 0.32) (Figure 3B). No significant differences were found between males and females in their preference for the partner (U = 121, p=0.81, effect size r = 0.05) or the stranger voles (U = 118, p=0.72, effect size r = 0.07), and there were also no significant differences in partner preference between the time periods when PPT tests were performed (48 and 72 hr) (U = 119, p=0.75, effect size r = 0.06) (Figure 3—figure supplement 1).

Figure 3 with 1 supplement see all
Partner preference test.

(A) Representative figure showing the design of the arena in which voles were tested for partner preference test. (B) Between 48 and 72 hr of cohabitation, partner preference was evaluated on each subject (N = 32). Partner preference index revealed a significant difference between the time spent on the incentive area related to the partner, with the incentive area related to the stranger vole. Boxplot graphs show whiskers with 10–90 percentiles; horizontal line inside the box shows data median, and ‘+’ represents data mean. (*) denotes significance at p<0.05.

Figure 3—source data 1

Values of parter preference index obtained from the proportion of time spent on social incentive areas related to the partner and stranger voles.

https://cdn.elifesciences.org/articles/55081/elife-55081-fig3-data1-v2.csv

Functional connectivity at 24 hr of cohabitation predicted social interaction during the partner preference test

During the partner preference test (PPT), the proportion of time spent on social incentive areas over the total recorded time was calculated for analysis in all subjects (see Materials and methods). The percentage of time spent on social incentive areas over the total time of the test had a mean of 79.70 ± 2.65%. Spearman correlation analyses showed that the amount of time spent on social incentive areas had no relationship with the partner preference index (rs(30) = 0.08912 p=0.6277) (Figure 4—figure supplement 1). However, we tested if functional connectivity 24 hr after the start of cohabitation had a relationship with the total time spent on social incentive areas during the test, regardless if the subject interacted with the partner or stranger stimulus vole (N = 32). Through NBS analysis, a significant positive linear relationship (p=0.013) was found with the following network: LS–NAcc–mPFC–MeA–VP–MOB–DG (Figure 4A). These findings suggest that the lower the connectivity between these regions 24 hr after the onset of cohabitation, the longer the subject would interact with a conspecific of the opposite sex during the PPT, which was evaluated between 48 and 72 hr after the onset of cohabitation. The correlation strength between each connection of such network and the percentage of time on social incentive areas was obtained with a posteriori Pearson’s correlations: LS–NAcc (r(30) = −0.369, p=0.037), NAcc–mPFC (r(30) = −0.312, p=0.081), mPFC–MeA (r(30) = −0.373, p=0.034), MeA–VP (r(30) = −0.376, p=0.033), VP–MOB (r(30) = −0.531, p=0.001), and MOB–DG (r(30) = −0.404, p=0.021) (Figure 4B–G).

Figure 4 with 1 supplement see all
Relationships between functional connectivity 24 hr after the onset of cohabitation and social interaction during partner preference test in male and female prairie voles.

(A) Representation of a prairie vole brain with regions (nodes) that constitute the network with a significant negative association with the amount of social interaction during the PPT. Scatter-plot graphs (B–G) of the connections in a with best line fit between baseline functional connectivity (Fisher z-transformed partial-correlation values) and time on social incentive areas during cohabitation (percentage). The lower the connectivity between these regions at 24 hr of cohabitation, the longer the time spent on social incentive areas during the PPT. DG: dentate gyrus. LS: lateral septum. MeA: medial amygdala. MOB: main olfactory bulb. mPFC: medial prefrontal cortex. NAcc: nucleus accumbens. VP: ventral pallidum.

Figure 4—source data 1

Values of percentage of time spent on social incentive areas related to the partner and stranger voles.

https://cdn.elifesciences.org/articles/55081/elife-55081-fig4-data1-v2.xlsx

Male and female voles share network-level changes related to social bonding in different time points

To determine whether cohabitation with mating induced changes in brain functional connectivity between sessions, longitudinal data was analyzed with linear mixed models (LMM), implemented via the Network-Based R-statistics package (NBR; Gracia-Tabuenca and Alcauter, 2020), which allows to implement LMM in an NBS framework (see Materials and methods). NBR analysis yielded significant session effects, that is changes between baseline, 24 hr, and 2 weeks of cohabitation (pFWE = 0.0482), in a network consisting of 10 regions: ACC, BLA, dHIP, LS, mPFC, NAcc, RSC, vHIP, VP, and VTA. No significant differences were found between male and female voles, a result that led us to run the analysis without the sex variable. The same network component was obtained (pFWE = 0.042) (Figure 5A).

Figure 5 with 2 supplements see all
Network-wide changes in time in brain functional connectivity of female and male prairie voles.

(A) NBR analysis via Linear mixed models (LMM) analysis results represented in a prairie vole brain with regions (nodes) comprising the brain network that undergoes significant changes in functional connectivity (Fisher z-transformed correlation values) after cohabitation with mating. Interregional connectivity (edges) is shown by color code. Red: increase of functional connectivity; blue: decrease of functional connectivity. ACC: anterior cingulate cortex. BLA: basolateral amygdala. dHIP: dorsal hippocampus. LS: lateral septum. mPFC: medial prefrontal cortex. NAcc: nucleus accumbens. RSC: retrosplenial cortex. vHIP: ventral hippocampus. VP: ventral pallidum. VTA: ventral tegmental area. (B–K) Functional connectivity values in violin plots showing full distribution of data and median. Connecting lines track longitudinal data of each subject between regions through specific MR acquisition time points (Session): Baseline, 24 hr, and 2 weeks of cohabitation. Color codes for data points and connecting lines distinguish male (cyan) from female subjects (pink). False discovery rate (FDR) post hoc significant differences are shown: *<0.05, **<0.01, ***<0.001.

Figure 5—source data 1

Functional connectivity values that changed significantly across sessions (Ses; 1, 2, and 3) in male and female prairie voles (females: F; males: M).

NA: data not available.

https://cdn.elifesciences.org/articles/55081/elife-55081-fig5-data1-v2.csv
Figure 5—source data 2

Estimated marginal means (EMM) post hoc comparisons values obtained from the longitudinal analysis of functional connectivity data (edge), via linear mixed models with NBR.

Values are corrected for a false discovery rate (FDR) of q < 0.05 (PH_pFDR).

https://cdn.elifesciences.org/articles/55081/elife-55081-fig5-data2-v2.csv

Post hoc analyses (false discovery rate [FDR] corrected) identified differential longitudinal changes among rsfMRI sessions. Specifically, four connections had significant changes 24 hr after the onset of cohabitation: ACC–vHIP (t(56) = −2.766; p=0.0077), LS–RSC (t(56) = −3.270; p=0.0018), and LS–VP (t(56) = −2.540; p=0.0138) showed increased functional connectivity, while LS–dHIP (t(56) = 3.004; p=0.0039) had decreased connectivity in the same period of time. The edges LS–VP, LS–dHIP, and LS–RSC had no differences between baseline (before cohabitation) and the third session (2 weeks after the onset of cohabitation), suggesting acute plastic changes related to 24 hr of cohabitation; however, LS–mPFC (t(56) = −2.856; p=0.0060) and VP–VTA (t(56) = −2.752; p=0.0079) have increased connectivity, while ACC–LS (t(56) = 3.227; p=0.0020) and NAcc–VTA (t(56) = 2.970; p=0.0043) had decreased connectivity between the first and third scanning sessions (baseline vs 2 weeks of cohabitation). In addition, four connections reflected differences between 24 hr and 2 weeks after the onset of cohabitation (sessions 2 and 3): ACC–LS (t(56) = 3.675; p=0.0005) and BLA–NAcc (t(56) = 2.593; p=0.0121) had decreased connectivity, but LS–mPFC (t(56) = −4.013; p=0.0001) and mPFC–dHIP (t(56) = −2.605; p=0.0117) exhibit increased connectivity after 2 weeks of cohabitation (third session), suggesting long-term functional changes related to cohabitation (Figure 5B–K). Overall, the node with the most changes in functional connectivity was the LS.

To control for potential confounders, we selected an independent set of brain structures not expected to be involved in social behavior and were tested with the same data and analysis methods than the ROIs mentioned previously. These nine regions were the primary auditory area (AUDp), the cerebellar cortex (CBX), forceps minor of the corpus callosum (fmi), laterodorsal thalamic nucleus (LD), primary motor area (MOp), motor-related medulla (MY), supplemental somatosensory area (SSs), primary visual area (VISp), and ventricle areas (Vent) (Figure 5—figure supplement 1). We used the same data and analysis methods to test an independent set of brain structures not expected to be involved in social behavior. Linear mixed models (LMM) implemented via the Network-Based R-statistics package (NBR; Gracia-Tabuenca and Alcauter, 2020), considering sex as a fixed variable, and session and intercept as random variables, found no significant networks with differences between sessions (pFWE = 0.1568), or sex*session interactions (pFWE = 0.619 and pFWE = 0.3554), suggesting that in both female and male prairie voles, cohabitation with mating and social bonding do not influence changes in functional connectivity in these regions between any session (baseline, 24 hr, and 2 weeks after the onset of cohabitation).

Group-independent component analysis (gICA) was also performed to address the exploration of large-scale functional brain networks and potential differences between sessions. The gICA revealed five components associated with sensory and motor cortices, putative default-mode, and salience networks, a striatum-centered component, and two components with relevant connectivity of the ventral hippocampi, with a degree of lateralization found in some of the latter (Figure 5—figure supplement 2). These results are strikingly similar to other networks reported previously in the male prairie vole (Ortiz et al., 2018), although this is a larger sample also including female voles and an optimized anesthesia protocol for the detection of rsfMRI networks in rodents (Grandjean et al., 2020). In order to evaluate whether networks found through gICA maps had significant changes between sessions, dual regression was applied and group differences across sessions were assessed using two-sample paired t-tests (see Materials and methods). However, the obtained FWE-corrected p statistics found no significant differences between sessions. Being a voxel-wise method, it may not be sensitive enough to detect punctual, region-specific changes in brain functional connectivity.

Specific network connections that undergo changes during cohabitation with mating correlate with partner preference in male and female voles

Although several regions were found to change significantly after cohabitation and social bonding, we tested if any specific connection of the detected network component (Figure 5A) could have a relationship with the partner preference index obtained between 48 and 72 hr after the onset of cohabitation, which was used to evaluate the level of pair-bonding in each subject. Two-tailed Pearson’s correlation tests revealed a significant positive relationship in LS–VP baseline functional connectivity (r(26) = 0.435, p=0.020) with the partner preference index, suggesting that the higher the baseline functional connectivity between these regions, the higher the partner preference index would be (Figure 6A). Also, significant negative relationships were found between the partner preference index and BLA–NAcc functional connectivity at 24 hr of cohabitation (session 2) (r(30) = −0.464, p=0.0074), and in LS–RSC functional connectivity at 2 weeks of cohabitation (session 3) (r(28) = −0.437, p=0.015), which may suggest that the lower the connectivity between these regions at the specified time points, the higher the partner preference index (Figure 6B,C).

Relationships between relevant network connections and partner preference index during cohabitation with mating in male and female prairie voles.

Scatterplots (A–C) that show significant correlations (with best line fit) between functional connectivity (Fisher z-transformed partial-correlation values) and partner preference index in network connections that undergo longitudinal changes. BLA: basolateral amygdala. LS: lateral septum. NAcc: nucleus accumbens. RSC: retrosplenial cortex. VP: ventral pallidum.

Figure 6—source data 1

Functional connectivity values of edges that correlate significantly with partner preference index of prairie voles in each session: session 1 (s01), session 2 (s02), and session 3 (s03), which correspond to baseline, 24 hr, and 2 weeks after the onset of cohabitation, respectively.

https://cdn.elifesciences.org/articles/55081/elife-55081-fig6-data1-v2.xlsx

Discussion

Several studies have described the relevance of different brain regions involved in pair-bond induction and maintenance in prairie voles (Johnson and Young, 2017; Walum and Young, 2018). However, longitudinal explorations of the brain before and after pair-bonding are scarce (Bales et al., 2007), especially from a network perspective. Here, by using rsfMRI, we were able to detect a brain network in which baseline functional connectivity (before cohabitation) predicted the latency for huddling behavior during the first hours of cohabitation, providing the potential neurofunctional substrate for the variability in affiliative behavior and further extending the recent findings that the corticostriatal electric activity modulates social bonding in prairie voles (Amadei et al., 2017). A relationship between functional connectivity and social interaction was also found, in which a network detected from data obtained 24 hr after the onset of cohabitation predicted the amount of social interaction during the PPT. Finally, our results reflected significant longitudinal changes in functional connectivity of prairie voles after pair-bonding. Post hoc analyses revealed differential short- and long-term connectivity changes mainly involving the lateral septum (LS), with three network connections correlating with the level of partner preference in different sessions. We further discuss the potential neurophysiological basis and implications of our findings.

Correlations between network functional connectivity and social behavior

Although each of the identified networks would likely act as a whole or possess emergent properties, in the following text, we mention previous evidence in rodents that may aid in their functional interpretation, and each component will be dissected in segments to better understand node relationships and their putative role in social bonding and related behavior.

Even though it has been reported that in prairie voles, 24 hr of cohabitation or 6 hr of ad libitum mating is sufficient for a pair-bond to be developed (Williams et al., 1992), a considerable amount of evidence has shown that other factors influence its development and maintenance. Specifically, AVP (Ophir et al., 2008) and OXT receptor gene expression and density (King et al., 2016), paternal nurturing (Ahern and Young, 2009), and neonatal isolation (Barrett et al., 2015) have shown to produce variability in the exhibition of prairie vole social behavior. Though they were under the same experimental conditions, subjects from both sexes in this study showed a wide behavioral variability not only during their sexual encounters in the first hours of cohabitation, but also in their bonding behavior evaluated between 48 and 72 hr after the onset of cohabitation. It is likely that the sum of previously mentioned factors gives each subject a distinctive brain network configuration that ultimately relates to bonding behavior. Hence, we hypothesized that there may be individual differences in functional connectivity that could explain the variability in socio-sexual behavior.

Correlation with huddling latencies

Indeed, we identified a network for which the functional connectivity at baseline was negatively related to huddling latencies during the first hours of cohabitation. In other words, baseline functional connectivity predicted how quickly subjects would begin affiliative huddling with an opposite-sex conspecific. Huddling is a measurable affiliative behavior in prairie voles and a useful indicator of social receptiveness (Salo et al., 1993). The network included the following connections: VP–MeA–LS–VTA–RSC–PVN, PVN–RSC–BLA–NAcc–ACC–NAcc, and BLA–DG (Figure 3A). This network has regions reported to be involved in social salience, social memory and recognition, spatial memory, and reward-seeking mechanisms. Therefore, it is possible that subjects that exhibited huddling at an earlier time found it more rewarding than those who did not.

The VP plays a major role in reward and motivation (Smith et al., 2009). Since MeA activity in rodents is necessary for social recognition (Ferguson et al., 2001) and responds to sex-specific chemosensory cues (Wang and Young, 1997; Yao et al., 2017), increased MeA-VP connectivity may improve the rewarding response of sex-related social interaction. Additionally, DA neurons are reported to modulate VP responses evoked by amygdala stimulation (Maslowski-Cobuzzi and Napier, 1994). The LS is a structure involved in social recognition and retrieval of relevant social information (Bielsky et al., 2005), in addition that it is known to be anatomically interconnected with regions involved in social behavior, including the hippocampus and the MeA (Risold and Swanson, 1997). In rats, it has been shown that OXT release in the LS is required for the maintenance of social memory and may also be modulated by the MeA according to the relevance of the social stimulus (Lukas et al., 2013). It is likely that MeA-LS connectivity further integrates social information. LS–VTA–RSC functional connections, through the VTA, may contribute to linking contextual information with the midbrain DA system and regulate motivational behaviors (Luo et al., 2011; Maeda and Mogenson, 1981), since the RSC is known to be involved in the processing of spatial and contextual memory in rodents (Todd and Bucci, 2015).

Connectivity between PVN–RSC–BLA may be related to an association of a spatial context with a socially salient stimuli. While OXT projections from PVN may probably influence the regulation of social salience in the whole network (Johnson et al., 2017), it is also reported to be sensitive to external stimuli (Anacker et al., 2014; Liu et al., 2001) that can impact OXT synthesis and release (Smith and Wang, 2014). The BLA may act as an associative site for stimulus-outcome representations (Cardinal et al., 2002), which would allow an appropriate response according to previous social encounters, information that possibly requires hippocampus-associated memory (BLA-DG) (Frey et al., 2001; Tashiro et al., 2007).

In a previous report, the BLA and the ACC were found relevant in coordinating brain activity when social interaction is initiated and in the formation of social recognition memory through gene expression (Tanimizu et al., 2017). The NAcc is known to translate reward-predictive information from the amygdala (BLA–NAcc) to promote cue-evoked, reward-seeking behavioral responses in rodents (Ambroggi et al., 2008), while the ACC has been considered an important region in the decision-making process between sensory perception, motivation, and final motor performance (Assadi et al., 2009). Thus, input from the ACC (NAcc–ACC) may be necessary for a social decision-making process. It is important to note that a recent study demonstrated this particular circuit may be unique in prairie voles (Horie et al., 2020), and in female prairie voles, the functional connectivity of the prefrontal cortex and NAcc after the first encounter predicts affiliative huddling toward a partner, and the activation of such circuit biases later preference toward a partner (Amadei et al., 2017). Our results extend such findings, showing that a larger network including similar corticostriatal connectivity, measured even before the exposure to a potential partner, predicts affiliative behavior. Moreover, such relation is consistent for both males and females, and the circuit includes the amygdala and hippocampus, as predicted by Amadei et al., 2017. Overall, the functional connectivity of this network may indicate the predisposition of a prairie vole to engage in prosocial behavior, which might be influenced by previous experience and other factors mentioned beforehand.

Correlation with amount of social interaction

A network component detected at 24 hr after the onset of cohabitation, LS–NAcc–mPFC–MeA–VP–MOB–DG, was negatively correlated with the amount of social incentive a subject would have during the partner preference test, a test in which it interacted with two conspecifics of the opposite sex, one being the partner and the other a stranger vole. Interestingly, this component includes five regions from the huddling network previously described, and all but two nodes (MeA–VP) are related differently when comparing it with the former network.

Regardless of the level of preference for the partner, prairie voles of both sexes had a differential interest in engaging socially. While the huddling network is related to pre-pair-bonding affiliative behavior, this network is necessarily related to post-bonding behavior, since it was only detected 24 hr after the onset of cohabitation. It is possible that functional connectivity of this component is associated with the modulation of a behavioral response, that is approach or avoidance, according to social olfactory cues that may be rewarding according to previous experience. Besides its social recognition role, in rodents the LS is reported to modulate reward response by decreasing neuronal activity in the NAcc through AVP release (Gárate-Pérez et al., 2021). The NAcc–mPFC pathway has been widely characterized as modulating reward-seeking, goal-oriented behavior (Gill et al., 2010), and the mPFC and the amygdala are extensively interconnected and tune the expression of fear and anxiety (Liu et al., 2020; Marek et al., 2013). The involvement of the MeA–VP hints modulation of sex-related social interaction, particularly via chemosensory stimuli (MOB-DG) (Castro et al., 2020; Liu et al., 2014).

Our data suggests that functional connectivity 24 hr after the onset of cohabitation has a relationship to the level of sociability, reflecting a phenotype independent of pre-bonding behavior and partner preference (Figure 4—figure supplement 1). Further investigation of this network would be necessary to understand the impact of sociability in other behaviors, such as mate-territory guarding and parental nurturing. Apparently, its change in modulation is on a short term, since a correlation at 2 weeks of cohabitation was not found and other regions or networks may play this role on a long-term basis.

Longitudinal changes in functional connectivity after pair-bonding

It has been proposed that pair-bonding results from the convergence of the mesolimbic DA reward circuit and social discrimination circuits (Walum and Young, 2018). The results here presented are consistent with this model, by demonstrating both short- and long-term changes in a brain network including regions largely associated with reward/motivation (ACC, VTA, NAcc, VP) and the social decision-making network, specifically related to sensory contextualization (BLA, LS, RSC, dHIP), saliency processing (LS, ACC, VTA), and memory formation and retrieval (vHIP, dHIP, mPFC, RSC).

Changes in functional connectivity after 24 hr after the onset of cohabitation may result from familiarization of a new spatial context that also implies a novel social context, that is exposure to new housing and to a novel, opposite-sex conspecific stranger. During this process, subjects are exposed to novel stimuli and engage in their first socio-sexual interactions, forming new memories of the partner and after-bonding behaviors such as mate and territorial guarding may appear. As mentioned earlier in this section, functional connectivity changes detected at 2 weeks after the onset of cohabitation may be related to long-term modulation of behavior as a result of social bonding, in which the partner and its associated cues become salient/rewarding and the pair-bond enters a maintenance phase. In general, these temporally dynamic functional connectivity changes may be related to the modulation of socio-sexual interactions with the partner and add a level of complexity that can be glimpsed due to the longitudinal analysis of the data. While the nodes here identified are interconnected into a larger network and their precise contribution to complex social behavior remains elusive to our methods, their potential role will be discussed based on previous results and the connectivity patterns here identified.

The lateral septum (LS) strikes as a relevant region in the modulation of social behavior, appearing as a hub that may integrate or relay spatial, contextual, and reward information between the ventral striatum and limbic regions, with cortical areas and hippocampal structures (Wirtshafter and Wilson, 2020). The hippocampus may also participate as an important relay region, being anatomically connected to the DG, RSC (Sugar et al., 2011), VTA (Gasbarri et al., 1994), and the prefrontal cortex (Preston and Eichenbaum, 2013). Twenty-four hours after the onset of cohabitation, the LS showed a decrease of functional connectivity with the dHIP; in mice, silencing this circuit decreases social aggression (Leroy et al., 2018). At the same period, an increase of connectivity between the LS and the RSC was also identified, which may be related to contextual memory acquisition (Opalka and Wang, 2020) and possibly associate social information with a spatial context. This change of connectivity may involve relay from the hippocampus, since anatomical projections between them are known (Risold and Swanson, 1997). In regard to the VP, it has been reported that its stimulation increases c-fos expression in regions including the LS and could be mediating reward processes (Panagis et al., 1997). The increase of connectivity between the latter regions seems to be on a relatively short term.

Social interaction would necessarily require new memory consolidation and, possibly, formation of the neural representation of the partner. Our data showed functional connectivity changes occurred in the hippocampus 24 hr after the onset of cohabitation (ACC–vHIP, LS–dHIP; Figure 5), which has demonstrated to be critical in encoding spatial and mnemonic information in rodents (Lin et al., 2017). ACC–vHIP activity has been associated with the regulation of fear in novel environments and in contextual fear generalization in rodents (Bian et al., 2019). However, the vHIP has also shown to contribute in social memory processing (Chiang et al., 2018), and the ACC is proposed to integrate new information to modify future behavior (Kolling et al., 2016). Increased ACC–vHIP connectivity could be related to social memory formation, but it may also be related to higher anxiety triggered by a new spatial and social environment.

In regard to functional connectivity changes 2 weeks after the onset of cohabitation, a long-term increase of connectivity was detected between the LS–mPFC, while a decrease was found in ACC–LS connectivity. It is enticing to propose that the LS is necessary for conspecific recognition and subsequent modulation of behavior, since AVP administration enhances pair-bond formation in prairie voles (Liu et al., 2001), and male voles exposed to females have increased Fos-immunoreactive cells in the LS (Wang et al., 1997). In rats, OXT release in the LS is also required for the maintenance of social memory, modulated by the relevance of the social stimulus (Lukas et al., 2013). Consequently, LS connectivity changes might be related to long-term social recognition and behavior modulation, mediated by afferents from prefrontal structures, that is the mPFC and the ACC (ACC–LS–mPFC). The aforementioned ACC may have a role in saliency regulation and in decision-making processes, while the mPFC is proposed to have a role in the initiation, maintenance, and modulation of social attachment and behavior in rodents (Ko, 2017; Tanimizu et al., 2017).

We also observed a decrease of connectivity between NAcc and VTA by the 2 week period. Considering the NAcc is densely innervated by VTA DA neurons (Steinberg et al., 2014), the after-bonding decrease of connectivity between these two regions may reflect the plasticity events necessary for pair-bond maintenance, in which D1-like receptor upregulation induces a long-term behavioral state in which stranger voles are attacked or avoided (Aragona et al., 2006). In mice, VP neuron activity drives positive reinforcement through projections to the VTA (Faget et al., 2018). The long-term increase of functional connectivity between these two regions suggests they may play a role in the maintenance of the bond through positive reinforcement.

Functional connectivity between the mPFC and the dHIP was increased. The dorsal region of CA2 has been characterized as a hub for sociocognitive memory processing in rodents, specifically for social memory encoding, consolidation, recall (Meira et al., 2018), and social recognition (Hitti and Siegelbaum, 2014). Also, the mPFC has been proposed as a critical region for remote memory retrieval and for memory consolidation, reportedly relying on the hippocampus (Euston et al., 2012). Therefore, their increased interaction could support the integration of new memories into pre-existing networks (Preston and Eichenbaum, 2013) and, potentially, in the incorporation of the partner into long-term memory representations. Indeed, OXT receptor activation in dHIP promotes the persistence of long-term social recognition memory in mice (Lin et al., 2018).

Lastly, we detected a decrease of BLA–NAcc connectivity between 24 hr and 2 weeks of cohabitation. It has been proposed that these regions regulate social reward processing in rats (Ambroggi et al., 2008) and may facilitate pair-bonding in prairie voles (Horie et al., 2020). Our results suggest that these regions may encode the reward-predictive cue, that is the partner, during bond establishment, but other regions could be involved in the maintenance of the bond, such as VP–VTA and NAcc–VTA.

Although both AVP and OXT have been proposed to have sex-differentiated roles in pair-bond formation (Young and Wang, 2004), the functional connectivity here explored showed no significant differences between sexes or significant sex-session interactions. Sexual dimorphism may still be present in a structural level, but the outcome of brain functional connectivity could remain similar between sexes.

Even though our results show ample individual variability in display of behavior as much as in brain functional connectivity, our analysis was able to capture group-wise consistent changes for both female and male prairie voles, suggesting these could be crucial in modulating social behavior after cohabitation with mating independently of the strength of the pair-bond. We propose that the long-term functional connectivity changes observed in the network are related to social bonding and may lead to pair-bonding induction and maintenance.

Concerning gICA, the obtained FWE-corrected P statistics found no significant differences between sessions. This being a voxel-wise type of analysis that implies multiple comparisons with its respective correction, gICA maps may not be sensitive enough to detect punctual, region-specific changes in brain functional connectivity. While ICA in rsfMRI is a valuable method with advantageous properties for identifying spatio-temporal signal sources that are not well-characterized or well-understood, different methods of analyses such as NBS seem more suitable for the aim of our study, which was to investigate the relationship between social behavior and functional connectivity, as well as the exploration of changes of connectivity in specific regions as a consequence of social bonding.

Correlation of brain functional connectivity with partner preference

While consistent longitudinal changes in connectivity were present in our data, we wanted to explore whether some of these nodes had a relationship with the level of partner preference, which captures the strength of the pair-bond (Williams et al., 1992), which was evaluated between 48 and 72 hr after the onset of cohabitation.

Positively, our results highlight the role of functional circuits in individual variability of social behavior, in which regions previously associated with reward processing are most relevant in determining bonding behavior. These correlations are dynamic in time, and each may contribute to the encoding and behavioral outcome of the pair-bond.

LS–VP baseline connectivity positively predicted the level of partner preference. It has been reported that AVP receptor overexpression in VP induces partner preference in socially promiscuous meadow voles (Lim et al., 2004). Twenty-four hours after the onset of cohabitation, a negative correlation was found between BLA–NAcc connectivity and partner preference. Variation in OXT receptor density in the NAcc relates to partner preference in prairie voles (King et al., 2016). Two weeks after the onset of cohabitation, LS–RSC negatively correlated with partner preference. Previous work demonstrated that AVP receptor expression in the RSC predicts sexual fidelity and territorial behavior in male prairie voles (Ophir et al., 2008) and may reflect different mating tactics, that is ‘wanderer’ or ‘resident’, in which resident voles maximize mating success through mate guarding and reduced territory space, in contrast to wanderers that overlap many home ranges to increase opportunistic mating (Getz et al., 2005). Our findings suggest that the baseline state of reward circuits (LS–VP) may influence initial bond development that will be encoded and established by BLA–NAcc, but long-term pair-bond maintenance may be modulated by social-reward, context-related regions (LS–RSC). In general, network functional connectivity variability may also reflect behavioral diversity influenced by genetic and environmental factors.

The prairie vole has proven to be a valuable model that has enabled the characterization of the neurobiological mechanisms behind complex socio-sexual behaviors, potentially useful to understand human social bonding and its alterations in psychological disorders. To our knowledge, this is the first study to demonstrate correlations between social bonding with brain functional connectivity, as well as time-dependent changes of multiple interacting brain regions (networks) of male and female prairie voles in cohabitation with mating. However, there are some limitations related to the use of rsfMRI. First, the animals were scanned under anesthesia, which may potentially alter the brain functional connectivity. Yet, we have used an anesthesia protocol that practically preserves the functional interactions as in the awake state (Grandjean et al., 2014; Paasonen et al., 2018) and that has not been reported to significantly alter rsfMRI and BOLD response in longitudinal sessions (Adamczak et al., 2010; Weber et al., 2006; Zhao et al., 2008). Although functional neuroimaging of awake prairie voles is possible, the acclimation process may induce significant stress and even increase the risk of physical harm (Yee et al., 2016), limiting the longitudinal design and the interpretation of the results.

Second, unlike male voles, female voles were subjected to surgical procedures and hormonal treatment before the MR scanning sessions. Although male voles are not precise controls for females, no significant differences were identified between male and female voles, behavior-wise or in rsfMRI connectivity. Also, since no significant longitudinal changes in functional connectivity were found in brain structures not expected to be involved in social behavior, the observed changes may be mostly related to social bonding and behavior. However, a control group with no behavioral protocols and same hormonal conditions as the ones reported in this study is desirable to confirm that there are no confounding effects in the longitudinal analysis.

Third, the precision in the definition of the ROIs is limited by the shape and the size of the voxels, being hard to assure that the defined regions specifically and uniquely include signal from the anatomical ROIs. Nevertheless, we are certain that the defined ROIs include the actual ROIs, and at the least, the changes here reported are related to those areas and their surrounding tissue. It is important to note that functionally connected nodes may not necessarily have direct axonal projections to each other, and the interaction between nodes could be mediated or relayed through other structures (Friston et al., 1993). Also, contrary to other electrophysiological or neuropharmacological methods, functional MRI captures an indirect measure of neuronal activity (Kim and Ogawa, 2012). However, it poses the great advantage that it allows the longitudinal exploration of the brain, with the best spatial resolution and wider coverage that any imaging method can achieve non-invasively. The consistency of our results with recent findings using direct electrophysiology readings (Amadei et al., 2017) further supports their relevance in identifying the neurophysiology of complex social behaviors. With evidence from previous reports as support, ROIs analyzed in this study constitute at least three different networks and have multiple roles in the development, establishment, and maintenance of social bonding. Additional research is required to increase understanding on their precise role in prairie vole behavior.

In conclusion, our findings suggest the existence of several functional brain networks involved in the modulation of social bonding in the prairie vole. Network-based statistics revealed a network involving cortical and striatal regions in which its functional connectivity predicted the onset of affiliative behavior, and another overlapping network was associated with the amount of social interaction after pair-bonding. Furthermore, a network with significant changes in time was revealed, with the lateral septum playing a role as a central hub that connects prefrontal and cortical regions with regions from the limbic system and the ventral striatum. Additionally, the retrosplenial cortex-lateral septum, lateral septum-ventral pallidum, and basolateral amygdala-nucleus accumbens functional connectivity correlate with the level of partner preference. In summary, brain functional connectivity allows exploring the mechanisms that underlie individual variability in the expression of socio-sexual behavior, enabling its prediction, and sexual experience and long-term cohabitation induce network-wide changes in socio-sexual relevant circuits. Overall, these results provide a novel approach and analyses to further investigate the neurophysiology of complex social behaviors displayed in the prairie vole.

Materials and methods

Key resources table
Reagent type
(species) or resource
DesignationSource or referenceIdentifiersAdditional information
Software, algorithmNetwork-Based Statistic ToolboxPMID:20600983RRID:SCR_002454
Software, algorithmNetwork-Based R-Statisticsdoi: 10.1101/2020.11.07.373019RRID:SCR_019114
Software, algorithmnlmeCRANRRID:SCR_015655
Software, algorithmFSLPMID:21979382RRID:SCR_002823
Software, algorithmANTS – Advanced Normalization ToolSPMID:20851191RRID:SCR_004757
Software, algorithmggplot2CRANRRID:SCR_014601
Software, algorithmGraphPad PrismGraphPadRRID:SCR_002798
Software, algorithmMATLABMathWorksRRID:SCR_001622
Chemical compound, drugDexmedetomidineZoetisPubChem CID:
5311068
s.c. bolus dose:0.05 mg/kg
Chemical compound, drugIsofluranePiSAPubChem CID: 3763
Chemical compound, drugβ-Estradiol 3-benzoateSigmaPubChem CID: 222757Dose: 0.5 mg/mL dissolved in corn oil
OtherDow Corning Silastic Laboratory TubingThermoFisher ScientificLength of tubing for capsule: 1.7 cm

Animals

Thirty-two 3-month-old sexually naïve female (N = 16) and male (N = 16) prairie voles (M. ochrogaster) were used in the study. The animals were housed in a temperature (23°C) and light (14:10 light–dark cycle) controlled room and provided with rabbit diet HF-5326 (LabDiet, St. Louis, MO) oat, sunflower seeds, and water ad libitum. These voles were previously weaned at 21 days, housed in same-sex cages, and were descendants of voles generously donated by Dr. Larry J. Young from his colony at Emory University. The number of subjects per group is comparable to the largest sample sizes using rsfMRI in rodents (Bajic et al., 2016; Christiaen et al., 2019; Grandjean et al., 2014). All surgical, experimental, and maintenance procedures were carried out in accordance with the ‘Reglamento de la Ley General de Salud en Materia de Investigación para la Salud’ (Health General Law on Health Research Regulation) of the Mexican Health Ministry that follows the National Institutes of Health’s ‘Guide for the Care and Use of Laboratory Animals’ (NIH Publications No. 8023, revised 1978). The animal research protocols were approved by the bioethics committee of the Instituto de Neurobiología, UNAM.

Surgical procedures

Request a detailed protocol

Fourteen days before the experimental protocol, female voles were bilaterally ovariectomized. After recovery, silastic capsules (Dow Corning Silastic Laboratory Tubing; ThermoFisher Scientific, Pittsburg, PA) containing estradiol benzoate (E2B; Sigma–Aldrich, St. Louis, MO) dissolved in corn oil (0.5 mg/mL of E2B) were implanted via s.c. to induce sexual receptivity 4 days before cohabitation protocol and remained implanted during the entire experimental protocol. Females in Microtus species are induced ovulators and do not show cyclic changes (Taylor et al., 1992), and for the purpose of this experiment, estrogen-in-oil capsule implantation allowed a stable and equal hormonal dose in all female subjects. The corresponding dose and procedure reliably induced sexual receptivity (Ingberg et al., 2012).

Anesthesia for image acquisition

Request a detailed protocol

Animals were anesthetized to avoid stress and excessive movement during scanning sessions. Isoflurane at 3% concentration in an oxygen mixture was used for induction and positioning in the scanner bed, in which the head was immobilized with a bite bar and the coil head holder. Once voles were securely placed in the scanner bed, isoflurane anesthesia (Sofloran; PiSA, Mexico) was adjusted at a 2% concentration and a single bolus of 0.05 mg/kg of dexmedetomidine (Dexdomitor; Zoetis, Mexico) was administered subcutaneously. Five minutes after the bolus injection, isoflurane anesthesia was lowered and maintained at 0.5%. MRI acquisition started when physiological readings were stabilized (~15 min after bolus injection). The use of both anesthetics has been reported as optimal for rsfMRI acquisition in rodents and closely resemble an awake condition (Grandjean et al., 2014; Paasonen et al., 2018), yet the mentioned combination of dose and administration route were previously standardized for this specific protocol in prairie voles (unpublished data). Body temperature was maintained with a circulating water heating pad within the scanner bed, respiration rate was monitored with an MR-compatible pneumatic pillow sensor, and blood oxygen saturation was measured with an MR-compatible infrared pulse-oximeter (SA Instruments Inc, Stony Brook, NY). After the scanning sessions, animals were monitored until fully recovered and transferred back to their housing.

Image acquisition

Request a detailed protocol

MRI acquisition was conducted with a Bruker Pharmascan 70/16US, 7 Tesla magnetic resonance scanner (Bruker, Ettlingen, Germany), using an MRI CryoProbe transmit/receive surface coil (Bruker, Ettlingen, Germany). Paravision-6 (Bruker, Ettlingen, Germany) was used to perform all imaging protocols. Before running the fMRI sequence, local field homogeneity was optimized within an ellipsoid covering the whole brain and skull using previously acquired field maps. rsfMRI was acquired using a spin-echo echo-planar imaging (SE-EPI) sequence: repetition time (TR) = 2000 ms, echo time (TE) = 19 ms, flip angle (FA) = 90°, field of view (FOV) = 18 × 16 mm2, matrix dimensions = 108 × 96, yielding an in-plane voxel dimensions of 0.167 × 0.167 mm2, and slice thickness of 0.7 mm, total volumes acquired = 305 (10 min and 10 s). EPI bandwidth was 288,461.5 Hz, the number of slices was 25, pi pulse (refocusing pulse) duration was 1.5455 ms with 2200.0 Hz bandwidth, and pi/2 (excitation pulse) was 1.9091 ms and 2200.0 Hz, respectively. After the rsfMRI sequence, an anatomical scan was obtained using a spin-echo rapid acquisition with refocused echoes (Turbo-RARE) sequence with the following parameters: TR = 1800 ms, TE = 38 ms, RARE factor = 16, number of averages (NA) = 2, FOV = 18 × 20 mm2, matrix dimensions = 144 × 160, slice thickness = 0.125 mm, resulting in isometric voxels of size 0.125 × 0.125 × 0.125 mm3.

Cohabitation and behavior analysis

Request a detailed protocol

Female and male voles unrelated to each other were randomly assigned as couples and placed for cohabitation in a new home cage with fresh bedding to promote ad libitum mating and social interaction. The first 6 hr of cohabitation was video recorded for subsequent analysis of social and mating behavior: mount, intromission, and ejaculation latencies from male voles; lordosis latency from females; and huddling latencies for each male and female pair. In male voles, mount latency was scored if it straddled the female from behind with pelvic thrusting, a continued mount behavior with repetitive thrusting was scored as intromission latency, and ejaculation latency was scored if after intromission and deeper thrusting an evident period of male inactivity was followed. Lordosis latency was scored if the female vole adopted an immobile posture with concave back flexion, neck extension, elevation of the hindquarters, and tail deviation to facilitate male mounting and intromission. Some females displayed lordosis reflex as pre-mounting behavior and were scored for lordosis latency. Since the onset of behaviors such as grooming or licking involve social contact and may overlap with huddling (Burkett et al., 2016), in this study huddling latency was only scored if there was continuous, side-to-side bodily contact for at least 10 s. Once joined, voles were housed in couples for the remaining of the experiment and were only separated for MRI scanning sessions and behavioral tests.

Partner preference test

Request a detailed protocol

Subjects underwent a 3 hr PPT to evaluate pair-bond formation. This protocol was based on a previously described test (Williams et al., 1992) and was performed in custom-built, three-chambered clear plastic arenas divided by perforated clear plastic barriers that allowed visual, auditory, and olfactory contact, but not physical interaction or mating behavior between subjects. In the central chamber, the vole being tested could roam freely, and time spent in the incentive areas at opposite sides of the chamber was recorded. The incentive areas were defined as the proximal space next to the chambers with its ‘partner’ or with an opposite-sex, novel ‘stranger’ vole. All stranger voles were unrelated to subjects in the test and had the same age and hormonal condition than the sexual partner. On each test, partner or stranger voles were randomly and alternately positioned on the opposite chambers of the arena. Each subject was tested only once in the PPT, either at 48 hr or at 72 hr after the onset of cohabitation. If a vole had PPT (assigned randomly and alternately) at the 48 hr period, its partner was tested at the 72 hr period to enable rest between tests and avoid excessive stress. PPT data analysis was performed with UMATracker software (Yamanaka and Takeuchi, 2018), which allowed quantification of the proportion of time spent with each of the stimulus voles. The percentage of time spent on the area related to the partner and of the area of the stranger were obtained. From this data, a partner preference index was calculated for each subject, consisting of the proportion of time spent on the area related to the partner divided by the proportion of time spent on both social incentive areas (time with partner plus time with stranger vole).

Imaging data pre-processing

Request a detailed protocol

Imaging data pre-processing was performed with FMRIB’s Software Libraries (FSL; Jenkinson et al., 2012) and Advanced Normalization Tools (ANTs; Avants et al., 2011) for spatial registration. To avoid initial signal instability, the first five volumes of each functional series were discarded. The slice-timing correction and motion correction were applied, using the first non-discarded volume as reference. The reference volume was also taken to determine the deformable transformation to the corresponding anatomical image. The resulting transformation was combined with a non-linear transformation to a prairie vole brain template obtained from previously published work (Ortiz et al., 2018Figure 1—figure supplement 3). Functional images were later warped to the brain template and resampled to a resolution of 0.16 × 0.16 × 0.7 mm3. To minimize physiological confounds, the first five eigenvectors (time series) within the combined non-gray matter mask were obtained (Behzadi et al., 2007), since recent findings have shown that vascular, ventricle, and white matter signal regression enhances functional connectivity specificity from rsfMRI data (Grandjean et al., 2020). These eigenvectors and the six motion parameters (three rotations, three displacements) were regressed out from each subject’s functional series. Datasets were band-pass filtered to retain frequencies between 0.01 and 0.1 Hz (Gorges et al., 2017). Finally, smoothing was applied with a box kernel with size of three voxels, using FSL.

Functional connectivity analysis

Request a detailed protocol

Regions of interest (ROIs) were manually defined on the anatomical prairie vole brain template, visually guided with the Allen Mouse Brain Atlas (Lein et al., 2007), to select the minimum possible number of voxels that included/covered the respective region. Connectivity matrices were calculated using MATLAB (Mathworks, Natick, MA) as follows: first, the average time series for each region were extracted from the pre-processed fMRI datasets and then, partial-correlation estimates were obtained for all possible pairs of ROIs, partialing out the remaining time series. Before analysis, correlation values were Fisher z-transformed. For each region, seed-based functional connectivity maps (Figure 1—figure supplement 4), as well as signal-to-noise ratio values, were generated for the assessment of data quality (Figure 1—figure supplement 5C).

Due to technical problems, two subjects missed the baseline MRI acquisition (session 1), and two subjects missed the 2 week cohabitation MRI acquisition (session 3). Additionally, two baseline datasets showed signal loss in the dorsal cortex (Figure 1—figure supplement 5). The final imaging sample consisted of 90 datasets, with only six subjects missing one session.

Correlations between behavioral data and functional connectivity were assessed with the NBS framework, using the NBS Toolbox (Zalesky et al., 2010), which estimates the statistical significance of clusters of connections (networks) by comparing their strength with a null distribution of such property obtained with 5000 permutations of the original data. Since the NBS Toolbox is restricted to general linear hypothesis testing, the longitudinal data was analyzed with linear mixed models (LMM), implemented via the Network-Based R-statistics package (NBR; Gracia-Tabuenca and Alcauter, 2020), which allows to implement LMM in an NBS framework. Also, NBR restricts the permutation of within-subject data based on the random variables. Specifically for functional connectivity longitudinal analysis, LMM were fitted for each connection (edge) considering sex (male, female) as a fixed variable, and session (baseline, 24 hr, and 2 weeks of cohabitation) and intercept as random variables. Those edges with significant (p<0.05) sex, session, or interaction effects at this level were used to identify sets of connected nodes, that is components or networks, and their statistical strength (sum of their statistical estimates, Z-values) were compared with a null distribution of the largest clusters from the statistical tests of the permutated data (5000 permutations), this being the NBS framework (Zalesky et al., 2010).

Group-independent component analysis

Request a detailed protocol

gICA was applied using FSL’s melodic (FSL; Jenkinson et al., 2012) on the same pre-processed images previously used in this study, only the female and male prairie voles that underwent all three MR scanning sessions were included on the analysis (n = 26), each subject having three different pre-processed images (n = 78). The number of components was set to 10, as described in previous work in prairie voles by our own group (Ortiz et al., 2018). To identify connected voxels over background, gICA maps were scaled to Z-scores and thresholded voxel-wise at Z ≥ 2.3 based on a Gaussian/Gamma mixture model and an alternative hypothesis testing approach. These gICA maps were visually inspected and labeled according to their anatomical distribution and location of their maximal regions, with additional visual guidance from the Allen Mouse Brain Atlas (Lein et al., 2007).

Statistical analysis

Request a detailed protocol

Data in the study is presented as mean ± standard error of the mean unless otherwise noted. Partner preference was explored with one-tailed Mann–Whitney U tests given that Shapiro–Wilk normality tests revealed that such data was not normally distributed. Behavioral data was analyzed with GraphPad Prism 5 (GraphPad Software, La Jolla, CA).

The NBS Toolbox was used with a level of significance of p<0.05 to identify networks (sets of connections) with significant linear associations with behavioral data. Correlation matrices were not thresholded for NBS analyses, but a fixed threshold of T > 1.7 was used to define the connections with significant associations, for both the actual test and the corresponding tests for the permuted data (5000 permutations). A posteriori Pearson’s two-tailed correlation tests were obtained for each edge within the significant networks and were depicted with R software package ggplot2 (ggplot2; Wickham, 2016) to better describe the relationships between specific connections and the associated behavior.

Implementation of LMM in NBR relies on the ‘nlme’ R package (nlme; Pinheiro et al., 2017) to fit LMM at the edge level and integrate them with the permutation tests required for the NBS framework. Since the degrees of freedom in LMM may vary for each variable (e.g., F(1,30) for Sex vs F(2,56) for Session), the edge-wise threshold was defined at p<0.05. Similar to the GLM-based NBS, no initial threshold was applied to the connectivity matrices prior to LMM-based NBS. To better describe the differences between sessions for each edge in the significant networks, estimated marginal means post hoc comparisons were used for LMM and were corrected for an FDR q < 0.05.

In order to evaluate whether networks found through gICA maps had significant changes between sessions, dual regression was applied on the spatial maps from the group-average analysis to generate subject-specific associated time series and the corresponding spatial maps (Beckmann et al., 2009; Nickerson et al., 2017). First, for each subject, the group-average set of spatial maps was regressed into the subject’s 4D space-time dataset, obtaining a set of subject-specific time series, one per group-level spatial map. Next, those time series were regressed into the same 4D dataset, resulting in a set of subject-specific spatial maps, one per group-level spatial map (n = 10). Group differences across sessions were assessed using FSL’s randomize permutation-testing tool (Winkler et al., 2014), in which two-sample paired t-tests were used (session 1 vs session 2, session 1 vs session 3, and session 2 vs session 3, respectively), given that an adequate anova version is not available in randomize for this experimental design.

Data availability

Data generated or analysed during this study are included in the manuscript and supporting files. Source data files have been provided for Figures 3 to 6 and supplementary source data has been provided for Figure 2. Code for Figure 5 is an R-based package available at https://cran.r-project.org/web/packages/NBR/index.html.

The following data sets were generated

References

  1. Software
    1. Pinheiro J
    2. Bates D
    3. DebRoy S
    4. Sarkar D
    5. Heisterkamp S
    6. Van Willigen B
    7. Maintainer R
    (2017)
    Package ‘Nlme’, Linear and Nonlinear Mixed Effects Models
    Package ‘Nlme’, Linear and Nonlinear Mixed Effects Models.
  2. Book
    1. Wickham H
    (2016)
    Ggplot2: Elegant Graphics for Data Analysis
    Springer-Verlag .

Decision letter

  1. Tali Kimchi
    Reviewing Editor; Weizmann Institute of Science, Israel
  2. Catherine Dulac
    Senior Editor; Harvard University, United States
  3. Noam Shemesh
    Reviewer; Champalimaud Centre for the Unknown, Portugal
  4. William Kenkel
    Reviewer; Georgia State University, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This work provides new understanding of the changes in brain functional connectivity that occur following pair bonding in the monogamous prairie vole, by analyzing resting state functional MRI (rsfMRI) before and after formation of pair bonding. The authors found that in both males and females partners there is a connectivity change in specific brain networks associated with pair bonding formation. Moreover, it was found that partner preference evaluated at 48 or 72 hours of cohabitation, may predict long-term functional brain connectivity.

Decision letter after peer review:

Thank you for submitting your article "Brain functional connectivity modulates social bonding in monogamous voles" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Catherine Dulac as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Noam Shemesh (Reviewer #1); William Kenkel (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, when editors judge that a submitted work as a whole belongs in eLife but that some conclusions require a modest amount of additional new data, as they do with your paper, we are asking that the manuscript be revised to either limit claims to those supported by data in hand, or to explicitly state that the relevant conclusions require additional supporting data.

Our expectation is that the authors will eventually carry out the additional experiments and report on how they affect the relevant conclusions either in a preprint on bioRxiv or medRxiv, or if appropriate, as a Research Advance in eLife, either of which would be linked to the original paper.

Summary:

The study examines how functional connectivity in the brain of prairie voles changes as a consequence of the formation of pair bonding by analyzing resting state functional MRI (rsfMRI) in anesthetized voles before and after formation of pair bonding. The authors found that in both males and females partners there is a connectivity changes in specific brain networks (including the ACC, VTA, LS and the hippocampus) associate with pair bonding formation. Moreover, they found that partner preference evaluated at 48 or 72 hours of cohabitation, predicted long-term functional connectivity between the medial amygdala and the ventral pallidum.

All reviewers agreed that the findings are interesting and novel. Reviewer #1 noted that "The authors went to great length to scan a relatively large number of animals and also did a very good job with analyses of resting-state data, which can be tricky at times."

Few important concerns were raised which require to be dealt in the revised submission. Generally the three reviewers agreed that the essential revision could be dealt with in the text.

The revisions should include additional analyses, greater level of methodological detail, additional clarity and explanation of the measured behaviors and statistics, thorough discussion to address the issue of lack of control for the temporal dimension or other control groups. Also, all reviewers agreed that the title should be revised as it is misleading.

The reviews are included in their entirety below.

Reviewer #1:

My main major comments regard (1) the lack of controls for the temporal dimension; and (2) the lack of figures showcasing the quality of the data acquired (and potential processing pipeline effects). These could probably be addressed without any specific new experiments but would require quite a bit of new analyses.

1) Control for the 2w temporal dimension

a) While I understand that the surgery and capsule implantation are commonplace for behavioral or neurophysiological experiments, I am worried that they may introduce some temporal confounds that – to my understanding of the Materials and methods – have not been accounted for. How do we know that the some of the changes in female rsfMRI patterns do not covary with recovery from surgery / hormonal cycles? In an ideal scenario, a control group with no behavior would have been added. However, this can also be mentioned as a confounding factor in the Discussion, and the authors could probably reason against this possibility using prior neurophysiological data.

b) Along the same line, the authors should also present data from ROIs that are NOT expected to be related with social behavior, to control for some of these potential confounders.

2) Quality of experiments.

While again commending the authors for the large number of animals scanned for this study, I would need to see figures with raw data to estimate the levels of artifacts in the original images and how they were addressed. For example, were ghosts detected along the time series? How were these dealt with?

To assess all these, I'd expect the authors to include figures with (a) representative raw data from the rsfMRI scans (to show the quality of the raw data arising from the cryogenic coil); (b) time series analyses in ROIs before/after preprocessing; (c) corresponding maps of rs connectivity, along with connectivity matrices for the different brain areas for the different conditions.

These could be presented as supplementary data, although for me it makes sense to include in the main document.

Reviewer #2:

1) Can the choice of huddling latency be justified in more detail? I don't know whether a reader would appreciate the meaning of this measure.

2) On a similar note: partner preference indexes, being calculated by dividing the time spent huddling with the Partner over the total time spent huddling, seem like they might miss out on differences in total social contact, which could be interesting. That is, a vole with 5 seconds spent huddling with its Partner and 10 seconds total huddling would have an equivalent preference index (0.5) as a vole that spent 5 minutes huddling with its Partner and 10 minutes total huddling (0.5), but these would strike me as very different phenotypes. Total time spent huddling could be an interesting variable to consider alongside the latency to huddle and preference index.

3) While I appreciate the careful consideration taken to select ROIs for this study, a strength of whole brain neuroimaging often comes from its ability to survey all brain regions simultaneously and thereby discover regions previously not known to be involved. Social neuroscience has a list of "usual suspects" and it may be the case that other brain regions that have not yet been implicated in pair bonding could still be important. An exploratory analysis could therefore be useful to the field.

4) I appreciated the presentation of the findings. The time-dependent changes, particularly the acute plastic changes were interesting. I also appreciated the figures' presentation of individual data points and longitudinal trends. That said, the effects do not appear very large, since the distributions overlap so much. Of course, there could still be large effects with overlapping distributions if there was just large individual variation that stayed consistent from time point to time point, but that doesn't appear to be the case for most brain regions. Perhaps violin plots would show the distributions and effect sizes more clearly than box plots? This is something I have struggled with in the past as well, and finding a reliable data visualization solution isn't trivial.

5) Subsection “Partner preference measured between 48h and 72h of cohabitation predicts long-term functional connectivity between MeA and VP” – What was the strength of the correlation between partner preference and MeA-VP connectivity at two weeks of cohabitation? I think the reader would assume all brain region correlations were examined, but were all time points considered as well? I imagine there would be difficulties considering the large number of tests being run, but it also seems like measures from the baseline scan pre-bonding would be the most relevant?

6) Each partner in a pair was tested as the subject in a PPT. Were their preferences for one another correlated? That is, if a male showed a high preference index for his female, was that reciprocated?

7) The Materials and methods state "correlation values were Fisher's z-transformed" so the reader does not have a sense of how strong the connectivity truly is between two regions, rather, what is shown is how strong the connectivity is between two regions relative to the connectivity between two randomly selected regions. Somewhere the range of correlation strengths should be noted.

This manuscript does an admirable job distilling a large number of measures and results down into manageable take-home messages for the reader.

To my knowledge, there has not been much work done distinguishing the phases of 24 hours post introduction to 2 weeks post introduction. I think the reader would benefit from knowing why this time point was chosen. The subsection “Longitudinal changes in functional connectivity after pair bonding”, mentions mate and territorial guarding but these behaviors would have been present by 24 hours.

8) As shown in Figure 3A, the mPFC showed consistent declines in functional connectivity. Beyond my facetious first reaction that this shows romance does indeed impair higher order cognitive processes, I think this is an interesting finding. However, the authors then state "A long-term increase of connectivity was detected between the LS with the mPFC" which doesn't appear to agree with Figure 3A unless I'm misunderstanding something.

Reviewer #3:

1) The first and most fundamental is that there is no control group, and so interpretation of temporal changes associated with bonding assumes that repeated anesthesia and time had no effect on patterns of functional connectivity. Are there any studies that explicitly address this concern? They cite two studies suggesting their protocol "practically preserves the functional interactions as in the awake state", but this falls short of addressing the possibility that repeated anesthesia might change their measures. I really wish they had allocated some of their effort to examining a suitable control-group (sib-housed, age-matched adults, in which control females were OVX + E implanted).

2) The second is that use of NBS and its integration with the LMM is insufficiently described. I cannot tell how they did their analyses based on the provided description (subsection “Statistical analysis”).

3) Apparently the connectivity matrix that serves as an input to NBS is based on the partial correlation matrix at each time point. What was the threshold the authors used to define connectivity? This is essential for the NBS statistic but I could not find it in the Materials and methods or Results.

4) Subsection “Functional connectivity analysis”, they say that linear mixed models were fit, but they are not explicit about what the dependent variables are for these models. The following sentence says that networks "with significant sex, session or interaction effects were identified based on the Network Based Statistics." It is not at all clear what this means. NBS are tests for longer-than-random spans of significant links. Did they simply use NBS to identify significant spans of linkages, and then test each region that went into a significant span with an LMM? The following sentence makes it sounds as though the connections themselves were defined by significant effects of sex or time, not by correlations. The "Statistics" and "Functional connectivity" subsections of the Materials and methods need to be revised in order to be appropriately evaluated. Please be explicit about the relationship between the LMM and the NBS.

5) If sex and time were part of the NBS analysis, I would also ask the authors to be explicit about how they executed the permutation test of the NBS statistic. I assume they permuted the assignment of both sex and time. Did they preserve the structure of between-individual variation?

6) I have similar concerns about the relationship between NBS and correlations with behavior. I found the methods of analysis hard to follow.

7) I had some concerns about mating and bonding in the study. Regarding the estradiol-in-oil implant method, how long do females remain in behavioral estrus in this condition? Preventing pregnancy and extending estrus may distort the dynamics of bonding.

8) The authors only recorded mating behavior over 6h. On a related note, the strength of partner preference seems surprisingly weak (counting each subject as independent across a total 32 tests, but p only ~0.03-0.04); the authors claim an n=32, but tested 32 subjects twice (48h and 72h). In the subsection “Statistical analysis”, the authors say they used a Mann Whitney U, but in Figure 2, they seem to show results of the 2 partner-preference tests lumped together. The authors should be explicit about this analysis to ensure they have not pseudoreplicated their sample size. (If they included 64 tests from 32 individuals in 16 pairs, that would certainly seem to be a problem with pseudoreplication.) Moreover, some pairs did not mate in the 6 hour window. It seems possible that some of the effects of the 2 week vs. 24h time period might be due to the lack of mating in the initial 24h period. I wonder how the results would change if they excluded animals who didn't mate and/or showed no partner preference?

https://doi.org/10.7554/eLife.55081.sa1

Author response

Reviewer #1:

My main major comments regard (1) the lack of controls for the temporal dimension; and (2) the lack of figures showcasing the quality of the data acquired (and potential processing pipeline effects). These could probably be addressed without any specific new experiments but would require quite a bit of new analyses.

1) Control for the 2w temporal dimension

a) While I understand that the surgery and capsule implantation are commonplace for behavioral or neurophysiological experiments, I am worried that they may introduce some temporal confounds that – to my understanding of the Materials and methods – have not been accounted for. How do we know that the some of the changes in female rsfMRI patterns do not covary with recovery from surgery / hormonal cycles? In an ideal scenario, a control group with no behavior would have been added. However, this can also be mentioned as a confounding factor in the Discussion, and the authors could probably reason against this possibility using prior neurophysiological data.

We thank the reviewer for pointing out the need to discuss the potential confounds of the capsule implantation. First, we have extended the description of female prairie voles hormonal conditions to highlight the need of such implantation. The following text has been added in the Materials and methods subsection “Surgical procedures”:

“Females in Microtus species are induced ovulators and do not show cyclic changes (Taylor et al., 1992), and for the purpose of this experiment, estrogen-in-oil capsule implantation allowed a stable and equal hormonal dose in all female subjects. The corresponding dose and procedure has reliably induced sexual receptivity in rodents (Ingberg et al., 2012).”

Although male voles are not precise controls for females, we did not identify any significant longitudinal changes between females and males, nor in functional connectivity nor in behavioral tests, which could have been potentially attributed to the surgery. As suggested by the reviewer we have included a brief discussion within the limitations of the study, including the reference to the ideal scenario of a control group. The following text is now included:

“Second, unlike male voles, female voles were subjected to surgical procedures and hormonal treatment before the MR scanning sessions. […] However, a control group with no behavioral protocols and same hormonal conditions as the ones reported in this study is desirable to confirm there are no confounding effects in the longitudinal analysis.”

Finally, in order to highlight the relevance of the brain-behavior associations here presented as suggested by the editor, we now present behavioral associations with functional connectivity as the main results, then we show the longitudinal changes, highlighting those with significant associations with social behaviors, and discussing the potential limitations of our findings.

b) Along the same line, the authors should also present data from ROIs that are NOT expected to be related with social behavior, to control for some of these potential confounders.

We appreciate this suggestion. We tested an independent set of brain structures not expected to be involved in social behavior using the same data, i.e., connectivity matrices, with the same analysis methods described previously. As expected, no significant differences were identified.

The following text has been added in the Results section:

“To control for potential confounders, we selected an independent set of brain structures not expected to be involved in social behavior and were tested with the same data and analysis methods than the ROIs mentioned previously. […] Linear mixed models (LMM) implemented via the Network-Based R-statistics package (NBR; Gracia-Tabuenca and Alcauter, 2020) considering sex as a fixed variable, and session and intercept as random variables, found no significant networks with differences between sessions (pFWE=0.1568), nor sex*session interactions (pFWE=0.619 and pFWE=0.3554), suggesting that in both female and male prairie voles, cohabitation with mating and social bonding do not influence changes in functional connectivity in these regions between any session (baseline, 24h and 2 weeks after the onset of cohabitation).”

2) Quality of experiments.

While again commending the authors for the large number of animals scanned for this study, I would need to see figures with raw data to estimate the levels of artifacts in the original images and how they were addressed. For example, were ghosts detected along the time series? How were these dealt with?

To assess all these, I'd expect the authors to include figures with (a) representative raw data from the rsfMRI scans (to show the quality of the raw data arising from the cryogenic coil); (b) time series analyses in ROIs before/after preprocessing; (c) corresponding maps of rs connectivity, along with connectivity matrices for the different brain areas for the different conditions.

These could be presented as supplementary data, although for me it makes sense to include in the main document.

We thank the reviewer for these suggestions. We have included new figure supplements for Figure 1 to keep the main focus of the manuscript on main results, allowing the interested readers to check the data quality.

Reviewer #2:

1) Can the choice of huddling latency be justified in more detail? I don't know whether a reader would appreciate the meaning of this measure.

Thank you for the suggestion. We have extended the description of the relevance of this measure in the study of social behaviors. The following text has been added in the Materials and methods subsection “Cohabitation and behavior analysis”:

“The amount and proneness to affiliative behavior has shown reliability in describing the level of sociability and mating systems in Microtus voles. […] Other studies in voles have also measured huddling latency to assess levels of affiliative behavior, since it may influence pair bond induction and formation (Shapiro et al., 1990; Amadei et al., 2017).”

2) On a similar note: partner preference indexes, being calculated by dividing the time spent huddling with the Partner over the total time spent huddling, seem like they might miss out on differences in total social contact, which could be interesting. That is, a vole with 5 seconds spent huddling with its Partner and 10 seconds total huddling would have an equivalent preference index (0.5) as a vole that spent 5 minutes huddling with its Partner and 10 minutes total huddling (0.5), but these would strike me as very different phenotypes. Total time spent huddling could be an interesting variable to consider alongside the latency to huddle and preference index.

This is a very interesting observation. As a reminder, as mentioned in the “Partner Preference Test” subsection in Materials and methods, due to the physical restrictions placed on the chamber, instead of measuring time spent huddling, we measured the percentage of time spent in the incentive areas related to the partner or the stranger vole. Therefore, the partner preference index was calculated by dividing the percentage of time spent on the partner social incentive area over the percentage of time spent on both social incentive areas.

To address the particular question of variability in total social interaction during the test, the proportion of time spent on social incentive areas over the total recorded time was calculated for all subjects. The average percentage of time spent on social incentive areas over the total time of the test was 79.70 ± 2.65 %. Spearman correlation analyses showed that the amount of time spent on social incentive areas had no relationship with the partner preference index (rs(30)=0.08912 p=0.6277), (see new Figure 4—figure supplement 1). However, we tested if functional connectivity 24h after the onset of cohabitation had a relationship with the total time spent on social incentive areas during the test, regardless if it interacted with the partner or stranger prairie vole (N=32). Using the NBS toolbox, a significant positive linear relationship (p ≤ 0.013) was found with the following network component: LS-NAcc-mPFC-MeA-VP-MOB-DG (Figure 4A).The correlation strength between each connected node was obtained with a posteriori Pearson correlations: LS-NAcc (r(30) = -0.369, p = 0.037), NAcc-mPFC (r(30) = -0.312, p = 0.081), mPFC-MeA (r(30) = -0.373, p = 0.034), MeA-VP (r(30) = -0.376, p = 0.033), VP-MOB (r(30) = -0.531, p = 0.001), and MOB-DG (r(30) = -0.404, p = 0.021) (Figure 4B-G).

This findings suggest that the lower the connectivity between these regions 24h after the onset of cohabitation, the longer the subject would interact with a conspecific of the opposite sex (known or unknown) during the PPT, which was evaluated between 48 and 72h after the onset of cohabitation. Hence, it reflects a phenotype independent both from pre-bonding behavior and partner preference, since regardless of the preference for the partner, prairie voles of both sexes had a differential interest in engaging socially. We have added this new result in the Results section:

“During the partner preference test (PPT), the proportion of time spent on social incentive areas over the total recorded time was calculated for analysis in all subjects (see Materials and methods). […] The correlation strength for each connection in such network was obtained with a posteriori Pearson correlations: LS–NAcc (r(30) = -0.369, p = 0.037), NAcc–mPFC (r(30) = -0.312, p = 0.081), mPFC–MeA (r(30) = -0.373, p = 0.034), MeA–VP (r(30) = -0.376, p = 0.033), VP–MOB (r(30) = -0.531, p = 0.001), and MOB–DG (r(30) = -0.404, p = 0.021) (Figure 4B-G.”

The results are commented in the Discussion section:

A network component detected at 24h after the onset of cohabitation, LS–NAcc–mPFC–MeA–VP–MOB–DG, was negatively correlated with the amount of social incentive a subject would have during the Partner Preference Test, a test in which it interacted with two conspecifics of the opposite sex, one being the partner and the other a stranger vole. […] Apparently, its change in modulation is on a short term since a correlation at 2 weeks of cohabitation was not found and other regions or networks may play this role on a long-term basis.”

3) While I appreciate the careful consideration taken to select ROIs for this study, a strength of whole brain neuroimaging often comes from its ability to survey all brain regions simultaneously and thereby discover regions previously not known to be involved. Social neuroscience has a list of "usual suspects" and it may be the case that other brain regions that have not yet been implicated in pair bonding could still be important. An exploratory analysis could therefore be useful to the field.

Undoubtedly, exploratory analyses are useful for the discovery of regions not previously suspected to be involved in a particular behavior. To address the exploration of large-scale functional brain networks, we have included a gICA in the following sections:

Results section:

“Group independent component analysis (gICA) was also performed to address the exploration of large-scale functional brain networks and potential differences between sessions. […] However, the obtained FWE-Corrected P statistics found no significant differences between sessions. Being this a voxel-wise method, it may not be sensitive enough to detect punctual, region-specific changes in brain functional connectivity.”

Discussion section:

“Group independent component analysis (gICA) was also performed to address the exploration of large-scale functional brain networks and potential differences between sessions. […] Being this a voxel-wise method, it may not be sensitive enough to detect punctual, region-specific changes in brain functional connectivity.”

Materials and methods:

“gICA was applied using FSL’s melodic (FSL; Jenkinson et al., 2012) on the same pre-processed images previously used in this study, only the female and male prairie voles that underwent all 3 MR scanning sessions were included on the analysis (n=26), each subject having 3 different pre-processed images (n=78). […] These gICA maps were visually inspected and labeled according to their anatomical distribution and location of their maximal regions, with additional visual guidance from the Allen Mouse Brain Atlas (Lein et al., 2007).”

Materials and methods subsection “Statistical analysis”:

“In order to evaluate if networks found through gICA maps had significant changes between sessions, dual regression was applied on the spatial maps from the group-average analysis to generate subject-specific associated timeseries and the corresponding spatial maps (Beckmann et al., 2009; Nickerson et al., 2017). […] Group differences across sessions were assessed using FSL's randomise permutation-testing tool (Winkler et al., 2014), in which two-sample paired t-tests were used (session 1 vs. session 2, session 1 vs. session 3, and session 2 vs. session 3, respectively), given that an adequate anova version is not available in randomise for this experimental design.”

4) I appreciated the presentation of the findings. The time-dependent changes, particularly the acute plastic changes were interesting. I also appreciated the figures' presentation of individual data points and longitudinal trends. That said, the effects do not appear very large, since the distributions overlap so much. Of course, there could still be large effects with overlapping distributions if there was just large individual variation that stayed consistent from time point to time point, but that doesn't appear to be the case for most brain regions. Perhaps violin plots would show the distributions and effect sizes more clearly than box plots? This is something I have struggled with in the past as well, and finding a reliable data visualization solution isn't trivial.

Certainly, our data reflects ample individual variability and the effects of these findings may not appear very large. Nonetheless, partial correlation analyses were used to assess significant functional connections between any two regions minimizing the effects of other regions in such relation. This type of analysis has high sensitivity in the detection of true functional connections, as it first regresses out all other nodes before correlating any two nodes (Smith et al., 2013). Hence, our analysis shows the significant functional connections that group-wise, were the strongest and most consistent in spite of underlying individual variability, the latter also allowing us to explain the observed variability in the expression of social behavior.

Following the advice of the reviewer, we are now presenting data of the mentioned figure (now Figure 5) in violin plots.

5) Subsection “Partner preference measured between 48h and 72h of cohabitation predicts long-term functional connectivity between MeA and VP” – What was the strength of the correlation between partner preference and MeA-VP connectivity at two weeks of cohabitation? I think the reader would assume all brain region correlations were examined, but were all time points considered as well? I imagine there would be difficulties considering the large number of tests being run, but it also seems like measures from the baseline scan pre-bonding would be the most relevant?

Under your suggestion, we have decided instead to analyze partner preference correlations with regions that underwent significant longitudinal changes after social bonding to better address such associations. We have added this new result in the Results section:

“Although several regions were found to change significantly after cohabitation and social bonding, we tested if any specific connection of the detected network component (Figure 5A) could have a relationship with the partner preference index obtained between 48 and 72 hours after the onset of cohabitation, which was used to evaluate the level of pair bonding in each subject. […] Also, significant negative relationships were found between the partner preference index and BLA–NAcc functional connectivity at 24h of cohabitation (session 2) (r(30) = -0.464, p = 0.0074), and in LS–RSC functional connectivity at 2 weeks of cohabitation (session 3) (r(28) = -0.437, p = 0.015 ), which may suggest that the lower the connectivity between these regions at the specified time points, the higher the partner preference index (Figure 6B-C).”

The results are further commented in the Discussion section:

“While consistent longitudinal changes in connectivity were present in our data, we wanted to explore if some of these nodes had a relationship with the level of partner preference, which captures the strength of the pair bond (Williams et al., 1992), which was evaluated between 48 and 72 hours after the onset of cohabitation.

[…] In general, network functional connectivity variability may also reflect behavioral diversity influenced by genetic and environmental factors.”

6) Each partner in a pair was tested as the subject in a PPT. Were their preferences for one another correlated? That is, if a male showed a high preference index for his female, was that reciprocated?

This is a very interesting question that we did not address before. Two-tailed Spearman correlation analysis revealed that partner preference was not reciprocated between male and female partners (rs(14)=-0.203, p=0.450), suggesting partner preference may not be majorly influenced by the partner’s perceived socio-sexual behavior (Figure 4—figure supplement 1).

7) The Materials and methods state "correlation values were Fisher's z-transformed" so the reader does not have a sense of how strong the connectivity truly is between two regions, rather, what is shown is how strong the connectivity is between two regions relative to the connectivity between two randomly selected regions. Somewhere the range of correlation strengths should be noted.

We have included the non-Fisher z-transformed partial-correlation connectivity matrices in the new Figure 1—figure supplement 2.

This manuscript does an admirable job distilling a large number of measures and results down into manageable take-home messages for the reader.

Thank you for your kind comment.

To my knowledge, there has not been much work done distinguishing the phases of 24 hours post introduction to 2 weeks post introduction. I think the reader would benefit from knowing why this time point was chosen. The subsection “Longitudinal changes in functional connectivity after pair bonding”, mentions mate and territorial guarding but these behaviors would have been present by 24 hours.

This time period was selected based on previous literature. Specifically, on studies that have shown neurophysiological and behavioral changes occurring by a two-weeks period that may be related to the maintenance of the pair bond. The following text has been added to include such information in the Introduction section:

“While prairie voles of both sexes already display changes in behavior by 24 hours of cohabitation with mating as a result of pair-bonding (Wang and Aragona, 2004; Williams et al., 1992), the NAcc, for example, has substantially more D1-like receptor binding 2 weeks after female exposure, relative to non–pair bonded males (Aragona et al., 2006), a phenomenon that’s also been observed in female voles (Resendez et al., 2016). […] While there is a possibility that these changes in behavior involve the interplay of broad networks, it has not been directly explored in vivo.”

Altogether, such evidence seemed relevant to explore changes in functional connectivity after cohabitation for two weeks.

8) As shown in Figure 3A, the mPFC showed consistent declines in functional connectivity. Beyond my facetious first reaction that this shows romance does indeed impair higher order cognitive processes, I think this is an interesting finding. However, the authors then state "A long-term increase of connectivity was detected between the LS with the mPFC" which doesn't appear to agree with Figure 3A unless I'm misunderstanding something.

This is an interesting observation and potential interpretation. However, between baseline and 24h after the onset of cohabitation, a non-significant decrease of connectivity is shown between LS-mPFC, but it increased significantly at the 2 weeks period (Session 3). Post-hoc analyses showed LS-mPFC connectivity had significantly increased when comparing Baseline with 2 weeks (t(56) = -2.856; p=0.0060), and also when comparing 24h vs. 2 weeks LS-mPFC (t(56) = -4.013; p=0.0001) (Figure 5H).

Reviewer #3:

1) The first and most fundamental is that there is no control group, and so interpretation of temporal changes associated with bonding assumes that repeated anesthesia and time had no effect on patterns of functional connectivity. Are there any studies that explicitly address this concern? They cite two studies suggesting their protocol "practically preserves the functional interactions as in the awake state", but this falls short of addressing the possibility that repeated anesthesia might change their measures. I really wish they had allocated some of their effort to examining a suitable control-group (sib-housed, age-matched adults, in which control females were OVX + E implanted).

Indeed, a number of studies have addressed the effects of sedation in fMRI longitudinal studies in rodents, specifically in forepaw stimulation. In rats, three fMRI sessions conducted in 5-day intervals, in which a 0.05 mg/kg medetomidine bolus injection plus a subcutaneous infusion of 0.1 mg/kg/h medetomidine were used, found no significant differences in BOLD signal change between sessions (Weber et al., 2006). In mice, repetitive medetomidine sedation with a dose of 0.3 mg/kg in combination with medetomidine infusion of 0.6 mg/kg/h in a 3 fMRI session protocol was also well tolerated and with reproducible BOLD response in forepaw stimulation (Adamczak et al., 2010).

While these studies suggest medetomidine effectiveness for fMRI longitudinal experiments and for rsfMRI (Zhao et al., 2008), the combined use of isoflurane with medetomidine has shown to have better cortical and partially preserved thalamo-cortical connectivity, having a higher resemblance to the awake condition in comparison to isoflurane and medetomidine alone (Grandjean et al., 2014, Paasonen et al., 2018). Even though a control group would certainly be more suitable for assessing the effects of anesthesia in our study, the relatively low dose used for our rodent model and the evidence shown on previous reports, support that our sedation protocol would not have a significant effect on the longitudinal follow-up. A brief sentence addressing previous studies on anesthesia has been included (subsection “Anesthesia for image acquisition”). In addition, we have now added the longitudinal analysis of regions not expected to be involved in social behavior, which had no significant changes between sessions.

However, as we now state in limitations (Discussion), a control group is desirable to confirm there are no confounding effects in the longitudinal changes two weeks after the onset of cohabitation. Accordingly, and also as suggested by the editor and reviewers, we now present the associations between functional connectivity and behavior as the main results, then we show the longitudinal changes, highlighting those with significant associations with social behaviors, discussing the potential limitations of our findings.

2) The second is that use of NBS and its integration with the LMM is insufficiently described. I cannot tell how they did their analyses based on the provided description (subsection “Statistical analysis”).

Edgewise LMM through NBS framework was reformulated as follows in the Materials and methods – subsection “Functional connectivity analysis”:

“Correlations between behavioral data and functional connectivity were assessed with the Network Based Statistics (NBS) framework, using the Network Based Statistic Toolbox (Zalesky et al., 2010), which estimates the statistical significance of clusters of connections (networks) by comparing their strength with a null distribution of such property obtained with 5000 permutations of the original data. […] Those edges with significant (p < 0.05) sex, session or interaction effects at this level were used to identify sets of connected nodes, i.e. components or networks, and their statistical strength (sum of their statistical estimates, Z-values) were compared with a null distribution of the largest clusters from the statistical tests of the permutated data (5000 permutations), this being the NBS framework (Zalesky et al., 2010).”

3) Apparently the connectivity matrix that serves as an input to NBS is based on the partial correlation matrix at each time point. What was the threshold the authors used to define connectivity? This is essential for the NBS statistic but I could not find it in the Materials and methods or Results.

As stated by Zalesky et al., 2010, one of the advantages of NBS is that no initial threshold is required for the connectivity matrices to be analyzed (although it is possible to do so), however it is critical to define a threshold to define the connections that surpass the statistical tests, which will be used to define the clusters of connections. Accordingly, we did not apply any connectivity threshold to the connectivity matrices, but we did apply a threshold at the sample inference level to define the connections that surpassed the statistical tests, and with these, define the clusters of connections or networks. For the behavioral correlates with NBS, using GLM, we used a threshold of T > 1.7. For the LMM-based NBR, we are able to define such a threshold in terms of the p-value, and it was set at p < 0.05 for each edge. Note that when using LMM the degrees of freedom between variables may vary (e.g., F(1,30) for Sex vs. F(2,56) for Session), for this reason we rely on significance instead of a static F-value in the LMM. Nevertheless, control of the family- wise error rate (FWER) of the defined networks (clusters of connections), is controlled irrespective of the threshold choice by the estimation of the null distributions.

The following information has been now included in the Materials and methods “Statistical analysis” subsection:

“The NBS toolbox was used with a level of significance of p < 0.05 to identify networks (sets of connections) with significant linear associations with behavioral data. […] To better describe the differences between sessions for each edge in the significant networks, estimated marginal means (EMM) post-hoc comparisons were used for LMM and were corrected for a false discovery rate (FDR) q < 0.05.

4) Subsection “Functional connectivity analysis”, they say that linear mixed models were fit, but they are not explicit about what the dependent variables are for these models. The following sentence says that networks "with significant sex, session or interaction effects were identified based on the Network Based Statistics." It is not at all clear what this means. NBS are tests for longer-than-random spans of significant links. Did they simply use NBS to identify significant spans of linkages, and then test each region that went into a significant span with an LMM? The following sentence makes it sounds as though the connections themselves were defined by significant effects of sex or time, not by correlations. The "Statistics" and "Functional connectivity" subsections of the Materials and methods need to be revised in order to be appropriately evaluated. Please be explicit about the relationship between the LMM and the NBS.

The dependent variable of the LMM models was the functional connectivity represented in the connectivity matrices (i.e., partial correlations between ROIs). As the reviewer correctly states, NBS tests for longer-than-random spans of significant links, but first we need to identify those significant links. That is done performing a statistical test in each possible connection (element of the matrix), be it using linear mixed models or the general linear model, as specified in any case. So, NBS first performs the statistical test in every element of the matrix, to identify the weighted or binary size of the clusters of links that surpass the test. Then, for every permutation of the original data, it performs the same tests, and identifies the size of the larger cluster to generate a null distribution. Finally, the size of clusters identified in the actual (real) case are compared with the null distribution of the largest random clusters, and the final family-wise p-value for each of them is estimated. So, depending on the data and the question being asked, the NBS framework can be applied with a variety of statistical tests. Shortly, associations between functional connectivity (from a single session) and behavioral data, are performed with the GLM; longitudinal data is analyzed with LMM. The great advantage of the NBS framework is that not only tests individual connections, but it identifies sets of connections, or networks.

We have reviewed the “Statistics” and “Functional connectivity” subsections that concern the NBS framework as follows:

“Correlations between behavioral data and functional connectivity were assessed with the Network Based Statistic (NBS) framework, using the Network Based Statistic Toolbox (Zalesky et al., 2010), which estimates the statistical significance of clusters of connections (networks) by comparing their strength with a null distribution of such property obtained with 5000 permutations of the original data. […] Those edges with significant (p < 0.05) sex, session or interaction effects at this level were used to identify sets of connected nodes, i.e. components or networks, and their statistical strength (sum of their statistical estimates, Z-values) were compared with a null distribution of the largest clusters from the statistical tests of the permutated data (5000 permutations), this being the NBS framework (Zalesky et al., 2010).”

“Statistical analysis. The NBS toolbox was used with a level of significance of p < 0.05 to identify networks (sets of connections) with significant linear associations with behavioral data. Correlation matrices were not thresholded for NBS analyses, but a fixed threshold of T > 1.7 was used to define the connections with significant associations, for both the actual test and the corresponding tests for the permuted data (5000 permutations). […] To better describe the differences between sessions for each edge in the significant networks, estimated marginal means (EMM) post-hoc comparisons were used for LMM and were corrected for a false discovery rate (FDR) q < 0.05.”

5) If sex and time were part of the NBS analysis, I would also ask the authors to be explicit about how they executed the permutation test of the NBS statistic. I assume they permuted the assignment of both sex and time. Did they preserve the structure of between-individual variation?

That is correct, the labels for sex and session were also permuted when they were included in the model, but preserving the data within-subjects to maintain the data structure.

6) I have similar concerns about the relationship between NBS and correlations with behavior. I found the methods of analysis hard to follow.

NBS takes as input the connectivity matrices from all observations and tests their relationship to linear covariates using a general linear model (GLM). The hypothesis tested at every connection was if connectivity strength was linearly associated with behavioral data. Relying on permutations, the NBS framework identifies clusters of connections for which the null hypothesis was rejected. The NBS enables rejection of the null hypothesis at the network level (i.e., at the level of those clusters), controlling the family-wise error rate (FWER). We have modified the description as referred in the previous points.

7) I had some concerns about mating and bonding in the study. Regarding the estradiol-in-oil implant method, how long do females remain in behavioral estrus in this condition? Preventing pregnancy and extending estrus may distort the dynamics of bonding.

Silastic implants were selected over estradiol s.c. administration because daily injections for two weeks may have been a potential stressor on the subjects and its partners. This method also enables more stable plasma estradiol levels than pellets (Mosquera, Shepherd and Torraro, 2015). A pilot study was performed to determine the efficiency of the estradiol implant and females remain sexually receptive for at least 20 days (non-published data).

Regarding female hormonal condition, Microtus species are induced ovulators and do not show cyclic changes. Female prairie voles are sexually inactive in the family nest; and stimuli from the environment, particularly male-related olfactory stimuli, are thought to induce endocrine changes, including increases in serum estrogen levels that promote sexual and behavioral receptivity (Taylor et al., 1992). In the presence of an unrelated male, females usually become sexually receptive at 48 h, will mate for 24h (Carter et al., 1987), and may remain in vaginal estrous for as long as 30 days when housed across a barrier from males (Richmond and Conaway, 1969). Therefore, the estradiol-capsule implantation enabled female sexual receptivity and reduced defensive aggression towards the male partner (Winslow et al., 1993; Insel et al., 1995). A brief sentence has been added to remind the reader the particular estrus condition in female prairie voles, in the Materials and methods subsection “Surgical procedures”:

“Females in Microtus species are induced ovulators and do not show cyclic changes (Taylor et al., 1992) , and for the purpose of this experiment, estrogen-in-oil capsule implantation allowed a stable and equal hormonal dose in all female subjects.

Additionally, mating is not essential for partner preference formation in female voles, though it has shown to accelerate bonding. Rather, cohabitation leads to pair bonding, and events associated with mating or estrus induction facilitate the onset of partner preference; and pair bonding may be established with at least 6 hours of cohabitation with mating or 24h of cohabitation without mating (Williams et al., 1992).

8) The authors only recorded mating behavior over 6h. On a related note, the strength of partner preference seems surprisingly weak (counting each subject as independent across a total 32 tests, but p only ~0.03-0.04); the authors claim an n=32, but tested 32 subjects twice (48h and 72h). In the subsection “Statistical analysis”, the authors say they used a Mann Whitney U, but in Figure 2, they seem to show results of the 2 partner-preference tests lumped together. The authors should be explicit about this analysis to ensure they have not pseudoreplicated their sample size. (If they included 64 tests from 32 individuals in 16 pairs, that would certainly seem to be a problem with pseudoreplication.) Moreover, some pairs did not mate in the 6 hour window. It seems possible that some of the effects of the 2 week vs. 24h time period might be due to the lack of mating in the initial 24h period. I wonder how the results would change if they excluded animals who didn't mate and/or showed no partner preference?

Each subject was tested only once in the partner preference test (PPT), within a 48h to 72h time frame. If a vole had PPT (assigned randomly and alternately) at the 48h period, its partner was tested at the 72h period. Specifically, each vole participated in two different PPT tests: once as the tested vole, and once as a stimulus vole in its partner’s PPT. Hence, we obtained a single partner preference measure for each subject. This 48h-72h testing period was chosen due to time constraints and to avoid excessive stress on the animals, which have no access to food and water during the 3-hour duration of the test. For better clarity, we have added the following text in the Materials and methods subsection “Partner Preference Test”:

“Each subject was tested only once in the partner preference test (PPT), either at 48 or 72h after the onset of cohabitation. If a vole underwent PPT (assigned randomly and alternately) at the 48h period, its partner was tested at the 72h period to enable rest between tests and avoid excessive stress.”

Regarding the strength of partner preference, as it has been previously characterized in this rodent model, we did find a statistically significant preference for the partner. However, as discussed in the manuscript, there is also a wide variability in partner preference in this species, which gave us the advantage to explore correlations with functional connectivity.

In respect to mating behavior, as mentioned before, mating is not essential for pair bonding in female voles, but it may facilitate its establishment (Williams et al., 1992). Though it has been reported that males may pair with diestrus females and without copulation (Shapiro and Dewsbury, 1990), mating in male voles has reported to be more relevant in pair bond formation than in females (Winslow et al., 1993; Insel et al., 1995). This trend is not entirely clear in our data: of 5 male voles with a partner preference index lower than 0.5, 2 of them did not mate during the 6 h window. It is important to consider that these subjects may have succeeded in mating both before the 24h scanning session or before the PPT. Further investigation would be necessary to fully understand the effects of mating in male vole pair bonding facilitation and its relationship with changes in brain functional connectivity, since only a small number of subjects were found in this particular circumstance.

References:

Adamczak, J. M., Farr, T. D., Seehafer, J. U., Kalthoff, D., and Hoehn, M. (2010). High field BOLD response to forepaw stimulation in the mouse. NeuroImage, 51(2), 704–712.

Carter, C. S., Witt, D. M., Auksi, T., and Casten, L. (1987). Estrogen and the induction of lordosis in female and male prairie voles (Microtus ochrogaster). Hormones and Behavior, 21(1), 65–73.

Grandjean, J., Canella, C., Anckaerts, C., Ayrancı, G., Bougacha, S., Bienert, T., … Gozzi, A. (2020). Common functional networks in the mouse brain revealed by multi-centre resting-state fMRI analysis. NeuroImage, 205(October 2019).

Insel, Thomas R., Preston, S., and Winslow, J. T. (1995). Mating in the monogamous male: Behavioral consequences. Physiology and Behavior, 57(4), 615–627.

Mosquera, L., Shepherd, L., and Torrado, A. I. (2015). Comparison of Two Methods of Estradiol Replacement: their Physiological and Behavioral Outcomes. Journal of Veterinary Science and Technology, 06(06).

Richmond, M. E., and Conaway, C. H. (1969). Induced ovulation and oestrus in Microtus ochrogaster. Journal of Reproduction and Fertility, Supplement(6), 357–376.

Shapiro, L. E., and Dewsbury, D. A. (1990). Differences in affiliative behavior, pair bonding, and vaginal cytology in two species of vole (Microtus ochrogaster and M. montanus). Journal of Comparative Psychology, 104(3), 268–274.

Smith, S. M., Vidaurre, D., Beckmann, C. F., Glasser, M. F., Jenkinson, M., Miller, K. L., … Van Essen, D. C. (2013). Functional connectomics from resting-state fMRI. Trends in Cognitive Sciences, 17(12), 666–682.

Winslow, J. T., Hastings, N., Carter, C. S., Harbaugh, C. R., and Insel, T. R. (1993). A role for central vasopressin in pair bonding in monogamous prairie voles. Nature, 365(6446), 545–548.

https://doi.org/10.7554/eLife.55081.sa2

Article and author information

Author details

  1. M Fernanda López-Gutiérrez

    Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Data curation, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9513-5607
  2. Zeus Gracia-Tabuenca

    Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Software, Formal analysis, Visualization, Methodology
    Competing interests
    No competing interests declared
  3. Juan J Ortiz

    Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Supervision, Methodology
    Competing interests
    No competing interests declared
  4. Francisco J Camacho

    Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Supervision, Methodology
    Competing interests
    No competing interests declared
  5. Larry J Young

    Silvio O Conte Center for Oxytocin and Social Cognition, Center for Translational Social Neuroscience, Yerkes National Primate Research Center, Department of Psychiatry and Behavioral Sciences, Emory University, Atlanta, United States
    Contribution
    Resources, Funding acquisition, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Raúl G Paredes

    1. Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    2. Escuela Nacional de Estudios Superiores, Unidad Juriquilla, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Resources, Validation, Writing - review and editing
    Competing interests
    No competing interests declared
  7. Néstor F Díaz

    Instituto Nacional de Perinatología Isidro Espinosa de los Reyes, Ciudad de México, Mexico
    Contribution
    Resources, Funding acquisition, Writing - review and editing
    For correspondence
    nfdiaz00@yahoo.com.mx
    Competing interests
    No competing interests declared
  8. Wendy Portillo

    Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    portillo@unam.mx
    Competing interests
    No competing interests declared
  9. Sarael Alcauter

    Instituto de Neurobiología, Universidad Nacional Autónoma de México, Querétaro, Mexico
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    alcauter@inb.unam.mx
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8182-6370

Funding

Consejo Nacional de Ciencia y Tecnología (252756)

  • Wendy Portillo

Consejo Nacional de Ciencia y Tecnología (253631)

  • Raúl G Paredes

Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México (IN202818)

  • Wendy Portillo

Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México (IN212219-3)

  • Sarael Alcauter

Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México (IN203518-3)

  • Raúl G Paredes

Instituto Nacional de Perinatología (2018-1-163)

  • Néstor F Díaz

National Institutes of Health (P50MH100023)

  • Larry J Young

Consejo Nacional de Ciencia y Tecnología (626152)

  • M Fernanda López-Gutiérrez

National Institutes of Health (P51OD011132)

  • Larry J Young

Consejo Nacional de Ciencia y Tecnología (600922)

  • M Fernanda López-Gutiérrez

Instituto Nacional de Perinatología (3230-21216-05-15)

  • Néstor F Díaz

Dirección General de Asuntos del Personal Académico, Universidad Nacional Autónoma de México (IN208221)

  • Wendy Portillo

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Dr. Fernando A Barrios for his valuable comments to the manuscript, and Deisy Gasca, Martín García, Alejandra Castilla, Leopoldo González Santos, Nuri Aranda López, and Ma. De Lourdes Lara for their technical assistance. The authors thankfully acknowledge the imaging resources and support provided by the Laboratorio Nacional de Imagenología por Resonancia Magnética (LANIREM), part of the CONACYT’s network of national laboratories. This research was supported by grants CONACYT 252756 to WP, 253631 to RGP, UNAM-DGAPA-PAPIIT IN202818 and IN208221 to WP, IN212219-3 to SA, IN203518-3 to RGP, INPER 2018-1-163 and 3230-21216-05-15 to NFD. LJY’s contribution was supported by NIH grants P50MH100023 to LJY and P51OD011132 to YNPRC. MF López-Gutiérrez was supported by a fellowship from CONACYT (fellowship #626152) for her Masters of Science studies in Neurobiology at UNAM (Maestría en Ciencias [Neurobiología], UNAM), is a doctoral student from the Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional Autónoma de México (UNAM), and has received CONACyT fellowship (2020-000026-02NACF-17340, cvu600922).

Ethics

Animal experimentation: All surgical, experimental and maintenance procedures were carried out in accordance with the "Reglamento de la Ley General de Salud en Materia de Investigación para la Salud" (Health General Law on Health Research Regulation) of the Mexican Health Ministry which follows the National Institutes of Health's "Guide for the Care and Use of Laboratory Animals" (NIH Publications No. 8023, revised 1978). The animal research protocols were approved by the bioethics committee of the Instituto de Neurobiología, UNAM (Protocol 072). All fMRI scanning sessions were performed under a mixture of isoflurane and dexmedetomidine anesthesia, and all surgery was performed under sevoflurane or a mixture of ketamine/xylazine/saline anesthesia, with every effort to minimize suffering.

Senior Editor

  1. Catherine Dulac, Harvard University, United States

Reviewing Editor

  1. Tali Kimchi, Weizmann Institute of Science, Israel

Reviewers

  1. Noam Shemesh, Champalimaud Centre for the Unknown, Portugal
  2. William Kenkel, Georgia State University, United States

Publication history

  1. Received: January 11, 2020
  2. Accepted: January 11, 2021
  3. Accepted Manuscript published: January 14, 2021 (version 1)
  4. Version of Record published: January 29, 2021 (version 2)

Copyright

© 2021, López-Gutiérrez et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,778
    Page views
  • 282
    Downloads
  • 12
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. M Fernanda López-Gutiérrez
  2. Zeus Gracia-Tabuenca
  3. Juan J Ortiz
  4. Francisco J Camacho
  5. Larry J Young
  6. Raúl G Paredes
  7. Néstor F Díaz
  8. Wendy Portillo
  9. Sarael Alcauter
(2021)
Brain functional networks associated with social bonding in monogamous voles
eLife 10:e55081.
https://doi.org/10.7554/eLife.55081

Further reading

    1. Neuroscience
    Xiaosha Wang, Bijun Wang, Yanchao Bi
    Research Article Updated

    One signature of the human brain is its ability to derive knowledge from language inputs, in addition to nonlinguistic sensory channels such as vision and touch. How does human language experience modulate the mechanism by which semantic knowledge is stored in the human brain? We investigated this question using a unique human model with varying amounts and qualities of early language exposure: early deaf adults who were born to hearing parents and had reduced early exposure and delayed acquisition of any natural human language (speech or sign), with early deaf adults who acquired sign language from birth as the control group that matches on nonlinguistic sensory experiences. Neural responses in a semantic judgment task with 90 written words that were familiar to both groups were measured using fMRI. The deaf group with reduced early language exposure, compared with the deaf control group, showed reduced semantic sensitivity, in both multivariate pattern (semantic structure encoding) and univariate (abstractness effect) analyses, in the left dorsal anterior temporal lobe (dATL). These results provide positive, causal evidence that language experience drives the neural semantic representation in the dATL, highlighting the roles of language in forming human neural semantic structures beyond nonverbal sensory experiences.

    1. Neuroscience
    Ayako Yamaguchi, Manon Peltier
    Research Article Updated

    Across phyla, males often produce species-specific vocalizations to attract females. Although understanding the neural mechanisms underlying behavior has been challenging in vertebrates, we previously identified two anatomically distinct central pattern generators (CPGs) that drive the fast and slow clicks of male Xenopus laevis, using an ex vivo preparation that produces fictive vocalizations. Here, we extended this approach to four additional species, X. amieti, X. cliivi, X. petersii, and X. tropicalis, by developing ex vivo brain preparation from which fictive vocalizations are elicited in response to a chemical or electrical stimulus. We found that even though the courtship calls are species-specific, the CPGs used to generate clicks are conserved across species. The fast CPGs, which critically rely on reciprocal connections between the parabrachial nucleus and the nucleus ambiguus, are conserved among fast-click species, and slow CPGs are shared among slow-click species. In addition, our results suggest that testosterone plays a role in organizing fast CPGs in fast-click species, but not in slow-click species. Moreover, fast CPGs are not inherited by all species but monopolized by fast-click species. The results suggest that species-specific calls of the genus Xenopus have evolved by utilizing conserved slow and/or fast CPGs inherited by each species.