Fatherhood alters gene expression within the MPOA

Female parenting is obligate in mammals, but fathering behavior among mammals is rare. Only 3–5% of mammalian species exhibit biparental care, including humans, and mechanisms of fathering behavior remain sparsely studied. However, in species where it does exist, paternal care is often crucial to the survivorship of offspring. The present study is the first to identify new gene targets linked to the experience of fathering behavior in a biparental species using RNA sequencing. In order to determine the pattern of gene expression within the medial preoptic area that is specifically associated with fathering behavior, we identified genes in male prairie voles (Microtus ochrogaster) that experienced one of three social conditions: virgin males, pair bonded males, and males with fathering experience. A list of genes exhibiting different expression patterns in each comparison (i.e. Virgin vs Paired, Virgin vs Fathers, and Paired vs Fathers) was evaluated using the gene ontology enrichment analysis, and Kyoto Encyclopedia of Genes and Genomes pathways analysis to reveal metabolic pathways associated with specific genes. Using these tools, we generated a filtered list of genes that exhibited altered patterns of expression in voles with different amounts of social experience. Finally, we used NanoString to quantify differences in the expression of these selected genes. These genes are involved in a variety of processes, with enrichment in genes associated with immune function, metabolism, synaptic plasticity, and the remodeling of dendritic spines. The identification of these genes and processes will lead to novel insights into the biological basis of fathering behavior.


Introduction
Biparental care, in which both mother and father contribute to the care of the offspring, is displayed by a minority of mammalian species -usually cited as 3-5% [1][2][3]. Female parenting is obligate because mammalian offspring need to nurse, but the participation of the male is seen only in our own and a limited number of other mammalian species [4][5][6][7][8]. In species where paternal care does exist, including humans, it is often crucial to the survivorship of offspring; or at least has significant long-term impacts on growth as well as neural, reproductive, and social development [9,10]. However, much is still unknown about the specific hormonal and neurobiological regulation of paternal care [6].
The vast majority of parenting research focuses on the mother, while the role of the father has mostly been considered in the context of paternal absence [9,[11][12][13]. Considering paternal care through the absence of the father in a biparental species has drawbacks, however, since it confounds the quantitative absence of another caregiving individual with the qualitative absence of the father in particular. Paternal absence is the most extreme situation. Individual variation in fathering can also have long-term effects on offspring [9], and in the context of nonhuman mammals is always carried out in a biparental care situation. In prairie voles, we have shown that natural variation in biparental parenting behavior predicts pup development and juvenile social behavior [14], exploratory behavior and pair-bonding, and adult aggression and stress responses [15][16][17]. It is not always possible in an intact biparental family to disentangle which outcomes in offspring are due to maternal care and which are due to paternal care. However, some very interesting roles for the father have been observed. For instance, males may compensate for poor maternal care (or allow mothers to expend less energy on non-nutritive tasks like carrying) [14,18]; or a paternal behavior such as retrievals (carrying pups back to the nest or territory) may be directly linked to offspring display of retrievals and aggression as an adult [19,20].
The hormonal mechanisms underlying fathering behavior have been much less studied than those underlying maternal behavior, although it has been hypothesized that similar neural circuits are responsible for both maternal and paternal behaviors [21]. While alterations in neural activity with parenthood appear to be hormonally regulated in females, hormonal manipulation in males has often resulted in outcomes that are either ambiguous or species-specific [22]. For instance, testosterone is inversely related to paternal care in most species [9,23,24], but is obligate for paternal care in California mice [25]. Prolactin, another leading candidate for the regulation of male parenting, decreases male parenting when administered, as well as when blocked [26]. These inconsistencies have led some to suggest that across species, paternal behavior depends upon nonhomologous neuroendocrine circuits [6], and has raised the question of what factors are involved in the generation of these behaviors.
Although the hormonal regulation of parenting may vary by sex, it is believed that the neural circuit governing parental behavior is similar in mothers and fathers [27]. The medial preoptic area of the hypothalamus (MPOA) is a central node in the neural circuit that regulates both maternal and paternal care and has long been recognized as playing a critical role in the generation and regulation of parental behavior (see [9,28] for reviews). In biparental California mice, paternal experience increases Fos immunoreactivity in the MPOA [29]. Virgin male prairie voles that were exposed to pups also showed an increase in Fos immunoreactivity within the MPOA [30]. Lesions of the MPOA disrupt both maternal and paternal behavior in California mice [31]. In California mouse males, aromatase levels within the MPOA vary in response to parental status [32], while another study in male California mice showed a decrease in progesterone receptor mRNA expression in the MPOA of fathers compared with virgin males [33]. Studies in mice, which are not parental in the wild but can show induced paternal care in the laboratory [34][35][36], have bolstered the view of a central role for the MPOA in paternal care.
The goal of this study was to identify novel gene targets and potential mechanisms that may contribute to the production and regulation of paternal behavior. We analyzed gene expression in three groups of adult male prairie voles: virgin males, males who had formed a pair bond with a female, and males who had both pair bonded and gained fathering experience. Samples were taken from the MPOA, a region that is central to the expression of both maternal and paternal behaviors [21,[37][38][39], and RNA was extracted and sequenced.

Subjects
Subjects were 18 adult male prairie voles. Animals were born and housed in the Psychology Department Vivarium at the University of California, Davis. These animals were descendants of a wild stock originally caught near Champaign, IL. The animals were weaned at 20 days of age and pair housed with an animal of the same sex (sibling if available, similarly aged nonsibling if not) in small laboratory cages (27 Â 16 Â 13 cm) in which food and water were available ad libitum. All animals were maintained at $70 F (21 C) on a 14:10 light/dark cycle with the lights on at 6 a.m. All experiments were performed under National Institutes of Health guidelines for the care of animals in research and were approved by the Institutional Animal Care and Use Committee of the University of California, Davis.
At postnatal day (P) 42-45 subjects were placed in one of three groups of age-matched males: (i) virgin males, (ii) sexually experienced, 'pair-bonded' males, or (iii) males with fathering experience. This was designed to dissociate alterations in gene expression that were related to pair bonding from alterations related to paternal behavior. Virgin males were housed with a male same-age conspecific for $20 days, and they were euthanized without engaging in sexual contact with females. Pairbonded males were housed with a same-age female conspecific for $20 days, after which the males were euthanized. Because mating and pregnancy strengthens pair bonds in prairie voles [40][41][42][43], we confirmed that females were pregnant. Pair-bonded males were euthanized before females gave birth, ensuring they had no contact with pups. The third group consisted of males that had 3 days of paternal experience. These males were also housed with female pair-mates with whom they presumably formed a pair-bond. The females gave birth, and the males were permitted three days of contact with pups before they were euthanized. Three days of parental experience was chosen to minimize age differences between subjects. Furthermore, prairie vole fathers already exhibit large amount of paternal care by postnatal day 3 [14,44].
Subjects were anesthetized using isoflurane and euthanized via cervical dislocation. Upon euthanasia, brains were removed and flash frozen. The brains were sliced on a cryostat into 120 mm sections and mounted on slides. Punches were taken from the MPOA using a 15.5-gauge blunt needle (Fig. 1) and were stored in a À80 C freezer until RNA extraction.

RNA Extraction
Total RNA was isolated with Qiazol reagent (Qiagen, Valencia, CA) and purified with an RNeasy V R Plus Micro Kit (74004; Qiagen) as well as the optional DNase digestion (Qiagen 129046). A Nanodrop TM Spectrophotometer was used to determine the quality and the quantity of the RNA. All samples had a 260/280 ratio >1.8.

RNA Sequencing
A total of 18 RNA-seq libraries were prepared from the RNA of the 18 male prairie voles (Table 1). RNA sequencing and library preparation was performed by the DNA Technologies and Expression Analysis Core in the Genome Center of the University of California, Davis. A total RNA analysis ng sensitivity (Eukaryotes) of all 18 samples resulted in a mean RIN of 9.2 (range 8.3-9.9). Barcoded RNA-seq libraries were generated from 1 lg total RNA each after poly-A enrichment using the Kapa Stranded RNA-seq kit (Kapa Biosystems, Cape Town, South Africa) following the instructions of the manufacturer. The libraries were generated on a Sciclone G3 liquid handler (Caliper Life Sciences, Alameda, CA). Quality was verified with the Bioanalyzer 2100 instrument (Agilent, Santa Clara, CA) and quantified by fluorometry on a Qubit instrument (Life Technologies, Carlsbad, CA) and pooled in equimolar ratios. The pooled library was then quantified by qPCR with a Kapa Library Quant kit (Kapa Biosystems) and sequenced on one lane of an Illumina HiSeq 4000 (Illumina, San Diego, CA) with paired-end 150 bp reads.
Raw sequencing data have been deposited at NCBI's sequence read archive under study accession number SRP128134.

Bioinformatic Analysis
Bioinformatic analysis was performed by the UC Davis Bioinformatics Core Facility, also in the Genome Center. Briefly, reads were trimmed for adapter contamination and quality using scythe (version c128b19) and sickle (version 7667f147e6), respectively. The reads were then aligned to the prairie vole genome (MicOch1.0) using bwa mem (version 0.7.13), after which featureCounts (version 1.5.0-p1) was used to create the raw gene expressions counts. Finally, R (version 3.3.2) with the edgeR and limma/voom packages were used to filter and transform (voom transformation), and test for statistical significances between groups. Briefly, data were prepared by first choosing to keep genes that achieved at least 0.5 count per million in at least five samples, normalization factors were calculated using trimmed mean of M-value, and the voom transformation was applied. A completely randomized design was implemented, comparisons of interest were extracted using contrasts, and moderated statistics were computed using the Figure 1: A schematic representing the area from which tissue samples were taken. The circumference of the tissue punch is delineated by a circle, and the MPOA is outlined in black. The tissue punches removed the entirety of the MPOA, as well as small portions of adjacent hypothalamic tissue empirical Bayes procedure eBayes. Finally, each gene was corrected for multiple testing using the Benjamini-Hochberg false discovery rate correction. Gene expression was directly compared between each pair of groups, resulting in three comparisons: Virgin males vs Paired males (V vs P), Virgin males vs Fathers (V vs F), and Paired males vs Fathers (P vs F). No genes reached a statistical significant threshold (adjusted P-value <0.05) in any of the three pairwise comparisons.

Gene Ontology Analysis
In order to capture the genes that were most likely to show common functional differentiation between comparison groups we performed gene set enrichment analysis with Gene Ontologies (cellular component, molecular function, or biological process). Enrichment testing was conducted using the Kolmogorov-Smirnov test as implemented in the Bioconductor package topGO [45]. We next examined the gene ontology (GO) annotations that were significantly enriched (P-value < 0.05) and filtered the GO annotations in each comparison that were related to the brain or behavior, excluding unrelated annotations (i.e. GO: 0003014, Renal system process or GO: 0008354, Germ cell migration). We then categorized the remaining annotations based on gross function within each comparison group.

Kyoto Encyclopedia of Genes and Genomes Pathways Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways offer insight into how genes interact within biological and metabolic processes. We identified individual genes associated with significantly enriched GO annotations and ran each individual gene through the KEGG pathways database (http://www.ge nome.jp/kegg/) to identify molecular signaling pathways associated with that gene. A single gene may be involved in a number of different pathways, so we then identified commonly recurring pathways associated with the individual genes. Pathways that were unrelated to brain function (e.g. those that were involved in kidney, liver, or heart metabolism) were not included in the analysis. When a specific gene was associated with multiple pathways of interest it was identified as a candidate gene for further analysis. For example, Grin2a was associated with six pathways that are involved in neural plasticity.
We identified 49 candidate genes across 9 biologically significant KEGG pathways and examined how their expression changed between conditions of social experience. For each gene, we standardized expression relative to virgin males. This allowed us to identify whether there were coordinated changes between genes that were involved in a specific pathway. We averaged gene expression across animals within each condition and transformed the data into ratios; the values we used for all analyses were the ratios of gene expression in each condition relative to the virgin condition. For each KEGG pathway, we determined the mean relative expression for all genes in each condition. The expression ratio of genes in virgin animals was set at 1, a value >1 indicated that genes were more expressed relative to virgins, and a value <1 indicated that genes were less expressed relative to virgins. Thus, we analyzed whether gene expression for each gene varied between paired animals and fathers. Effect size was measured using Cohen's d.

Assessment of Gene Interaction Networks
After generating a list of candidate genes, we analyzed the connectivity of the gene network using the STRING Database (string-db.org) [46]. The STRING database identifies protein-protein interactions between members of a gene set, which allows the user to build a network of functional gene interactions. STRING also measures the functional and interaction enrichments of the gene network, calling upon GO annotations, KEGG pathways, and connections between nodes.

NanoString Analysis
Following the identification of candidate genes, we performed a quantitative analysis of the expression of 33 genes (30 target genes and 3 housekeeping genes) using the nCounter SPRINT profiler (NanoString Technologies, Seattle, WA). Genes were chosen to be included in the NanoString analysis based on their log fold change values as determined by the expression data, as well as their functional significance. One additional gene, Bdnf, was chosen due to previous studies indicating that it plays a significant role in plasticity and parenting [47,48]. The nCounter analysis assay was conducted using RNA that remained after the completion of the sequencing experiment. Briefly, NanoString is a medium-throughput method that can analyze many genes within a single sample with comparable sensitivity and accuracy to quantitative real-time RT-PCR [49]. NanoString designed and manufactured custom probes corresponding to the 33 genes we identified for quantitative analysis, consisting of 30 target genes and 3 housekeeping genes (Gusb, Pgk1, and Eif4a2). A code set specific to a 100-base region of the target mRNA was designed using a 3' biotinylated capture probe and a 5' reporter probe tagged with a specific fluorescent barcode. Data were collected using the nCounter Digital Analyzer by counting the number of individual barcodes.
Each transcript of interest was recognized by a capture probe and a reporter probe, each containing 30-50 bases complementary to the target mRNA. To minimize assay variability, the code sets also included negative and positive control reporter probes that were developed by the External RNA Control Consortium (ERCC). Six positive control reporter probes (ERCC-selected mRNA targets) were pre-mixed with (Spike-Ins) the code set at a concentration range (0.125-128 fM), a range corresponding to the expression levels of most mRNA of interest, to control for overall efficiency of probe hybridization and determine the detection range for transcripts of interest in each assay. A scaling factor was calculated for each sample, and a scaling factor outside the range of 0.3-3 indicated suboptimal hybridization. In our samples, the scaling factor always fell within the optimal range and was thus applied to all counts in the sample.
Quantitative expression data from the nCounter were downloaded and analyzed using the nSolver software package (NanoString Technologies). The raw counts for all transcripts were multiplied by the scaling factor to produce the adjusted counts. The relative expression was determined for each comparison group, and the effect size of the difference between expression values was determined using Cohen's d. Expression was also compared using t-tests, and P-values were adjusted for multiple comparisons in nSolver.

GO Analysis
Individual genes are associated with GO annotations in order to describe the various functions of a particular gene product. The cellular component analysis describes the locations of gene expression, at the levels of subcellular structures. The molecular function analysis describes the function that each gene product performs within the cell. The biological process analysis describes a recognized series of events or collection of molecular functions associated with a gene or gene product. Each analysis was completed for all genes with differential patterns of expression between the three comparison groups, V vs P, V vs F, and P vs F. Because each GO annotation references many genes, in some instances the same GO annotation was present in multiple comparison groups.
The initial GO enrichment analysis returned 209 GO annotations in the V vs P comparison, 222 annotations in the V vs F comparison, and 264 annotations in the P vs F comparison that were significantly enriched. Upon selecting the GO annotations in each comparison that were related to the brain or behavior, we were left with 47 GO annotations in the V vs P comparison, 47 annotations in the V vs F comparison, and 61 annotations in the P vs F comparison (Tables 2-4). We then categorized these annotations based on gross function (Fig. 2).
The functional categories of GO annotations were differentially distributed across the three comparison groups. Annotations related to Neuropeptide activity were only found in the V vs P comparison, whereas immune function annotations were most predominant in the V vs F comparison. The P vs F comparison contained the greatest number of annotations related to Plasticity, DNA/RNA/Transcription, and Axon/Dendrite/Synapse.

KEGG Pathways Analysis
In the V vs P comparison group, the commonly recurring pathways included: protein export, protein processing in the endoplasmic reticulum, thyroid hormone synthesis, antigen processing and presentation, Ras signaling, Rap1 signaling, neuroactive ligand-receptor pathway, calcium signaling, and regulation of the actin cytoskeleton. In the V vs F comparison group, the commonly recurring pathways included: protein processing in the endoplasmic reticulum, regulation of the actin cytoskeleton, Ras signaling, metabolic pathways, axon guidance, protein processing in the endoplasmic reticulum, thyroid hormone synthesis, and antigen processing and presentation. In the P vs F comparison group, the commonly recurring pathways included: Ras signaling, Rap1 signaling, neuroactive ligand-receptor pathway, calcium signaling, MAPK signaling, LTP, glutamatergic pathways, dopaminergic pathways, and regulation of the actin cytoskeleton.
Using the list of genes generated from the GO annotations analysis, we next identified genes that were associated with multiple KEGG pathways. By excluding genes that were not associated with any KEGG pathways, or were associated with pathways that were not related to brain function, we further narrowed the range of genes of interest. Ultimately, in each comparison group we identified genes with patterns of expression that differed across social experience and that were linked to biological pathways within the brain (Table 5). We standardized the expression of each gene relative to its expression in virgin males then grouped genes that were associated with nine commonly recurring KEGG pathways and compared the expression of those genes across groups. Since this was an exploratory study, we did not perform statistical tests, and instead used Cohen's d as a measure of effect size (Table 6). Figure 3 shows the changes in gene expression in paired males (left) and fathers (right) compared with virgin males (dashed line). Overall, we observed only moderate changes in gene expression in paired males, but fathers exhibited an overall decrease in gene expression, especially in genes that were associated with long-term potentiation (LTP) and long-term depression (LTD) (d ¼ 1.072), neurotransmitters (d ¼ 0.911), and Ca 2þ signaling (d ¼ 0.877). We found medium effects of differential patterns of expression in genes that were associated with oxytocin signaling (d ¼ 0.787), protein processing in the endoplasmic reticulum (d ¼ 0.599), and Ras/Rap1 signaling (0.578). These results suggest that the genes associated with these KEGG pathways undergo coordinated changes in expression patterns that are related to social experience. Furthermore, different social experiences can result in dramatically different patterns of gene expression, i.e. Ca 2þ signaling.

STRING Database Analysis
We used the STRING database to assess the network connectivity between the genes in each comparison group that were identified as having both differential patterns of expression and functional significance.
Fatherhood alters gene expression within the MPOA | 5 The 11 genes from the V vs P comparison group produced a network with 11 nodes and 11 edges, and a protein-protein interaction (PPI) enrichment P-value of 5.86 Â 10 À7 (Fig. 4A). Thus, the proteins expressed by these genes have significantly more interactions than would be expected by chance, as defined as a random set of similarly sized proteins selected from the genome. There was one cluster of seven interacting proteins, and the functions of these gene products were primarily related to functions of the endoplasmic reticulum, as well as the cellular response to stimulation.
The 33 genes from the V vs F comparison group produced a network with 32 nodes and 29 edges, and a PPI enrichment P-value of 6.99 Â 10 À15 (Fig. 4B), indicating that the proteins expressed by these genes have significantly more interactions than would be expected by chance. These gene products produced one large cluster of nine interacting proteins, one medium cluster of five interacting proteins, and three separate small clusters of two interacting proteins. The large cluster was predominantly involved with the function of the endoplasmic reticulum. The medium cluster was involved with process of neural plasticity, including signaling pathways and modification of the actin cytoskeleton. The three small clusters were involved with the elongation of fatty acid chains, the formation of cholinergic receptors, and GPI-anchor synthesis.
The 31 genes from the P vs F comparison group produced a network with 31 nodes and 27 edges, and a PPI enrichment P-value of 1.36 Â 10 À12 (Fig. 4C), indicating that the proteins expressed by these genes have significantly more interactions than would be expected by chance. These gene products produced 1 large network consisting of 20 interacting proteins. The genes in this network were involved in a variety of functions, including synaptic plasticity and neural transmission, ion transmembrane transport, the cellular response to stimulus, and the structure of the synapse and dendrite.

NanoString Analysis
A total of 33 genes (30 target genes and 3 housekeeping genes) were selected for quantitative analysis using NanoString. The housekeeping genes (Gusb, Pgk1, and Eif4a2) did not show different levels of expression across conditions, confirming that these genes can serve as a good baseline in prairie voles. A heat map analysis revealed that 23 of our 30 target genes had lower expression levels in fathers than in either virgins or paired males (Fig. 5). Six genes had lower expression levels in virgins, and no gene in any group appeared to show inordinately high levels of expression. A regression analysis revealed similar levels of gene expression across all experimental conditions (Fig. 6A). Expression data for each individual gene were compared across groups using t-tests, which were run and P-values were adjusted for multiple comparisons using nSolver software. Of the 30 target genes, 11 genes showed significant differential expression between comparison groups (P < 0.05; Cckbr, Rgs14, Itpr1, Ddn, Baiap, Gabrd, Chrm1, Kcnj4, Ngef, Prkcg, and Cacna2d3;    Table 7). We also calculated the effect sizes using Cohen's d, examining differential expression of each gene across groups (Fig. 6B-O; Table 7). In the V vs P group, we saw a large effect (defined as 0.8 < d < 1.2) in Tiam1. In the P vs F

Discussion
In this experiment, we compared gene expression in the MPOA of virgin male prairie voles, males that had formed a pair bond, and males with fathering experience. We found that these groups differed in gene expression. Distinct patterns were revealed using a series of analyses, including GO annotation enrichment, KEGG pathways, STRING network analysis, and   quantitative assessment using NanoString. Males with fathering experience showed a relative decrease in gene expression compared with virgins and paired males, and many of the genes that exhibited decreased expression were involved in synaptic transmission and plasticity. These results suggest that fathers may exhibit a decreased amount of synaptic plasticity within the MPOA.
The transition to fatherhood is associated with a variety of potential environmental and behavioral changes, such as the   presence of infants, changes in energetic requirements and feeding behavior, and stress responsiveness [50][51][52]. In this experiment, we sought to identify alterations in central nervous system gene expression that are associated with fathering experience. To our knowledge, this is the first time that RNA sequencing has been performed on prairie vole brain tissue. As such, we faced several technical challenges over the course of this study. For instance, while the prairie vole genome has been sequenced [53,54], its annotation is incomplete, leaving us to rely on the annotated mouse genome (Mus musculus) for many of our analyses. In addition, it is important to consider the consequences associated with working in an outbred rodent-like prairie voles. The individual differences associated with an outbred population may have masked additional target genes associated with the onset of paternity. Furthermore, in this study we concentrated on traditional analyses. By examining protein coding transcripts, and restricting our analysis to pathways and genes that were related to the brain and behavior, it is possible that we overlooked some nontraditional candidate genes. Regardless, we still observed significantly altered expression on both the individual gene and system level. These results suggest that paternity engages similar physiological mechanisms across prairie vole males despite genetic diversity. Biparental care is rare in mammals, but prairie voles are not the only rodents who exhibit this behavior. The males of several species of Peromyscus, including Peromyscus californicus and Peromyscus polionotus, exhibit paternal care, while other species, including P. maniculatus, do not. This behavioral distinction allowed Bendesky et al. to investigate genetic differences between P. polionotus and P. maniculatus that are linked to parenting behavior [55]. In a series of experiments, they identified several quantitative trait loci that were linked to specific behaviors of interest, including nest building. Further analysis revealed that the gene for arginine vasopression (AVP) was directly related to nest building, and when AVP was administered intracerebroventricularly there was a significant decrease in the quality of nest building [55]. Unlike, the study by Bendesky et al., we did not find changes implicating AVP. However, there are several differences between the two experiments. In this study, we specifically examined gene expression within one hypothalamic nucleus, the MPOA. Our study was in a different species and used males that had very specific social experiences: virgin males, pair bonded males, and males with fathering experience.
RNA sequencing is a powerful technique that allows us to identify alterations in gene expression that are associated with behavioral and other phenotypic changes [56]. The greatest challenge with this technique, however, is the large amount of data it produces. There is no one agreed upon analysis that most effectively identifies specific genes of interest [57,58]. Thus, in this study we used several techniques to reveal novel gene targets to further our understanding of paternal behavior. We believe that this is a strength rather than a weakness. The ultimate goal of this experiment was to increase our understanding of the alterations that occur within the MPOA following exposure to different social contexts in male prairie voles. As such, we have identified a set of genes and their associated pathways that we can use to further explore male parenting behavior. Our quantitative assessment of gene expression in the MPOA revealed an overall decrease in the expression of many genes in fathers relative to both virgins and pair-bonded males. The specific genes of interest that we identified were involved in a range of physiological processes, including metabolism, stress responsiveness, and plasticity. However, most of the genes that showed different patterns of expression between groups, and specifically decreased expression in fathers, were associated with synaptic transmission and dendritic spine motility (Table 8). For example, several genes involved in the production and maintenance of receptors (including Cckbr, Chrm1, Gabrd, Grin2b, and Itpr1) and ion channels (including Cacna2d3, Kcnj4, and P2rx3) were significantly downregulated. These results suggest that GABA, glutamate, and cholinergic systems are all affected by fathering experience, as are calcium and potassium channels. Other genes that exhibited significant downregulation in fathers were involved with the actin cytoskeleton, dendritic spine motility, and other components of the physical plasticity of dendrites. We emphasize that this is not an exhaustive list of differentially expressed genes; however, these results suggest that synaptic plasticity may be diminished in the MPOA of male prairie voles with fathering experience.
We were surprised by the lack of differential expression of oxytocin and vasopressin-related genes; however, this finding is not unique within the literature. In a series of experiments, Kenkel et al. examined the neuroendocrine correlates of pup exposure in male prairie voles that were virgins or had fathering experience [52,59]. They saw changes in OT immunoreactivity in PVN/BNST, but there were no changes to OT/AVP in the MPOA. Another study examined OT immunoreactive cells in male prairie voles that were virgins, had established pair bonds, or had fathering experience [60]. They saw an increase in the number of OT immunoreactive cells in the MPOA of paired males and fathers compared with virgin males, but there was a In future studies we hope to examine patterns of gene expression in additional brain regions. Fatherhood also seems to be associated with structural alterations in neural plasticity, as measured by changes in the number and density of dendritic spines. Mice with fathering experience show increased survival of newborn neurons and increased dendritic spine density within the hippocampus [61,62]. Male marmosets show an increase in dendritic spine density in the prefrontal cortex after fathering experience [63]. However, other studies have shown reductions in the survival of adult-generated neurons in the amygdala of the prairie vole and hippocampus of California mice [64,65]. The effects of fatherhood clearly vary across brain regions, but we do not yet know what is causing these changes in neural plasticity.
The lower gene expression related to dendritic spines, associated with fatherhood in the present study, is evocative of similar changes seen in a recent study of the MPOA of mother rats [66]. Rem2, a gene associated with reduction of dendritic branching but increases in spine density [67,68] was increased in the MPOA of high licking/grooming rats, but only in lactating mothers (not in virgins). This increase was accompanied by decreased dendritic complexity. Rem2 is involved with GTPase activity and GTP binding. While we did not see alterations in Rem2 expression in this study, we found altered expression in several genes that are involved in Ras and Rap1 signaling. Both Ras and Rap1 are GTPases that work in concert to modulate cellular growth and plasticity [69][70][71]. Ras relays NMDA receptor signaling that drives the delivery of AMPA receptors during LTP, while Rap1 is involved in the NMDA receptor-dependent removal of AMPA receptors during LTD [69,72]. The altered pattern of expression of these genes in fathers suggests that this extremely salient social experience triggers a molecular cascade that is involved in neuronal plasticity.
The down-regulation of genes associated with dendritic complexity in the present study, as well as the study by Parent et al., is similar to what one would expect in an animal that had experienced high amounts of stress. It is well established that stress, mediated by corticotropin-releasing hormone, results in a loss of dendritic spines [73][74][75][76][77]. Additionally, rat mothers show a decrease in the number and density of dendritic spines in the amygdala and stria terminalis 4 days after birth [78], and an increase in dendritic spine density in the hippocampus during the postpartum period [79]. This suggests that alterations in dendritic spine density in mothers may be both transient and region specific, perhaps linked to the peak in corticosterone that occurs during parturition [80]. More studies must be done to determine if the same holds true for vole fathers.
In many species, the transition to fatherhood is associated with a suite of behavioral and hormonal changes, including those indicative of stress. In California mice (P. californicus), fathers exhibit attenuated anxiety-like behavior $2 weeks after pups are born [61,62]. Human males show a peak in cortisol levels during the transition to fatherhood [81]. Prairie vole fathers show increased anxiety-like behavior, and chronic pup exposure (in this case, 20 min of unrelated pup exposure per day for 10 days) results in an increase in basal CORT levels [64]. In an open field test, fathers spent more time in corner squares, and in an elevated plus maze, fathers spent less time in open arms. In forced swim tests, fatherhood decreased the latency to immobility, and increased the number and duration of immobility bouts [64]. In the long-term, fatherhood may be beneficial for male health, but the transition to fatherhood is a tremendously stressful period [82].
In male voles with fathering experience, we also see the upregulation of genes related to protein processing in the endoplasmic reticulum. The endoplasmic reticulum is instrumental in managing the protein folding process, including disposing of misfolded proteins [83]. Homeostatic imbalances, including stress, can alter the functioning of the endoplasmic reticulum, leading to the initiation of the unfolded-protein response, which can in turn lead to apoptosis [84,85]. This may be one mechanism by which physiological stress can result in homeostatic perturbations [86,87], including some of the changes that are evident in vole fathers, such as weight loss [50,59].
While many of the changes we saw in gene expression may be partially attributable to stress, there are likely many other additional factors at play. Fathers in many species show systematic endocrine changes [22]. Environmental factors, including changes in the types and amount of sensory stimulation, or the amount of parental care they received, may play a role as well [88,89]. Much more work must be done to tease apart these many factors.
In this study, we saw the most varied and interesting differences between the paired males and males with fathering experience. This was surprising, as we expected that the greatest differences would be between the virgin males and fathers. However, examination of the quantitative results begins to clarify these findings (Figs 3 and 6). The expression of genes of interest is slightly elevated in paired animals relative to virgins, but the expression in fathers is decreased relative to virgins. Thus, while the expression levels of some genes do not significantly differ between virgins and paired males, and virgins and fathers, we found significant differences between paired males and fathers. This may suggest that the experience of fathering is functionally distinct from any other type of social interactions that these animals have encountered.

Conclusions
The purpose of this study was to explore how gene expression changed across the transition to fatherhood, and to identify novel targets to allow for deeper investigation of male parenting behavior. The use of RNA sequencing confirmed that there are differences in gene expression between voles that had different social experiences, including virgin males, males that had formed a pair bond with a female, and males with parenting experience. The genes identified in this study suggest novel processes that are related to paternal behavior and offer new targets for the further exploration of fathering behavior.