All subjects were laboratory-bred prairie voles (Microtus ochrogaster) derived from wild-caught stock near Champaign, Illinois treated as previously described (Perkeybile et al., 2019). Subjects are at least four generations removed from wild-caught stock. Breeding pairs were housed in polycarbonate cages (44 cm × 22 cm × 16 cm). Same sex sibling pairs were housed in polycarbonate cages (27 cm × 16 cm × 16 cm) after weaning on postnatal day 20 (PND20). Animals were given food (high-fiber Purina rabbit chow) and water ad libitum, cotton nestlets for nesting material, and were maintained on a 14:10 light:dark cycle.

Four nucleus accumbens RNA samples (two male and two female, aged between postnatal days 21 and 24) from our previous study of early life experience were pooled for 5′ RACE (Perkeybile et al., 2019). Uterine samples were taken from two adults: one virgin (115 days old) and one postpartum day 1 (multiparous and over a year old) from our colony at the University of Virginia. Uterine samples were collected from animals that were deeply anaesthetized with an overdose of Euthasol (Covetrus) and sacrificed via cervical dislocation. An incision was made on the abdomen and the uterus was removed and flash frozen on dry ice, then stored at −80°C until nucleic acid isolation. RNA was isolated using an AllPrep DNA/RNA Mini kit (Qiagen) according to manufacturer instructions. All procedures for collection of brain tissue were approved by the IACUC at the University of California, Davis. All procedures for collection of uterine tissue were approved by the IACUC at the University of Virginia.

5′ RACE was completed using the GeneRacer Kit with Superscript III RT (Invitrogen) according to manufacturer instructions. For uterine samples, 3 μg of RNA was used as input. For nucleus accumbens, 2.5 μg, which was pooled from four different samples, was used as input. Briefly, RNA was dephosphorylated and decapped and the GeneRacer oligo was ligated to the mRNA. mRNA was transcribed into cDNA using Superscript III RT using random primers. Touchdown PCR was performed with Platinum Taq DNA Polymerase High Fidelity (Invitrogen) according to manufacturer instructions with an annealing temperature of 65°C, with the provided forward GeneRacer primer, and the reverse gene-specific primer 5′-CGT​CTC​TCC​GGG​CCT​GCT​GCC​CT-3′. A nested PCR was then performed using 1 μL of initial PCR as input, the provided nested forward GeneRacer primer, and the nested reverse gene-specific primer 5′- CCG​GGC​CTG​CTG​CCC​TTC​AGG​T-3′ according to manufacturer instructions. 20 μL of nested PCR product was run on a 1.5% agarose gel with ethidium bromide. Single bands denoted in Figure 1 were cut from the gel and DNA was extracted using the QIAquick Gel Extraction Kit (Qiagen) and resuspended in 10 μL of Buffer EB. Nonspecific PCR bands were identified by: 1) PCR controls missing either template, forward primer, or reverse primer; 2) an experiment where we denatured PCR products and slowly cooled them to prevent PCR product hybrids; and 3) sequencing. We performed PCR and downstream cloning and sequencing using 5′ RACE products on three separate occasions and present sequencing data from all three runs. However, data are compiled into one figure for visualization purposes.

Each gel-extracted band represents a PCR product that was cloned using the TOPO TA Cloning Kit (Invitrogen) according to manufacturer instructions. Briefly, 4 μL of the purified PCR product was cloned into the TOPO vector. One Shot TOP10 E. coli was transformed using the heat shock method. Transformed E. coli were spread on LB-agar plates containing 50 μg/mL ampicillin for selection and incubated at 37°C for 24 h. Three single colonies from each PCR product were picked and grown overnight in liquid LB with 50 μg/mL ampicillin at 37°C shaking at 200 rpm. Cells were pelleted by spinning at 4000 rcf for 10 min and plasmids containing 5′ end clones were isolated using the PureLink Quick Plasmid Miniprep Kit (Invitrogen). Plasmids were sent to Genewiz for sequencing. Sequences were aligned to Oxtr using NCBI’s blastn. Sequence information for each transcript is provided in Supplementary Data.

Methods for generation of samples for the induction of labor study and measurement of DNA methylation in these samples were previously reported and are summarized below (Kenkel et al., 2019). Prairie voles were paired and allowed to mate using a timed-mating paradigm which enabled prediction of the day of birth. On the expected day of birth, pregnant dams were intraperitoneally injected with either saline (n = 9 dams), low dose oxytocin (0.125 mg/kg, n = 6 dams), medium dose oxytocin (0.25 mg/kg, n = 8 dams), or high dose oxytocin (0.5 mg/kg, n = 6 dams). An additional group of pregnant dams were not injected with anything and used as a no treatment control (n = 9 dams). 90 min following injection, pregnant dams were anesthetized with isoflurane and fetuses were removed and euthanized. Fetal brains were extracted and sectioned into thirds (which roughly correspond to forebrain, midbrain, and hindbrain), then they were flash frozen on dry ice. DNA and RNA were extracted from the midbrain of one male and one female fetus per mother using AllPrep mini kits (Qiagen). All procedures for the induction labor study were approved by the IACUC at Northeastern University.

Offspring used for gene expression analysis were restricted to animals with G/G genotype at SNP KLW2. This genotype which has previously been associated with expression of Oxtr-H (Danoff et al., 2021). The number of offspring included in each group following filtering for genotype are: saline treatment: 7 males, 7 females; low dose oxytocin: 6 males, 6 females; medium dose oxytocin: 6 males, 7 females; high dose oxytocin: 5 males, 6 females; no treatment control: 7 males, 8 females. 1 μg RNA from fetal midbrain was processed for cDNA synthesis using the iScript cDNA synthesis kit (Bio-Rad) according to manufacturer instruction. Reverse transcriptase quantitative PCR (RT-qPCR) was completed using the CFX966 system (Bio-Rad) with Power SYBR Green chemistry (Applied Biosystems). Primers were designed to capture distinct transcripts or sets of transcripts. We note that we attempted to generate an assay specific to transcripts Oxtr-K and Oxtr-L, and while PCR was successful and generated the correct amplicon, the assay could not be made efficient enough for quantitation, likely because of low abundance of these transcripts in brain tissue. Expression of all transcripts were normalized to Gapdh using the comparative Ct method, 2^−ΔCt. See Supplementary Table S1 for all primers and PCR conditions. Specificity of products was confirmed using melt curve analysis. All RT-qPCR reactions were completed in triplicate, and the average Ct value of replicates was used for further analysis.

5 ng of DNA was used as template for PCR using a Pyromark PCR kit (Qiagen, Valencia, CA) according to manufacturer instructions using 0.2 μM forward primer 5′-GGA​CCC​AAG​GCC​TAT​GTC​A-3′ and biotinylated reverse primer 5′-ATG​AGC​TTG​ACG​CTA​CTG​ACT​CG-3’. The following cycling conditions were used for amplification of the target fragment: [Step 1: (95°C/15 min)/1 cycle, Step 2: (94°C/30 s, 58°C/30 s, 72°C/30 s)/45 cycles, Step 3: (72°C/10 min)/1 cycle, Step 4: 4°C hold. Pyrosequencing was completed using Pyromark Gold Q24 reagents (Qiagen) according to manufacturer instructions with sequencing primer 5′-TCA​AGA​TCT​GGC​AGA​AC-3’.

Methods for measurement of DNA methylation in induction of labor samples have previously been reported (Kenkel et al., 2019). Briefly, 200 ng genomic DNA were subject to bisulfite treatment (Kit MECOV50, Invitrogen, Carlsbad, CA) per manufacturer instructions. 12 ng of bisulfite converted DNA was used as a template for PCR using a Pyromark PCR kit (Qiagen) and 0.2 uM of primers 5′-GGG​GAT​AGG​ATG​GTT​AGT​TAG​TAT​T-3′ and 5′- CCA​ACA​ACC​TCA​AAA​CTC​TAC​T-3’. The following cycling conditions were used for amplification of the target Oxtr fragment, which included CpG sites −934_1, −934_2, −924, and −901: [Step 1: (95°C/15 min)/1 cycle, Step 2: (94°C/30 s, 58°C/30 s, 72°C/30 s)/50 cycles, Step 3: (72°C/10 min)/1 cycle, Step 4: 4°C hold. Pyrosequencing was performed using two primers: 5′-GAG​GGA​AGG​TTT​TGG​AGT​TTT​TTA​TAT-3′, which covers CpG sites −934_1, −934_2, and −924, and 5′-AGG​GAT​TGA​AAA​GTG​A-3′, which covers CpG site −901, on a Pyromark Q24 using PyroMark Gold Q24 Reagents (Qiagen) per the manufacturer protocol. Epigenotypes reported are an average of the three replicates.

Statistical analysis was completed in R statistical software. RT-qPCR data was analyzed using mixed-effect linear models with fixed effects of a treatment by sex interaction and random effects of litter ID to account for genetic relatedness using the afex package in R (Team, 2018; Singmann et al., 2022). Post-hoc tests were completed used the emmeans R package and a modified Dunnett correction (“dunnetx”) for multiple comparisons (Lenth, 2022). Comparison of RT-qPCR data from the no treatment and saline groups were completed using the same model structure. Since no sex or treatment by sex effects were associated with expression of any transcript, the relationship between DNA methylation and Oxtr transcript expression was examined using mixed-effect linear models with fixed effects of DNA methylation at a specific CpG site and treatment and random effects of litter ID. Separate models were constructed for each CpG site measured to avoid collinearity.

Using 5′ RACE, we identified seven Oxtr transcripts with unique 5′ transcription start sites in prairie vole nucleus accumbens and uterine tissue (Figures 1A, B). Though we are interested in transcriptional products of Oxtr in the brain, we also used two uterine samples, a virgin uterus and a postpartum uterus (day of birth), because of the abundance of Oxtr expression in this tissue. Four of these transcripts were previously identified in mouse hippocampus (Towers et al., 2018) and are named according to their homology. These include Oxtr-A, which is considered the canonical transcript, consisting of four exons and three introns. Oxtr-B consists of a 5′ extended exon 3 spliced to exon 4. Oxtr-E consists of a truncated exon 3 spliced to exon 4. Oxtr-H consists of a novel exon originating from the third intron of the gene spliced to exon 4. Three novel transcripts were also identified. Oxtr-J which originates in the middle of exon 3 and is spliced to exon 4. Oxtr-K and Oxtr-L both consist of two novel exons about 4,000 basepairs upstream of the canonical transcription start site, which get spliced to exon 2. Oxtr-L differs from Oxtr-K in that first exon is 35 bp shorter and the intron between exon 2 and 3 is retained, though we cannot rule out the possibility that Oxtr-L represents a transcript similar to Oxtr-K that has not been fully spliced. Sequences of each transcript are provided in Supplementary Data.

We predicted putative translation start sites using the ATG in the transcript sequence that produces the longest peptide sequence, which was always in frame. We compared this predicted sequence to the prairie vole OXTR protein sequence (Uniprot E0V872) to generate schematics of the predicted protein. The predicted protein products of these transcripts are shown in Figure 1C. Oxtr-A, Oxtr-B, Oxtr-K, and Oxtr-L would all make a full length OXTR protein if translated. The other transcripts could be translated to truncated OXTR proteins as shown in Figure 1C. Importantly, we do not yet have evidence that these transcripts are all translated; it is possible that some transcripts are not translated and possess gene regulatory functions or do not have any biological function.

We investigated how transcriptional dynamics of Oxtr isoforms in the fetal brain are impacted by biological context using a model of labor induction (Kenkel et al., 2019). In this study, we injected pregnant dams with saline, low dose oxytocin (0.125 mg/kg), medium dose oxytocin (0.25 mg/kg), or high dose oxytocin (0.5 mg/kg) on the expected day of birth. An additional group of pregnant dams were not injected and used as a no treatment control. 90 min after the injection, fetuses were extracted and the fetal brain was removed. Our previous results indicate that maternal oxytocin treatment leads to increased Oxtr DNA methylation and impacts Oxtr expression (Kenkel et al., 2019). We suspected that a relationship between maternal oxytocin exposure and Oxtr gene expression in the fetal brain might be further clarified by examining single transcript specific expression. We investigated this using reverse transcriptase quantitative PCR (RT-qPCR) with primers that can detect distinct alternative transcripts (Oxtr-A, Oxtr-B, Oxtr-H) or groups of transcripts (Oxtr-All, Oxtr-All except H) (Figure 2A).

Our results indicate that there are significant group differences in expression when using primers that capture all Oxtr transcripts, but post hoc tests do not indicate that any specific group differs from saline [F_(3,22.5) = 3.66, p = 0.028, Figure 2B]. When using primers that capture all Oxtr transcripts except Oxtr-H, there is again a significant group difference, with the medium group having increased expression compared to saline [F_(3,22.4) = 4.03, p = 0.020, Figure 2C]. We then examined if expression of only specific Oxtr transcripts is altered by maternally administered oxytocin. None of the treatment group expression differences in Oxtr transcripts can be attributed to Oxtr-A, which does not change expression levels with maternal oxytocin administration [F_(3,22.5) = 1.30, p = 0.229, Figure 2D]. However, these changes appear to be driven by changes in expression of Oxtr-B, which does significantly differ by maternal oxytocin administration, and is elevated in the medium dose and high dose groups [F_(3,21.7) = 3.50, p = 0.033, Figure 2E]. We also identify group differences in expression of Oxtr-H, where expression is elevated in the high dose group [F_(3,42) = 4.49, p = 0.008, Figure 2F]. Importantly, expression of each of these transcripts or groups of transcripts does not differ between offspring of saline treated dams and offspring of dams receiving no injection (Supplementary Figure S1). These results indicate that maternal oxytocin treatment alters expression of specific Oxtr transcripts in the fetal brain in a dose-dependent manner and suggests the importance of clarifying expression of these specific transcripts in previous literature in the vole and the mouse. Moreover, these results underscore the importance of precise quantification of specific transcripts when examining how endogenous or exogenous factors impact Oxtr expression in the brain.

DNA methylation is an epigenetic modification that occurs when a methyl group is covalently added to cytosine, typically in a cytosine-phosphate-guanine (CpG) context (Bird, 2011). DNA methylation in gene promoters is generally associated with transcriptional repression (Bird, 2011). DNA methylation in a regulatory region of Oxtr called MT2, which is located within exon 1 and intron 1 of Oxtr, is associated with reduced transcription of OXTR in both humans and prairie voles (Kusui et al., 2001; Danoff et al., 2021). We tested if DNA methylation in MT2 regulates specific Oxtr transcripts. We examined DNA methylation at four CpG sites in MT2 (CpGs −934_1, −934_2, −924, and −901) which we have previously shown to be sensitive to early life experience and are among the CpG sites where methylation is most highly associated with Oxtr expression (Danoff et al., 2021). Results of all DNA methylation models are in Supplementary Table S2 and indicate that all CpG sites measured show similar relationships with Oxtr transcript expression. For simplicity, results from CpG site −924 are provided below and in Figure 3. DNA methylation at CpG site −924 is negatively associated with expression of Oxtr when the assay does not distinguish between individual transcripts [F_(1,51.6) = 12.55, p < 0.001; Figure 3A] and Oxtr expression when the assay measures all transcripts except H [F_(1,58.7) = 11.28, p = 0.001; Figure 3B]. We next measured the expression of specific transcripts to identify which transcripts are associated with MT2 DNA methylation. In the fetal brain, the relationship between Oxtr gene expression and DNA methylation cannot be attributed to regulation of Oxtr-A, which is not associated with DNA methylation at CpG −924 [F_(1,59.0) = 1.29, p = 0.26; Figure 3C]. However, DNA methylation at CpG −924 is negatively associated with expression of Oxtr-B [F_(1,54.4) = 14.15, p < 0.001; Figure 2D]. DNA methylation at CpG −924 is not associated with expression of Oxtr-H [F_(1,59.0) = 2.44, p = 0.123; Figure 2E]. These results indicate that DNA methylation in MT2 impacts expression of Oxtr-B in the fetal brain.