Skip to main content
Download PDF

Genetic, epigenetic, and environmental factors controlling oxytocin receptor gene expression

Clinical Epigenetics volume 13, Article number: 23 (2021) Cite this article

Abstract

Background

The neuropeptide oxytocin regulates mammalian social behavior. Disruptions in oxytocin signaling are a feature of many psychopathologies. One commonly studied biomarker for oxytocin involvement in psychiatric diseases is DNA methylation at the oxytocin receptor gene (OXTR). Such studies focus on DNA methylation in two regions of OXTR, exon 3 and a region termed MT2 which overlaps exon 1 and intron 1. However, the relative contribution of exon 3 and MT2 in regulating OXTR gene expression in the brain is currently unknown.

Results

Here, we use the prairie vole as a translational animal model to investigate genetic, epigenetic, and environmental factors affecting Oxtr gene expression in a region of the brain that has been shown to drive Oxtr related behavior in the vole, the nucleus accumbens. We show that the genetic structure of Oxtr in prairie voles resembles human OXTR. We then studied the effects of early life experience on DNA methylation in two regions of a CpG island surrounding the Oxtr promoter: MT2 and exon 3. We show that early nurture in the form of parental care results in DNA hypomethylation of Oxtr in both MT2 and exon 3, but only DNA methylation in MT2 is associated with Oxtr gene expression. Network analyses indicate that CpG sites in the 3′ portion of MT2 are most highly associated with Oxtr gene expression. We also identify two novel SNPs in exon 3 of Oxtr in prairie voles and a novel alternative transcript originating from the third intron of the gene. Expression of the novel alternative transcript is associated with genotype at SNP KLW2.

Conclusions

These results identify putative regulatory features of Oxtr in prairie voles which inform future studies examining OXTR in human social behaviors and disorders. These studies indicate that in prairie voles, DNA methylation in MT2, particularly in the 3′ portion, is more predictive of Oxtr gene expression than DNA methylation in exon 3. Similarly, in human temporal cortex, we find that DNA methylation in the 3′ portion of MT2 is associated with OXTR expression. Together, these results suggest that among the CpG sites studied, DNA methylation of MT2 may be the most reliable indicator of OXTR gene expression. We also identify novel features of prairie vole Oxtr, including SNPs and an alternative transcript, which further develop the prairie vole as a translational model for studies of OXTR.

Background

A major goal of translational neuroscience is to identify biomarkers of psychiatric disorders which can be used to inform diagnosis, predict possible courses of disease progression, or predict likelihood of response to a particular treatment [1]. A common target of efforts to identify biomarkers of psychiatric disorders is the oxytocin system [2]. The oxytocin system is an attractive target for biomarkers because of its involvement in many psychiatric disorders [3], modulation of a wide range of physiological processes [4], and potential use as medication via intranasal administration [5,6,7]. While many studies use serum oxytocin levels as a biomarker, proper collection and measurement of oxytocin is challenging, and interpretation may vary based on methodology [8].

Epigenetic modifications are promising biomarkers for disease since they are easily measured, variable in the population, and can be influenced by the environment [9]. A prominent epigenetic biomarker is DNA methylation, which is the addition of a methyl group to cytosine residues. DNA methylation typically occurs when cytosine is followed by guanine and generally leads to decreased gene expression [10]. One epigenetic biomarker for the oxytocin system is DNA methylation of the oxytocin receptor gene, OXTR. OXTR structurally consists of four exons and three introns in humans [11]. DNA methylation in MT2, a 405 bp region overlapping exon 1 and intron 1 of OXTR, regulates transcription both in vivo and in vitro [12]. OXTR is expressed throughout the human brain, with highest levels of expression in striatum, thalamus, and olfactory regions [13]. Reduced expression of OXTR has been documented in the superior temporal gyrus in autism [14] and in the posterior medial temporal cortex in schizophrenia [15]. Increased expression of OXTR was reported in the prefrontal cortex of patients with depression and bipolar disorder [16].

Dysregulation of OXTR expression in multiple psychopathologies has led many to study the epigenetic state of OXTR in both typical and clinical groups. These studies have generally focused on two regions of the OXTR gene: MT2 (covering much of exon 1 and some of intron 1) and exon 3. For example, dysregulated DNA methylation in MT2 has been associated with autism [14], callous–unemotional traits [17, 18], anorexia nervosa [19, 20], schizoaffective disorders [21, 22], post-partum depression [23], attachment anxiety [24], obsessive compulsive disorder [25], and anxiety and depression [26]. DNA methylation in MT2 is associated with neural endophenotypes in typical populations: DNA methylation is positively correlated with brain activity during animacy perception [27], emotion perception [28, 29], and social attention [30]. DNA methylation in exon 3 has also been associated with several psychopathologies, including anxiety and depression [26, 31], social anxiety disorder [32], obsessive compulsive disorder [25, 33], post-traumatic stress disorder [34], callous–unemotional traits [17], and social communication deficits [35]. Despite the many studies connecting epigenetic markers at OXTR and psychological outcomes, several questions remain. Most studies measure DNA methylation in peripheral tissues such as blood or saliva, and the association between DNA methylation in these tissues and the brain is not established for all CpG sites. Furthermore, even if DNA methylation is measured in the brain, it is not yet known what the relative contributions of DNA methylation in the MT2 region and exon 3 region are in controlling OXTR expression.

In order to investigate epigenetic mechanisms controlling OXTR expression in the brain we use prairie voles as an animal model. Prairie voles are highly social rodents with characteristics resembling human social behavior including pair bonding and biparental care of offspring [36]. We have previously shown that prairie voles are also a good model for studies of OXTR: the region of MT2 containing CpG sites -901, -924, and -934 is highly conserved in prairie voles and humans, but not in mice or rats [37]. We have also shown that DNA methylation in this region is malleable to early life experience, which ultimately leads to changes in gene expression in the juvenile period [37] and adulthood [38]. In this study, we clarify the genetic structure of Oxtr in prairie voles and compare the gene structure to human OXTR. We then investigate the relationships between DNA methylation in MT2, DNA methylation in exon 3, and Oxtr expression in the prairie vole brain. Finally, we describe novel findings related to prairie vole Oxtr including previously unidentified SNPs in exon 3 and a novel alternative transcript present in the nucleus accumbens.

Results

Oxtr gene structure in prairie voles resembles human OXTR gene structure with high homology in CpG-rich regulatory regions

Previously, the structure of the Oxtr gene in prairie voles was reported to have two exons which are homologous to exons 3 and 4 in humans (MicOch1.0, GCA_000317375.1, Ensembl release 100). However, this gene in other rodents has four exons, similar to humans, including two exons at the 5′ end of the gene which are not translated [39]. In order to determine the full structure of Oxtr in prairie voles, we completed RNA-sequencing on RNA from the nucleus accumbens. Viewing the read alignment in the Integrated Genomes Viewer (IGV) reveals that this gene has four exons and three introns, much like the human gene (Fig. 1a, b). In addition, the 3′ UTR extends further than previously reported.

Fig. 1
figure 1

Oxtr gene structure and homology in prairie voles and humans. a RNA-seq evidence of Oxtr gene structure in prairie voles. The region shown is micOch1.0 scaffold JH996431.1: 26,340,626–26,358,364. b Gene schematic of OXTR in humans and Oxtr in prairie voles. Exons are indicated by boxes (coding, white; untranslated, gray). Introns are shown with solid lines. The black boxes below the transcripts delineate a CpG island identified using the CpG island function of the sequence manipulation suite [40]. Gray boxes above the transcripts denote conserved regions of the CpG island labeled as MT2 and exon 3. These conserved regions were assayed for DNA methylation

Full size image

We assessed the homology between prairie voles and human for two CpG-rich regions in OXTR where many human studies associated DNA methylation and neuropsychological outcomes, MT2 and exon 3. Homology was determined using PRSS (DNA:DNA). For the MT2 region, this analysis identified a conserved region between human and prairie voles with 64.3% shared identity (Additional file 1: Fig. S1A; Smith-Waterman score: 566; 64.3% similar; Z-score: 120.2; bits: 120.2; E(1000): 1 × 10–31). For the exon 3 region, this analysis identified a conserved region between human and prairie voles with 88.4% identity (Additional file 1: Fig. S1B; Smith-Waterman score: 2275; 88.4% similar; Z-score: 1966.2; bits: 373.7; E(1000): 1 × 10–107).

DNA methylation in MT2 is sensitive to early life experience and associated with Oxtr gene expression in prairie vole and human brains

We have previously shown that early life experience in the form of nurture prevents de novo DNA methylation on four sites in MT2 (-934_1, -934_2, -924, and -901) and leads to increased Oxtr gene expression as juveniles [37]. We next aimed to characterize how early life experience affects DNA methylation on other sites in MT2. We used a mixed model analysis to account for correlation between DNA methylation values at these sites (Additional file 1: Fig. S2A). Across all sites in MT2, offspring that received less parental care (MAN0) have increased DNA methylation (Fig. 2a, type II Wald F test F(1,24) = 8.386, p = 0.008). Several of the sites are significantly correlated with expression of Oxtr (Fig. 2b, Additional file 1: Table S1). In order to determine if this model is instructive for studying humans, we compared OXTR DNA methylation to gene expression from human samples collected from Brodmann Area 41/42 [14]. In human brain tissue, there is also a significant correlation between DNA methylation at CpG sites -934, -924, and -901 and OXTR gene expression (Fig. 2c, Additional file 1: Table S2, -934: ρ = − 0.78, p = 0.007; − 924: ρ = − 0.68, p = 0.02; − 901: ρ = − 0.66, p = 0.028).

Fig. 2
figure 2

DNA methylation in MT2 is impacted by early life experience and associated with Oxtr expression. a Offspring that received less parental care (MAN0) have increased DNA methylation across MT2 in the nucleus accumbens (13 male–female sibling pairs; MAN0: n = 7 males, 7 females; MAN1: n = 6 males, 6 females; mixed effect model with fixed effect of handling condition and random effects of individual and CpG site; Type II Wald F test F(1,24) = 8.386, p = 0.008). Follow-up t tests between handling conditions at each site were conducted with Bonferroni correction; *p < 0.05, **p < 0.01. b DNA methylation in MT2 is associated with Oxtr expression in prairie vole nucleus accumbens. Spearman’s rank correlation was used to determine correlation. Sites with significant correlations (not corrected for multiple comparisons) are shown in color while nonsignificant sites are in gray. Spearman’s rho and p values for each site are in Additional file 1: Table S1. c DNA methylation in MT2 is associated with OXTR expression in human temporal cortex. Spearman’s rank correlation was used to determine correlation. Sites with significant correlations (not corrected for multiple comparisons) are shown in color while nonsignificant sites are in gray. Spearman’s rho and p values for each site are in Additional file 1: Table S2

Full size image

CpG sites -934_1, -934_2, -924, and -901 are most highly associated with Oxtr gene expression

DNA methylation at CpG sites in MT2 are highly correlated with each other (Additional file 1: Fig. S2A). In order to further understand the structure of MT2, we employed exploratory graph analysis (EGA), a dimensionality reduction technique and network structure analysis tool [41, 42]. This tool allows us to cluster CpG sites in an unsupervised manner based on the covariance of DNA methylation values into distinct groups which provide non-redundant information. The results from EGA suggest three distinct communities of CpG sites, largely reflecting the physical structure of MT2: community 1 contains 5′ sites, community 2 contains sites in the middle, and community 3 contains 3′ sites (Fig. 3a, Additional file 1: Fig. S2B-C). The reliability of the community structure was determined using 1000 bootstrapped replications [43]. The probability of each CpG site remaining in their community is shown in Fig. 3b. Probabilities greater than 0.8 are considered stable. Notably, CpG sites -934_1, -934_2, -924, and -901 all remain in their community nearly 100% of the time, and this stability may indicate that these CpG sites are most tightly co-regulated.

Fig. 3
figure 3

Network analysis of MT2 DNA methylation reveals importance of 3′ region. a Community structure network of CpG sites in MT2 containing three communities: community 1 (green) which contains 5′ CpGs, community 2 (orange) which contains mostly CpGs in the middle of MT2, and community 3 (blue) which contains mostly 3′ CpGs. The weight of edges indicates higher covariance between CpG sites. b Item stability scores for each CpG site. CpGs with scores greater than 0.8 are considered stable. c Animals raised by low care parents (MAN0) have higher standardized network community scores, which represent DNA methylation within the community (13 male–female sibling pairs; MAN0: n = 7 males, 7 females; MAN1: n = 6 males, 6 females; 2 way handling condition x community ANOVA, main effect of handling condition, F(1, 72) = 22.41, p = 1.1e-5; no effect of community, F(2, 72) = 0.001, p = .99. Follow-up t tests between handling conditions at each community were conducted with Bonferroni correction (*p < 0.05, **p < 0.01). Cohen’s d is provided to measure effect size.) The effect size of handling condition on standardized network community score is largest in community 3. d Standardized network community scores for each community are correlated with Oxtr gene expression using Spearman’s rank correlation; community 1: ρ = -0.29, p = 0.086; community 1: ρ = -0.37, p = 0.03; community 3: ρ = -0.47, p = 0.003

Full size image

In order to determine which community of CpG sites are most sensitive to early nurture, we compared standardized network community scores [44] for each community between the two handling groups. The standardized network community scores are compressed representations of the original DNA methylation data from CpG sites within the community and are analogous to component vectors from principal components analysis. In this way, the network scores represent weighted DNA methylation values from CpG sites within each network community. The standardized network community scores differed by handling condition but not by community (Fig. 3c, two-way community x handling condition ANOVA, main effect of handling condition F(1, 72) = 22.41, p = 1.1e-5, no effect of community F(2, 72) = 0.001, p = 0.99). We then performed follow-up t tests to determine which community was most affected by early nurture. The largest effect was found in community 3, which contains sites -934_1, -934_2, -924, and -901 (Welch’s two sample t test with Bonferroni correction and Cohen’s d; community 1: t(17.3) = 1.88, p = 0.229, d = 0.77; community 2: t(22.7) = 2.72, p = 0.037, d = 1.08; community 3: t(23.4) = 3.70, p = 0.003, d = 1.42), indicating that DNA methylation at this cluster of sites is particularly sensitive to early nurture.

We then correlated standardized network community score with Oxtr expression. Community 3 most correlates with gene expression (Fig. 3d; community 1: ρ = -0.29, p = 0.086; community 1: ρ = -0.37, p = 0.03; community 3: ρ = -0.47, p = 0.003). Notably, community 3 contains five of the eight sites in MT2 where DNA methylation significantly correlates with Oxtr expression (Fig. 2b). This is also the community where DNA methylation is most sensitive to early life experience, indicating that CpG sites in this community may be functionally responsible for regulation of gene expression.

DNA methylation in exon 3 is sensitive to early life experience but is not associated with gene expression

In human studies, DNA methylation at CpG sites within exon 3 of Oxtr have been associated with several neuropsychological outcomes [17, 26, 31,32,33,34,35]. In prairie voles there is a highly conserved region within exon 3 that contains 42 CpG sites (Fig. 1b). We analyzed DNA methylation across this region and found that DNA methylation in exon 3 is sensitive to early life experience (Fig. 4a, fixed main effect of handling condition, type II Wald F test F(1,24) = 8.386). Follow-up t tests were completed at each site to determine the sites most sensitive to early life experience, but no sites remained significant after Bonferroni correction. DNA methylation in exon 3 was only associated with Oxtr gene expression at CpG 36 (Fig. 4b, Additional file 1: Table S3). Thus, while early nurture affects DNA methylation in both MT2 and exon 3 of Oxtr, DNA methylation levels in MT2 appear to be an important indicator of expression of the gene but the same is not true of exon 3. Studies in humans identifying associations between DNA methylation and neuropsychological outcomes likely find such results because DNA methylation in exon 3 is highly correlated with DNA methylation in the 3′ portion of MT2 (Fig. 5).

Fig. 4
figure 4

Genetic and epigenetic features of exon 3 do not control Oxtr expression. a Offspring receiving less parental care (MAN0) have decreased DNA methylation across exon 3 in the nucleus accumbens (13 male–female sibling pairs; MAN0: n = 7 males, 7 females; MAN1: n = 6 males, 6 females; mixed effect model with fixed effect of handling condition and random effects of individual and CpG site; Type II Wald F test F(1,24) = 8.386, p = 0.018. CpG sites are numbered according to Additional file 1: Fig S1B. b DNA methylation in exon 3 is not robustly associated with Oxtr expression in prairie vole nucleus accumbens. Spearman’s rank correlation was used to determine correlation. CpG 36, the only site that has a significant correlation with Oxtr expression (not corrected for multiple comparisons) is shown in black. Sites that do not have a significant correlation with Oxtr expression are in gray. Spearman’s rho and p values for each site are in Additional file 1: Table S3

Full size image
Fig. 5
figure 5

Correlation of DNA methylation in MT2 CpG sites and exon 3 CpG sites. Pearson’s R correlation of DNA methylation in MT2 and exon 3 CpG sites. The red box around MT2 CpGs 18 through -901 indicates that these CpGs have the highest correlation with exon 3 sites

Full size image

Identification of two single nucleotide polymorphisms (SNPs) within exon 3 which are not associated with Oxtr gene expression

Within exon 3 there are two CpG sites with patterns of DNA methylation with a cluster of prairie voles with DNA methylation between 50–75%, a cluster between 25–30%, and a few individuals near 0% methylation. We hypothesized that these two sites are polymorphic with one allele preserving the CpG and the other disrupting the CpG such that the cytosine is not methylated. In the same cohort of MAN0 and MAN1 animals, pyrosequencing analysis revealed that these two sites, KLW1 (CpG 25, JH996431.1:26,356,275, micOch1.0) and KLW2 (CpG 33, JH996431.1:26,356,107, micOch1.0), are indeed single nucleotide polymorphisms (SNPs). KLW1 is an A/G SNP where A is the minor allele in this population (frequency = 0.33). KLW2 is also an A/G SNP where A is the minor allele in this population (frequency = 0.27). Both SNPs are in Hardy–Weinberg equilibrium in this population (Haldane’s exact test [45]; KLW1: D = -0.22, p = 0.84; KLW2: D = -0.11, p = 0.82). We also tested if these SNPs are in linkage disequilibrium with each other or with another SNP in this gene, NT213739 [46]. All three SNPs are in linkage disequilibrium (Table 1).

Table 1 Linkage of SNPs KLW1, KLW2, and NT213739
Full size table

However, a larger sample is necessary for more accurate quantitation of linkage disequilibrium. We next examined the effect of genotype at KLW1 and KLW2 on expression of Oxtr. Neither KLW1 genotype nor KLW2 genotype impacted Oxtr expression (Fig. 6a, b; KLW1: t(23.9) = -0.06, p = 0.95; KLW2: t(20.4) = -0.89, p = 0.38).

Fig. 6
figure 6

Genetic x epigenetic interactions affecting Oxtr expression. a Genotype at KLW1 does not impact Oxtr expression (A carriers: n = 14; G/G: n = 12; Welch two sample t test, t23.9 = -0.06, p = 0.95). Animals are grouped by the presence of minor allele (A) because of low abundance of A/A homozygotes. The DNA sequence surrounding KLW1 is shown in the below graph. b Genotype at KLW2 does not impact Oxtr expression (A carriers: n = 12; G/G: n = 14; Welch two sample t test, t20.4 = -0.89, p = 0.38). Animals are grouped by presence of minor allele (A) because of low abundance of A/A homozygotes. The DNA sequence surrounding KLW2 is shown in the below graph. c The relationship between MT2 DNA methylation (represented by community 3 standardized network score) and Oxtr expression does not differ between KLW1 genotypes (A carriers: n = 14, R = -0.49, p = 0.076; G/G: n = 12, R = -0.51, p = 0.091), but D) does differ between KLW2 genotypes (A carriers: n = 12, R = -0.84, p = 0.00065; G/G: n = 14, R = -0.35, p = 0.22) and E) NT213739 genotypes (C/T: n = 9, R = -0.17, p = 0.66; T/T: n = 17, R = -0.61, p = 0.01)

Full size image

An interaction between DNA methylation of OXTR and SNP rs53576 genotype has been reported in relation to neuropsychological outcomes in humans [23, 31, 35]. This led us to hypothesize that there are genetic x epigenetic interactions affecting Oxtr expression. In order to study this in the prairie vole sample, we asked if the relationship between the standardized network community score of community 3, which represents DNA methylation at the 3′ end of MT2, and Oxtr expression differs by genotype at SNPs KLW1, KLW2, or NT213739. There are no significant interactions, likely because of the small sample size, but interesting patterns emerge. The correlation between community 3 standardized network score and Oxtr expression is very similar in KLW1 A carriers and G/G voles (ρ = -0.56 and -0.48, respectively; Fig. 6c). ‘A’ carriers (A/A and A/G) of KLW2 have a stronger correlation between these measures compared to KLW2 G/G voles (ρ = -0.88 and -0.28, respectively; Fig. 6d). Voles with T/T genotype at NT213739 also have a stronger correlation between these two measures compared to C/T individuals (ρ = -0.67 and -0.20, respectively; Fig. 6e). These results suggest a genetic by epigenetic interaction may impact Oxtr expression, but further studies with larger samples from several prairie vole colonies are needed to test this hypothesis.

Identification of an alternative transcript of Oxtr beginning in intron 3 which is associated with KLW2 genotype

In addition to Oxtr gene structure, RNA-seq data revealed a possible alternative transcript in Oxtr consisting of a small exon in the third intron of the main transcript and the fourth exon of the main transcript (Fig. 7a, reads highlighted in yellow). We then measured this transcript in samples from the early nurture experiment using RT-PCR. The transcript is present in all samples, and there are no differences in expression based on nurture (Fig. 7b, W = 82, p = 0.85) or sex (Fig. 7c, W = 63, p = 0.44). Expression of the alternative transcript is associated with the KLW2 A allele (Fig. 7d, W = 156, p = 3.85e-7). Alternative transcript expression is not associated with either KLW1 genotype (Additional file 1: Fig. S3A, W = 77, p = 0.98) or NT213739 genotype (Additional file 1: Fig. S3B, W = 42, p = 0.10).

Fig. 7
figure 7

Identification of an alternative transcript of Oxtr that is associated with KLW2 genotype. a RNA-seq evidence of a novel alternative transcript of Oxtr originating in the third intron. Reads providing evidence of this transcript are highlighted in yellow. b Expression of the alternative transcript is not impacted by sex (males: n = 13, females: n = 13; Wilcoxon rank-sum test, W = 63, p = 0.44). c Expression of the alternative transcript is not impacted by handling condition (MAN0: n = 14, MAN1: n = 12; Wilcoxon rank-sum test, W = 82, p = 0.85). d Expression of the alternative transcript is associated with the A allele of KLW2 (A carrier: n = 12, G/G: n = 14; Wilcoxon rank sum test, W = 156, p = 3.85e-7). e Expression of the alternative transcript is positively correlated with expression of the main Oxtr transcript in animals carrying the A allele of KLW2 (2 way KLW2 genotype x Oxtr main transcript expression ANOVA, interaction effect F(1,21) = 23.723, p = 8.14e-5, adjusted R2 = 0.734; A carriers: Spearman’s ρ = 0.81, p = 0.0024; G/G: Spearman’s ρ = 0.48, p = 0.097)

Full size image

We next asked if the main transcript and this alternative transcript are co-regulated, meaning that animals with high expression of the main transcript could also have high levels of alternative transcript expression. A two-way ANOVA was conducted to determine the influence of main transcript expression while accounting for KLW2 genotype. There was a significant interaction between the expression of the main Oxtr transcript and KLW2 genotype (Fig. 7e, F(1,21) = 23.723, p = 8.14e-5) such that Oxtr main transcript expression is positively correlated with alternative transcript expression in KLW2 A carriers (ρ = 0.81, p = 0.0024) but there is no correlation in animals with KLW2 G/G genotype (ρ = 0.48, p = 0.097). Thus, it appears that the main transcript of Oxtr is co-regulated with the alternative transcript via a mechanism that is associated with the A allele of KLW2.

Discussion

The results from these studies clarify the contributions of genetic, epigenetic, and environmental factors in determining Oxtr expression in the prairie vole brain. We show that the structure of Oxtr in prairie voles resembles human OXTR. Within the CpG island spanning the first 3 exons of OXTR, two regions have been associated with human neuropsychological outcomes: MT2 and exon 3, both of which are highly conserved in the prairie vole. Using the prairie vole model, we show that Oxtr expression in the nucleus accumbens is significantly associated with DNA methylation in MT2 but not DNA methylation in exon 3. Furthermore, network analysis implemented in EGA indicates that DNA methylation in the 3′ end of MT2, which contains CpG sites -901, -924, -934_1, and -934_2 that we have previously reported on [37, 38], is most strongly associated with Oxtr expression. Notably, we previously showed that DNA methylation of these CpG sites in whole blood samples are correlated with gene expression in the brain in prairie voles [37]. Additionally, we show that CpG sites -901, -924, and -934 are strongly associated with OXTR expression in human temporal cortex while other sites within MT2 are not. Previous studies finding associations between DNA methylation in exon 3 and neuropsychological outcomes in humans likely find this because DNA methylation in exon 3 reports on DNA methylation at CpG sites in the 3′ portion of MT2. Collectively, these results indicate that human studies measuring OXTR DNA methylation should focus measurements on MT2, where DNA methylation in both central and peripheral tissues serve as a biomarker for gene expression in the brain.

Our findings expand upon our previous studies examining the relationship between early life experience and Oxtr DNA methylation. Here, we show that increased parental care results in hypomethylation across MT2 and exon 3 in the prairie vole brain. These results mirror human studies which find that DNA methylation of OXTR is impacted by early life experience in both MT2 [47, 48] and exon 3 [49]. Additionally, childhood abuse in humans interacts with OXTR DNA methylation in both MT2 and exon 3 to predict anxiety and depression [26]. Our results further establish prairie voles as an excellent rodent model for studying how early life experience impacts regulation of the oxytocin system and social behavior.

We identified two novel SNPs in prairie vole Oxtr within exon 3. Both are A/G SNPs where A is the minor allele and disrupts a CpG site. These SNPs, in addition to the previously identified SNP NT213739 [46], offer an opportunity to model the common association of OXTR SNPs with behavioral and psychiatric outcomes in humans [50]. Our results show that neither KLW1 nor KLW2 genotype is associated with Oxtr expression in the nucleus accumbens. It was previously reported that the C allele of NT213739 was associated with increased Oxtr expression in the prairie vole nucleus accumbens [46], although our group found no association [37]. However, we found that there is a stronger negative association between Oxtr DNA methylation and gene expression in animals with A/A and A/G genotypes at KLW2 compared to those with G/G genotypes. Similarly, we found that animals with T/T genotype at NT213739 show a strong negative association between Oxtr DNA methylation and gene expression while there was no association in animals with the C/T genotype. Although this analysis is underpowered, these results suggest that the KLW2 G allele and NT213739 C allele lead to epigenetic dysregulation of Oxtr, perhaps by impeding the binding of factors that regulate the gene. Further studies with appropriate sample size will be necessary to study this phenomenon. These studies may shed light on mechanisms underlying the common finding in human studies that OXTR DNA methylation interacts with genotype at OXTR SNP rs53576 in certain psychopathologies [23, 31, 35].

The RNA-sequencing results led to the identification of a novel transcript of Oxtr. This transcript, which is homologous to mouse Oxtr isoform H (GenBank KU686801.1) [39], originates in the third intron of Oxtr. Expression of this transcript is associated with the A allele of KLW2, though further work is necessary to determine if this SNP is functionally related to transcription of this alternative transcript. In animals with a KLW2 A allele, expression of the main transcript and the alternative transcript of Oxtr is positively correlated. The function of the novel transcript is unknown. One possibility is that the alternative transcript regulates expression of the main transcript. For example, it has been reported that non-coding RNAs bind DNMT1 to prevent DNA methylation in a gene-specific manner, which in turn increases expression of that gene [51]. If the novel transcript is translated, it would result in a peptide containing most of the seventh transmembrane helix and the C-terminus of OXTR. This peptide could alter the activity of the OXTR protein. Oxtr has been shown to dimerize with type 2 dopamine receptors, though the precise dimerization sites are unknown [52, 53]. It is possible that this OXTR fragment could dimerize with D2Rs, affecting the function of oxytocin–dopamine interactions in the nucleus accumbens. This novel OXTR fragment could also disrupt desensitization of OXTR after oxytocin stimulation, which is dependent on β-arrestin binding to serine triplets in the C-terminus of OXTR [54], by sequestering available β-arrestin. Further studies are necessary to determine the function of this alternative transcript. While this transcript is present in mouse brain, we could not find any evidence of this transcript in human nucleus accumbens using publicly available RNA-seq data from the GTEx project [55]. Further characterization of this transcript in humans will require appropriate, high-quality human brain samples.

Taken together, the results from these experiments offer insight into mechanisms of regulation of Oxtr expression in the brain. These results will inform future studies of OXTR involvement in human psychopathologies. Further work should focus on the 3′ region of MT2, where DNA methylation is most sensitive to early life experience and most strongly correlates with OXTR expression in both prairie vole and human brain. Additional studies in prairie voles focused on the SNPs KLW2 and NT213739 can elucidate genetic by epigenetic interactions in OXTR which mediate behavioral outcomes. The identification of a novel transcript of OXTR, if also present in humans, has implications for OXTR gene regulation, protein function, and behavioral outcomes.

Conclusions

These results provide strong evidence that DNA methylation in MT2, particularly the 3′ portion, is a better indicator of OXTR gene expression than DNA methylation in exon 3 in prairie voles. We also provide evidence that in human temporal cortex, DNA methylation in the 3′ portion of MT2 is associated with OXTR gene expression while DNA methylation in other MT2 CpG sites is not. Though associations exist between DNA methylation in OXTR exon 3 and neuropsychological outcomes in humans, such results may reflect that DNA methylation in exon 3 is correlated with DNA methylation in MT2. Our findings suggest that future studies of OXTR DNA methylation in humans should focus on MT2, particularly CpG sites -934, -924, and -901, as these sites are most related to gene expression in both prairie vole and human brains.

Methods

Human brain DNA and RNA samples

Genomic DNA and RNA isolated from the temporal cortex (BA 41/42) was obtained from the Maryland National Institute of Child Health and Human Development Brain Tissue Center and the Harvard Brain Tissue Resource Center (n = 11, male, 11-30yrs). Genomic DNA was analyzed for OXTR DNA methylation at CpG sites in the MT2 region of the promoter and gene expression was evaluated as previously described (Gregory et al., 2009).

Animal model

Subjects were laboratory-bred prairie voles (Microtus ochrogaster), descendants of a wild-caught stock captured near Champaign, Illinois. Breeding pairs were housed in large polycarbonate cages (44 cm x 22 cm × 16 cm) and same sex offspring pairs were housed in smaller polycarbonate cages (27 cm × 16 cm × 16 cm) after weaning on postnatal day (PND) 20 (date of birth = PND0). Animals were given food (high-fiber Purina rabbit chow) and water ad libitum, cotton nestlets for nesting material in breeding cages, and were maintained on a 14:10 light/dark cycle. All procedures involved in generating tissue for the analysis of DNA methylation and gene expression following handling were reviewed and approved by the IACUC at the University of California, Davis.

Handling manipulation

Within 24 h of giving birth, breeding pairs underwent a single handling manipulation as previously described [37, 56]. Following this manipulation, parental care is heightened in the MAN1 group compared to the MAN0 group [37]. On PND24 offspring were anesthetized with isoflurane and euthanized via cervical dislocation and rapid decapitation for tissue collection. Brains were extracted and flash frozen in liquid nitrogen and stored at -80 °C. One male and one female from each litter were included in the analysis of Oxtr DNA methylation and Oxtr expression (MAN0: n = 7 per sex; MAN1: n = 6 per sex). An additional five sibling pairs were included in the EGA analysis in order to achieve adequate sample size such that the correlation matrix was full rank. These animals did not undergo MAN0 or MAN1 manipulations.

Sectioning of the nucleus accumbens

Whole brains were stored at -80 °C and equilibrated to -20 °C for two hours prior to sectioning. Following sectioning, nucleus accumbens tissue was placed in a DNAse/RNAse free microcentrifuge tube and flashed frozen with liquid nitrogen. Brain tissue was crushed using a mortar and pestle in preparation for DNA/RNA isolation.

Identifying conserved MT2 and exonic regions in prairie voles

The structure of Oxtr in prairie voles was determined by viewing RNA-sequencing reads near Oxtr. The sequence was taken for the whole region where reads were mapped and used to determine conserved regions of the gene. The MT2 region of OXTR in humans as identified by Kusui et al., 2001 (GRCh38, chr3:8,769,033–8,769,438) was compared to the Oxtr sequence in prairie voles and with the NCBI BLASTn program, and a region in prairie voles was identified as having a highly similar sequence of DNA. To assess the significance of this similarity of the MT2 region between human and prairie vole we used the UVa FASTA server (http://fasta.bioch.virginia.edu/) and PRSS (DNA: DNA) to shuffle the prairie vole sequence 200 times and estimate the statistical significance of the shuffled scores.

Oxtr DNA methylation analysis

Extraction of DNA from nucleus accumbens tissue was done using the Qiagen AllPrep DNA/RNA Mini Kit (Qiagen, Valencia, CA) following manufacturer instructions. Two hundred nanograms of DNA was subject to bisulfite treatment (Kit MECOV50, Invitrogen, Carlsbad, CA), which allows for the detection of methylated cytosines by sequencing. Twelve nanograms of bisulfite converted DNA was used as template for PCR reactions. PCR was performed using a Pyromark PCR kit (Qiagen, Valencia, CA), and each PCR reaction was amplified in triplicate on three identical PCR machines (S1000 Thermal Cycler, Bio-Rad, Hercules, CA.). PCR primers, PCR cycling conditions, and pyrosequencing primers are reported in Additional file 2. All reactions were completed according to manufacturer instructions, except for the following which all included 3.5 nM MgCl2 (provide in the Pyromark PCR kit): MT2, exon 3 Sects. 2–4, exon 3 Sects. 6–8. Standard controls of 0% and 100% methylated DNA, as well as a no DNA control standard, were included for each PCR plate. The median DNA methylation value of the 0% control was 2% (interquartile range: 1%–2%). The median DNA methylation value of the 100% control was 88% (interquartile range: 85%-90%). Following PCR of exon 3 Sect. 1, all samples were run on a 2% agarose gel and a 375 bp product was extracted using the QIAquick Gel Extraction Kit (Qiagen, Valencia, CA) which was subsequently used for pyrosequencing. All samples were randomized for pyrosequencing to account for plate and run variability. Pyrosequencing was performed on a Pyromark Q24 using Pyromark Gold Q24 Reagents (Qiagen, Valencia, CA) per the manufacturer's protocol. Epigenotypes reported are an average of three replicates. The mean deviation from the average ranged from 0.93% (MT2 CpG 7) to 2.29% (MT2 CpG 14).

Oxtr sequencing and SNP analysis

Ten nanograms of DNA was used as a template for PCR using a Pyromark PCR kit (Qiagen, Valencia, CA). In order to determine the alleles at polymorphic sites, pyrosequencing was completed using unknown base calling on a subset of samples using a Pyromark Q24 with PyroMark Gold Q24 Reagents (Qiagen, Valencia, CA) per the manufacturer protocol. Once alleles were identified, pyrosequencing was performed to measure genotypes of each sample. PCR primers, PCR cycling conditions, and pyrosequencing primers are reported in Additional file 2. The same primers were used for both unknown base calling and genotyping assays. SNPs KLW1 and KLW2 were tested for Hardy–Weinberg equilibrium using the HWExact function in the HardyWeinberg R package [45]. Linkage disequilibrium statistics for KLW1, KLW2, and NT213739 were determined using the LD function in the genetics R package [57].

Oxtr expression analysis

Extraction of RNA was done using the Qiagen AllPrep DNA/RNA Mini Kit (Qiagen, Valencia, CA) following manufacturer instructions. RNA was processed for cDNA synthesis following the protocol provided in the iScript cDNA Synthesis kit (Bio-Rad, Hercules, CA). Real-time PCR for the Oxtr main transcript was conducted using a 7500 Fast Real-Time PCR System (Applied Biosystems) using Power SYBR Green (Applied Biosystems No. 4367659). Real-time PCR for the Oxtr alternative transcript was conducted using the CFX96 System (Bio-Rad) using Power SYBR Green (Applied Biosystems). Real-time PCR for Pgk1 was completed on both Real-time PCR systems. See Additional file 2 for all RT-PCR primers and cycling conditions. All reactions were run in triplicate (replicate standard deviation was < 0.05) and their specificity verified by melting curve analysis and separation on a 2% agarose gel. Primer performance was evaluated using standard serial dilution and all primer sets performed within acceptable range for efficiency (See Additional file 2). Relative gene expression is presented using the comparative Ct method, 2−ΔCt, comparing target expression to Pgk1 expression measured on the same real-time PCR system. Pgk1 was chosen as a reference based on data in mouse brain showing its reliability across brain regions and developmental time points [58].

RNA-sequencing and alignment

To identify the structure of Oxtr in prairie voles, we performed transcriptome analysis using RNA-sequencing on a single RNA sample derived from the nucleus accumbens of a female MAN1 vole in our cohort. RNA quality was assessed using an Agilent Tape Station. The RIN score was 9.0. The library was generated from 500 ng RNA using the NEBNext Ultra Directional RNA Library Prep Kit with mRNA magnetic isolation. The UVA Genome Analysis and Technology Core performed paired-end sequencing (2 X 75 bp paired-end run) on the Illumina NextSeq 500 platform. A total of 125 million reads were generated. The raw sequencing data was subjected to preprocessing steps of adapter removal and quality-based trimming using TrimGalore [59] with the removal of auto-detected Illumina adapters and trimming of low-quality ends up to a threshold of Q20. Reads that became shorter than 35 bp due to either adapter removal or quality trimming were discarded. Quality control was completed with MultiQC [60]. Novel transcript identification and quantification was performed using a standard analysis approach outlined elsewhere [61]. Briefly, sequencing data were aligned to the Microtus ochrogaster genome using a splice-aware aligner STAR2 [62], followed by assembly and quantification using Stringtie [63]. Alignments were viewed in the Integrated Genomics Viewer [64].

Statistical analysis

Statistical computing and graphics were generated with R statistical software [65]. Graphics were generated using the ggplot2 and ggpubr packages unless otherwise stated [66, 67]. For each analysis, p < 0.05 was regarded as statistically significant with respect to multiple comparison procedures as appropriate. The effect of handling on DNA methylation in MT2 and exon 3 was determined using a mixed effects model with fixed effect of handling condition and random effects of intercepts for individual and CpG site using the lme4 package [68]. The type II Wald F test in the car package was used to determine significance of fixed effects [69]. The effect of handling condition on DNA methylation at individual CpG sites was determined using t tests followed by Bonferroni correction. In both humans and in prairie voles, Spearman’s rank correlation was used to determine the relationship between DNA methylation and gene expression and uncorrected p values are reported because the distribution of DNA methylation at some CpG sites in human MT2, prairie vole MT2, and prairie vole exon 3 are non-normal. For visualization purposes, we plotted trend lines from linear regression models predicting DNA methylation from expression. The community structure of the MT2 region was determined using exploratory graph analysis (EGA) implemented in the EGAnet package [70]. We estimated the EGA network using the Triangulated Maximally Filtered Graph (TMFG) model [71] as opposed to the graphical LASSO (GLASSO) model, because the GLASSO procedure typically requires n (sample size) >  > p (CpG sites) to accurately compute partial correlations for each site. TMFG skirts this issue by modeling the zero-order correlation matrix, rather than the partial correlation matrix. The effect of handling condition on standardized network community score was determined using a two-way community x handling condition ANOVA, followed by t tests within each community and Bonferroni correction. Spearman’s rho was used to determine the relationship between standardized network community score and Oxtr expression because the distribution of standardized network community scores is non-normal. The impact of Oxtr polymorphisms KLW1 and KLW2 (both A/G SNPs) on gene expression was determined using a t test comparing A carriers (A/A and A/G) to G/G homozygotes. The impact of sex, handling condition, and Oxtr polymorphisms KLW1, KLW2, and NT213739 on expression of the alternative transcript were determined using the Wilcoxon Rank Sum test as expression of the alternative transcript was non-normal.

Availability of data and materials

RNA-sequencing data from the prairie vole nucleus accumbens are available upon reasonable request. All other data generated and analyzed in this study are present in the published article.

References

  1. Lozupone M, La Montagna M, D’Urso F, Daniele A, Greco A, Seripa D, et al. The role of biomarkers in psychiatry. In: Guest PC, editor. Advances in experimental medicine and biology. Cham: Springer; 2019. p. 135–62.

    Google Scholar 

  2. Rutigliano G, Rocchetti M, Paloyelis Y, Gilleen J, Sardella A, Cappucciati M, et al. Peripheral oxytocin and vasopressin: biomarkers of psychiatric disorders? A comprehensive systematic review and preliminary meta-analysis. Psychiatry Res. 2016;241:207–20.

    Article  CAS  PubMed  Google Scholar 

  3. Romano A, Tempesta B, Di Bonaventura MVM, Gaetani S. From autism to eating disorders and more: The role of oxytocin in neuropsychiatric disorders. Front Neurosci. 2016;9:497.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Carter CS, Kenkel WM, MacLean EL, Wilson SR, Perkeybile AM, Yee JR, et al. Is Oxytocin “nature’s medicine”? Pharmacol Rev. 2020;72:829–61.

    Article  PubMed  CAS  Google Scholar 

  5. Zik JB, Roberts DL. The many faces of oxytocin: Implications for psychiatry. Psychiatry Res. 2015;226:31–7.

    Article  CAS  PubMed  Google Scholar 

  6. Jones C, Barrera I, Brothers S, Ring R, Wahlestedt C. Oxytocin and social functioning. Dialogues Clin Neurosci. 2017;19:193–201.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Andari E, Hurlemann R, Young LJ. A precision medicine approach to oxytocin trials. In: Hurlemann R, Grinevich V, editors. Current topics in behavioral neurosciences. Cham: Springer; 2018. p. 559–90.

    Google Scholar 

  8. MacLean EL, Wilson SR, Martin WL, Davis JM, Nazarloo HP, Carter CS. Challenges for measuring oxytocin: the blind men and the elephant? Psychoneuroendocrinology. 2019;107:225–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Berdasco M, Esteller M. Clinical epigenetics: seizing opportunities for translation. Nat Rev Genet. 2019;20:109–27.

    Article  CAS  PubMed  Google Scholar 

  10. Bird A. The dinucleotide CG as a genomic signalling module. J Mol Biol. 2011;409:47–53.

    Article  CAS  PubMed  Google Scholar 

  11. Inoue T, Kimura T, Azuma C, Inazawa J, Takemura M, Kikuchi T, et al. Structural organization of the human oxytocin receptor gene. J Biol Chem. 1994;269:32451–6.

    Article  CAS  PubMed  Google Scholar 

  12. Kusui C, Kimura T, Ogita K, Nakamura H, Matsumura Y, Koyama M, et al. DNA methylation of the human oxytocin receptor gene promoter regulates tissue-specific gene suppression. Biochem Biophys Res Commun. 2001;289:681–6.

    Article  CAS  PubMed  Google Scholar 

  13. Quintana DS, Rokicki J, van der Meer D, Alnæs D, Kaufmann T, Córdova-Palomera A, et al. Oxytocin pathway gene networks in the human brain. Nat Commun. 2019;10:668.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Gregory SG, Connelly JJ, Towers AJ, Johnson J, Biscocho D, Markunas CA, et al. Genomic and epigenetic evidence for oxytocin receptor deficiency in autism. BMC Med. 2009;7:62.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Uhrig S, Hirth N, Broccoli L, von Wilmsdorff M, Bauer M, Sommer C, et al. Reduced oxytocin receptor gene expression and binding sites in different brain regions in schizophrenia: a post-mortem study. Schizophr Res. 2016;177:59–66.

    Article  PubMed  Google Scholar 

  16. Lee MR, Sheskier MB, Farokhnia M, Feng N, Marenco S, Lipska BK, et al. Oxytocin receptor mRNA expression in dorsolateral prefrontal cortex in major psychiatric disorders: a human post-mortem study. Psychoneuroendocrinology. 2018;96:143–7.

    Article  CAS  PubMed  Google Scholar 

  17. Cecil C, Lysenko LJ, Jaffee SR, Pingault J-B, Smith RG, Relton CL, et al. Environmental risk, Oxytocin Receptor Gene (OXTR) methylation and youth callous-unemotional traits: a 13-year longitudinal study. Mol Psychiatry. 2014;19:1071–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Dadds MR, Moul C, Cauchi A, Dobson-Stone C, Hawes DJ, Brennan J, et al. Methylation of the oxytocin receptor gene and oxytocin blood levels in the development of psychopathy. Dev Psychopathol. 2014;26:33–40.

    Article  PubMed  Google Scholar 

  19. Kim YR, Kim JH, Kim MJ, Treasure J. Differential methylation of the oxytocin receptor gene in patients with anorexia nervosa: a pilot study. PLoS ONE. 2014;9:e88673.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Thaler L, Brassard S, Booij L, Kahan E, McGregor K, Labbe A, et al. Methylation of the OXTR gene in women with anorexia nervosa: Relationship to social behavior. Eur Eat Disord Rev. 2020;28:79–86.

    Article  PubMed  Google Scholar 

  21. Rubin LH, Connelly JJ, Reilly JL, Carter CS, Drogos LL, Pournajafi-Nazarloo H, et al. Sex and diagnosis-specific associations between DNA methylation of the oxytocin receptor gene with emotion processing and temporal-limbic and prefrontal brain volumes in psychotic disorders. Biol Psychiatry Cogn Neurosci Neuroimaging. 2016;1:141–51.

    PubMed  Google Scholar 

  22. Bang M, Kang JI, Kim SJ, Park JY, Kim KR, Lee SY, et al. Reduced DNA methylation of the oxytocin receptor gene is associated with anhedonia-asociality in women with recent-onset schizophrenia and ultra-high risk for psychosis. Schizophr Bull. 2019;45:1279–90.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Bell AF, Carter CS, Steer CD, Golding J, Davis JM, Steffen AD, et al. Interaction between oxytocin receptor DNA methylation and genotype is associated with risk of postpartum depression in women without depression in pregnancy. Front Genet. 2015;6:243.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  24. Ebner NC, Lin T, Muradoglu M, Weir DH, Plasencia GM, Lillard TS, et al. Associations between oxytocin receptor gene (OXTR) methylation, plasma oxytocin, and attachment across adulthood. Int J Psychophysiol. 2019;136:22–32.

    Article  PubMed  Google Scholar 

  25. Park CII, Kim HW, Jeon S, Kang JI, Kim SJ. Reduced DNA methylation of the oxytocin receptor gene is associated with obsessive-compulsive disorder. Clin Epigenetics. 2020;12:101.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Smearman EL, Almli LM, Conneely KN, Brody GH, Sales JM, Bradley B, et al. Oxytocin receptor genetic and epigenetic variations: association with child abuse and adult psychiatric symptoms. Child Dev. 2016;87:122–34.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Jack A, Connelly JJ, Morris JP. DNA methylation of the oxytocin receptor gene predicts neural response to ambiguous social stimuli. Front Hum Neurosci. 2012;6:280.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Puglia MH, Lillard TS, Morris JP, Connelly JJ. Epigenetic modification of the oxytocin receptor gene influences the perception of anger and fear in the human brain. Proc Natl Acad Sci. 2015;112:3308–13.

    Article  CAS  PubMed  Google Scholar 

  29. Krol KM, Puglia MH, Morris JP, Connelly JJ, Grossmann T. Epigenetic modification of the oxytocin receptor gene is associated with emotion processing in the infant brain. Dev Cogn Neurosci. 2019;37:100648.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Puglia MH, Connelly JJ, Morris JP. Epigenetic regulation of the oxytocin receptor is associated with neural response during selective social attention. Transl Psychiatry. 2018;8:116.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  31. Chagnon YC, Potvin O, Hudon C, Préville M. DNA methylation and single nucleotide variants in the brain-derived neurotrophic factor (BDNF) and oxytocin receptor (OXTR) genes are associated with anxiety/depression in older women. Front Genet. 2015;6:230.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Ziegler C, Dannlowski U, Bräuer D, Stevens S, Laeger I, Wittmann H, et al. Oxytocin receptor gene methylation: converging multilevel evidence for a role in social anxiety. Neuropsychopharmacology. 2015;40:1528–38.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Cappi C, Belo Diniz J, Requena GL, Lourenço T, Cristina Garcia Lisboa B, Camargo Batistuzzo M, et al. Epigenetic evidence for involvement of the oxytocin receptor gene in obsessive-compulsive disorder. BMC Neurosci. 2016;17:79.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Nawijn L, Krzyzewska IM, van Zuiden M, Henneman P, Koch SBJ, Mul AN, et al. Oxytocin receptor gene methylation in male and female PTSD patients and trauma-exposed controls. Eur Neuropsychopharmacol. 2019;29:147–55.

    Article  CAS  PubMed  Google Scholar 

  35. Rijlaarsdam J, van IJzendoorn MH, Verhulst FC, Jaddoe VWV, Felix JF, Tiemeier H, et al. Prenatal stress exposure, oxytocin receptor gene (OXTR) methylation, and child autistic traits: the moderating role of OXTR rs53576 genotype. Autism Res. 2017;10:430–8.

    Article  PubMed  Google Scholar 

  36. Tabbaa M, Paedae B, Liu Y, Wang Z. Neuropeptide regulation of social attachment: the prairie vole model. Compr Physiol. 2017;7:81–104.

    Google Scholar 

  37. Perkeybile AM, Carter CS, Wroblewski KL, Puglia MH, Kenkel WM, Lillard TS, et al. Early nurture epigenetically tunes the oxytocin receptor. Psychoneuroendocrinology. 2019;99:128–36.

    Article  CAS  PubMed  Google Scholar 

  38. Kenkel WM, Perkeybile AM, Yee JR, Pournajafi-Nazarloo H, Lillard TS, Ferguson EF, et al. Behavioral and epigenetic consequences of oxytocin treatment at birth. Sci Adv. 2019;5:eaav2244.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Towers AJ, Tremblay MW, Chung L, Li X, Bey AL, Zhang W, et al. Epigenetic dysregulation of Oxtr in Tet1-deficient mice has implications for neuropsychiatric disorders. JCI Insight. 2018;6:e120592.

    Article  Google Scholar 

  40. Stothard P. The sequence manipulation suite: JavaScript programs for analyzing and formatting protein and DNA sequences. Biotechniques. 2000;28(6):1102.

    Article  CAS  PubMed  Google Scholar 

  41. Golino HF, Epskamp S. Exploratory graph analysis: a new approach for estimating the number of dimensions in psychological research. PLoS ONE. 2017;12:e0174035.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. Golino H, Shi D, Christensen AP, Garrido LE, Nieto MD, Sadana R, et al. Investigating the performance of exploratory graph analysis and traditional techniques to identify the number of latent factors: a simulation and tutorial. Psychol Methods. 2020;25:292.

    Article  PubMed  Google Scholar 

  43. Christensen A, Golino H. Estimating the stability of the number of factors via Bootstrap Exploratory Graph Analysis: a tutorial. PsyArXiv; 2019. https://doi.org/10.31234/osf.io/9deay.

  44. Christensen AP, Golino H, Silvia PJ. A psychometric network perspective on the validity and validation of personality trait questionnaires. Eur J Pers. 2020;34(6):1095–108.

    Article  Google Scholar 

  45. Graffelman J. Exploring diallelic genetic markers: the HardyWeinberg package. J Stat Softw. 2015;64:1–23.

    Article  Google Scholar 

  46. King LB, Walum H, Inoue K, Eyrich NW, Young LJ. Variation in the oxytocin receptor gene predicts brain region-specific expression and social attachment. Biol Psychiatry. 2016;80:160–9.

    Article  CAS  PubMed  Google Scholar 

  47. Gouin JP, Zhou QQ, Booij L, Boivin M, Côté SM, Hébert M, et al. Associations among oxytocin receptor gene (OXTR) DNA methylation in adulthood, exposure to early life adversity, and childhood trajectories of anxiousness. Sci Rep. 2017;7:1–14.

    Article  CAS  Google Scholar 

  48. Krol KM, Moulder RG, Lillard TS, Grossmann T, Connelly JJ. Epigenetic dynamics in infancy and the impact of maternal engagement. Sci Adv. 2019;5:eaay0680.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Unternaehrer E, Meyer AH, Burkhardt SCA, Dempster E, Staehli S, Theill N, et al. Childhood maternal care is associated with DNA methylation of the genes for brain-derived neurotrophic factor (BDNF) and oxytocin receptor (OXTR) in peripheral blood cells in adult men and women. Stress. 2015;18:451–61.

    Article  PubMed  CAS  Google Scholar 

  50. Feldman R, Monakhov M, Pratt M, Ebstein RP. Oxytocin pathway genes: evolutionary ancient system impacting on human affiliation, sociality, and psychopathology. Biol Psychiat. 2016;79:174–84.

    Article  CAS  PubMed  Google Scholar 

  51. Di Ruscio A, Ebralidze AK, Benoukraf T, Amabile G, Goff LA, Terragni J, et al. DNMT1-interacting RNAs block gene-specific DNA methylation. Nature. 2013;503:371–6.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  52. de la Mora MP, Pérez-Carrera D, Crespo-Ramírez M, Tarakanov A, Fuxe K, Borroto-Escuela DO. Signaling in dopamine D2 receptor-oxytocin receptor heterocomplexes and its relevance for the anxiolytic effects of dopamine and oxytocin interactions in the amygdala of the rat. Biochim Biophys Acta - Mol Basis Dis. 2016;1862:2075–85.

    Article  CAS  Google Scholar 

  53. Romero-Fernandez W, Borroto-Escuela DO, Agnati LF, Fuxe K. Evidence for the existence of dopamine d2-oxytocin receptor heteromers in the ventral and dorsal striatum with facilitatory receptor-receptor interactions. Molecular Psychiatry. 2013;18:849–50.

    Article  CAS  PubMed  Google Scholar 

  54. Oakley RH, Laporte SA, Holt JA, Barak LS, Caron MG. Molecular determinants underlying the formation of stable intracellular G protein-coupled receptor-β-arrestin complexes after receptor endocytosis. J Biol Chem. 2001;276:19452–60.

    Article  CAS  PubMed  Google Scholar 

  55. Lonsdale J, Thomas J, Salvatore M, Phillips R, Lo E, Shad S, et al. The genotype-tissue expression (GTEx) project. Nat Genet. 2013;45:580–5.

    Article  CAS  Google Scholar 

  56. Bales KL, Lewis-Reese AD, Pfeifer LA, Kramer KM, Sue CC. Early experience affects the traits of monogamy in a sexually dimorphic manner. Dev Psychobiol. 2007;49:335–42.

    Article  PubMed  Google Scholar 

  57. Warnes G. genetics: Population Genetics. R package version 1.3.8.1.2. 2019. https://CRAN.R-project.org/package=genetics.

  58. Boda E, Pini A, Hoxha E, Parolisi R, Tempia F. Selection of reference genes for quantitative real-time RT-PCR studies in mouse brain. J Mol Neurosci. 2009;37:238–53.

    Article  CAS  PubMed  Google Scholar 

  59. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011;17(1):10–2.

    Article  Google Scholar 

  60. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT. StringTie and Ballgown Nat Protoc. 2016;11:1650–67.

    Article  CAS  PubMed  Google Scholar 

  62. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  Google Scholar 

  63. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): High-performance genomics data visualization and exploration. Brief Bioinform. 2013;14:178–92.

    Article  PubMed  CAS  Google Scholar 

  65. R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2018. https://www.R-project.org/.

  66. Wickham H. ggplot2: Elegant Graphics for Data Analysis. R package version 3.3.0. 2016. https://ggplot2.tidyverse.org/.

  67. Kassambara A. ggpubr: “ggplot2” Based Publication Ready Plots. R package version 0.1.7. https://CRANR-project.org/package=ggpubr. 2018.

  68. Bates D, Mächler M, Bolker BM, Walker SC. Fitting linear mixed-effects models using lme4. J Stat Softw. 2015;67:1–48.

    Article  Google Scholar 

  69. Fox J, Weisberg S. An R Companion to Applied Regression, (top). Thousand Oaks: Sage; 2019.

    Google Scholar 

  70. Golino H, Christensen A. EGAnet: Exploratory Graph Analysis—A framework for estimating the number of dimensions in multivariate data using network psychometrics. 2020.

  71. Massara GP, Di Matteo T, Aste T. Network filtering for big data: Triangulated maximally filtered graph. J Complex Networks. 2016;5:161–78.

    Google Scholar 

Download references

Acknowledgements

The authors thank the animal care staff at University of California, Davis, for animal care and the staff at the University of Virginia Genome Analysis and Technology Core for completion of RNA-sequencing. We also thank Drs. Krassimira Botcheva and Kathleen Krol for helpful comments on drafts of this manuscript.

Funding

This research was supported by Autism Speaks grant #7110 to J.J.C. and C.S.C., NIH grant HD075750 to C.S.C. and J.J.C., National Alliance for Autism Research grant to C.S.C., NIH grant MH073022 to C.S.C. and K.L.B., NIH grant HD060117 to K.L.B., A.M., and NSF grant 0437523 to K.L.B.

Author information

Author notes

  1. Joshua S. Danoff and Kelly L. Wroblewski have contributed equally to this work

Authors and Affiliations

  1. Department of Psychology, University of Virginia, 102 Gilmer Hall, P.O. Box 400400, Charlottesville, VA, 22904, USA

    Joshua S. Danoff, Kelly L. Wroblewski, Andrew J. Graves, Graham C. Quinn, Travis S. Lillard, Hudson F. Golino & Jessica J. Connelly

  2. The Kinsey Institute, Indiana University, 150 S Woodlawn Avenue, Bloomington, IN, 47405, USA

    Allison M. Perkeybile, William M. Kenkel & C. Sue Carter

  3. Department of Psychological and Brain Sciences, University of Delaware, 105 The Green, Newark, DE, 19716, USA

    William M. Kenkel

  4. Division of Infectious Diseases and International Health, University of Virginia, 345 Crispell Drive, Charlottesville, VA, 22908, USA

    Hardik I. Parikh

  5. Research Computing, University of Virginia, 560 Ray C. Hunt Drive, Charlottesville, VA, 22903, USA

    Hardik I. Parikh

  6. Duke Molecular Physiology Institute, Duke University School of Medicine, 300 N Duke St, Durham, NC, 27701, USA

    Simon G. Gregory

  7. Department of Psychology, University of California, Davis, One Shields Avenue, Davis, CA, 95616, USA

    Karen L. Bales

Authors
  1. Joshua S. Danoff
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Kelly L. Wroblewski
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Andrew J. Graves
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Graham C. Quinn
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Allison M. Perkeybile
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. William M. Kenkel
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Travis S. Lillard
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Hardik I. Parikh
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Hudson F. Golino
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Simon G. Gregory
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. C. Sue Carter
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Karen L. Bales
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Jessica J. Connelly
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

J.S.D., K.L.W., A.M.P., W.M.K., S.G.G., C.S.C., K.L.B., and J.J.C. were responsible for designing research. J.S.D., K.L.W., G.C.Q., T.S.L., A.M.P and J.J.C. conducted the studies. J.S.D., K.L.W., A.J.G., H.I.P., H.F.G., and J.J.C. analyzed and interpreted the data. J.S.D. and J.J.C. wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Jessica J. Connelly.

Ethics declarations

Ethics approval and consent to participate

All procedures involved in generating tissue for the analysis of DNA methylation and gene expression following handling were reviewed and approved by the IACUC at the University of California, Davis.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. Supplementary figures and tables.

Additional file 2

. Primers and PCR conditions for DNA methylation analysis and RT-PCR.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Danoff, J.S., Wroblewski, K.L., Graves, A.J. et al. Genetic, epigenetic, and environmental factors controlling oxytocin receptor gene expression. Clin Epigenet 13, 23 (2021). https://doi.org/10.1186/s13148-021-01017-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13148-021-01017-5

Keywords

Clinical Epigenetics

ISSN: 1868-7083

Contact us