Experience-dependent neuroplasticity of the developing hypothalamus: integrative epigenomic approaches

ABSTRACT

Augmented maternal care during the first postnatal week promotes life-long stress resilience and improved memory compared with the outcome of routine rearing conditions. Recent evidence suggests that this programming commences with altered synaptic connectivity of stress sensitive hypothalamic neurons. However, the epigenomic basis of the long-lived consequences is not well understood. Here, we employed whole-genome bisulfite sequencing (WGBS), RNA-sequencing (RNA-seq), and a multiplex microRNA (miRNA) assay to examine the effects of augmented maternal care on DNA cytosine methylation, gene expression, and miRNA expression. A total of 9,439 differentially methylated regions (DMRs) associated with augmented maternal care were identified in male offspring hypothalamus, as well as a modest but significant decrease in global DNA methylation. Differentially methylated and expressed genes were enriched for functions in neurotransmission, neurodevelopment, protein synthesis, and oxidative phosphorylation, as well as known stress response genes. Twenty prioritized genes were identified as highly relevant to the stress resiliency phenotype. This combined unbiased approach enabled the discovery of novel genes and gene pathways that advance our understanding of the epigenomic mechanisms underlying the effects of maternal care on the developing brain.

KEYWORDS:

• Augmented maternal care
• epigenomics
• DNA methylation
• early-life stress
• hypothalamus
• resilience

Introduction

Early postnatal care in mammals is critical for survival as well as neurodevelopment. Specifically, early maternal care can epigenetically program the behavior of offspring for their entire lifetime. Rodent models that manipulate maternal care and maternal stress have been used to examine the molecular and cellular mechanisms underlying the long-term behavioral changes observed in the offspring [1–3]. In the augmented maternal care (AMC) paradigm, daily separation of rat dams from their pups for 15 minutes from postnatal day 2 (PND2) to postnatal day 8 (PND8) triggers the increased maternal behaviors of arched-back nursing, as well as elevated licking and grooming of the pups by the dam [3,4]. Subsequently, it was found that multiple, continuous days of pup handling within the first postnatal week of life were required for the AMC response, suggesting that AMC initiates changes during this critical period of neurodevelopment in the offspring [5,6]. The repeated, daily sessions of enhanced maternal care triggered by the AMC paradigm result in offspring that show reduced plasma corticosterone levels in response to an acute restraint stress [7–11], reduced anxiety related behaviors [11,12], resilience to depressive-like behaviors [11], and enhanced performance on several cognitive tasks as adults [12,13].

AMC appears to alter brain function for the lifetime of the animal [14]. Adult males experiencing early-life AMC are stress resilient, which is manifested as reduced corticosterone response to stress [6,13] and improved cognitive function [6]. The current study focused on epigenetic processes that underlie a well-established phenotypic outcome of AMC. To date, our group has established the phenotype in males only [6,13]. Therefore, we focused our investigations to male hypothalamus to understand the molecular and epigenetic mechanisms that drive the cognitive and emotional phenotype observed in males. Previous work has demonstrated that this mechanism depends upon the temporally coordinated regulation of both corticotrophin releasing hormone (CRH) in the paraventricular nucleus (PVN) of the hypothalamus [5,13,16] and glucocorticoid receptor (GR, encoded by Nr3c1) [6] in the hippocampus [13,15]. Specifically, Crh transcripts in the hypothalamus are significantly reduced immediately following AMC at PND9 [5,6] and remain repressed throughout adulthood [9]. However, increased Nr3c1 expression in the hippocampus was delayed until PND23 [9], suggesting that early gene expression changes in the PVN of the hypothalamus drive the initial molecular mechanisms underlying AMC. Together, the elevated hippocampal Nr3c1 expression and reduced hypothalamic Crh expression [2,5,11,16] in AMC adult offspring suggest that increased negative feedback signaling may underlie the observed hypo-responsive HPA axis. These studies have raised an important question about how hypothalamic reprogramming is initiated and maintained by AMC induced early-life gene expression changes. AMC initiates transient structural changes in CRH-expressing cells in the hypothalamus of pups [12,16]. By PND9 there is decreased excitatory synaptic drive onto hypothalamic Crh expressing neurons in AMC offspring, that returns to baseline by adulthood (PND30-45) [16], suggesting that the AMC induced reduction in Crh expression is maintained long-term by additional mechanisms [14]. In vitro studies have suggested that the synaptic changes initiate transcriptional regulation in these CRH positive cells [12]. Specifically, recent work has implicated the neuron restrictive silencing factor (NRSF) in establishing long-term epigenetic gene expression patterns resulting from AMC [14,17]. NRSF expression and binding to the Crh gene is increased in hypothalamus of AMC offspring at PND9 but is no longer present in adulthood [12,16]. Thus, there is little information about the comprehensive epigenomic processes by which augmented sensory signals from the mother result in enduring, large-scale reprogramming of gene expression in the hypothalamus. In hippocampus, prior studies have described epigenetic DNA methylation alterations and H3K9 acetylation in response to differences in maternal care at Nr3c1 [18]. The search for additional genes underlying AMC effects on adult hippocampus revealed alterations to gene expression of more than 900 genes related to cellular metabolism and energy production, signal transduction, protein synthesis (mainly ribosomal genes), and neurodevelopment [19]. Targeted analysis of a seven megabase (Mb) region of chromosome 18 flanking Nr3c1 revealed that the clustered protocadherins (Pcdh) showed the highest differential response in DNA methylation [20].

Together, these combined results indicate that NRSF binding in PND9 hypothalamus [12,16] as well as differential gene expression and DNA methylation in adult hippocampus are involved in regulating long-term behavioral changes induced by AMC. Whereas NRSF ChIP-seq has identified gene sets contributing to AMC-induced cellular changes [12], a multi-methodology genome-wide approach has not been undertaken in the hypothalamus [21]. To address this gap in knowledge we performed an integrated, genome-wide characterization of epigenetic and gene expression differences in hypothalamus from PND9 pups in response to AMC. These experiments identified that the response to AMC occurs not only in the genes previously identified, such as Crh, but also in thousands of genes throughout the genome. In addition to a modest but significant global reduction in hypothalamic DNA methylation levels, AMC induces differential methylation and expression in mRNA transcripts related to neurodevelopment, synaptic signaling, ribosome function, and cellular stress. The results provide strong evidence for the role of epigenomic modifications in regulating gene expression profiles related to resilience and offer a suite of prioritized genes and gene pathways for the development of diagnostic biomarkers and therapies for neuropsychiatric disorders in humans.

Methods

Detailed methods are provided in the Supplemental Information section. All animal studies were performed in accordance with UC Davis IACUC protocol #18745. Augmented Maternal Care (AMC) was performed as described previously [13,16]. Briefly, pups were gathered at birth and five male and five female pups were randomly assigned to dams. Pups were removed for 15 min at 9 AM (PST) each day from PND2-PND9 to a separate cage in a different room equipped with a heating pad. Levels of maternal care varied between dams with three dams averaging 1069 ± 129 (SEM) minutes licking and grooming their pups in AMC litters compared with three control dams that averaged 482.66 ± 35 minutes licking and grooming their pups. Molecular analyses were performed on samples from one AMC and one control litter. Since the stress-resilient phenotype has been robustly established for males only in the literature [6,12,13,17], only male hypothalami were evaluated for genome-wide scans. Further details can be found in the Supplemental Information (SI) Materials and Methods. For WGBS DNA methylation analysis, genomic DNA was converted by bisulfite treatment using a Zymo kit followed by Illumina sequencing library preparation. WGBS libraries were sequenced on an Illumina Hi-Seq 2500. DMR identification and PMD analysis was performed as detailed in Schroeder [22] and Dunaway [23]. The SI Materials and Methods section contain an in depth description of WGBS methods. Bisulfite converted genomic DNA was amplified using custom, gene specific primers; amplicons were analyzed by pyrosequencing on a Biotage PSQ-96MA Pyrosequencer. A detailed description of pyrosequencing methods is contained in the SI Materials and Methods. For RNA-seq total RNA was prepared from tissue using a Qiagen RNeasy Kit (Qiagen). Total RNA was depleted of ribosomal RNA for creation of stranded RNA-seq libraries (See SI Materials and Methods). Bar coded RNA-seq libraries were sequenced on an Illumina Hi-Seq 2500 platform. Analysis of RNA-seq reads is described in SI Materials and Methods. Validation for differential RNA-seq identified genes was performed on samples from the same cohort of PND9 male and female hypothalami (SI Materials and Methods). For miRNA transcript analysis a multiplex miRNA assay of miRNA transcripts was performed using an nCounter Rat miRNA expression assay (NanoString). A detailed description of miRNA analysis is available in SI Materials and Methods. All sequencing data will be available through the NIH GEO database accession number GSE99282.

Results

AMC induces large-scale hypomethylation in the PND9 hypothalamus

To investigate genome-wide DNA methylation changes associated with AMC, hypothalami from AMC and control male rat pups were chosen for whole-genome bisulfite sequencing (WGBS) based on significant (P = 0.0042) differential AMC (total licking and grooming) at postnatal day 9 (Supplemental Figure 1). Principal component analysis (PCA) was performed on 20,000 base (20-Kb) windows of percent DNA methylation with CpG islands masked (Figure 1 A) to explore genome-wide patterns of methylation between treatment conditions in the individual rats. Methylation levels in individual PND9 hypothalami were separated by AMC condition for principal component 1 (30.2% of variance explained), indicating differences in methylation levels between conditions. Unexpectedly, hypothalami from AMC pups displayed significant global hypomethylation (P = 0.015) at PND 9 (Figure 1 B). Specifically, global CpG methylation levels for augmented hypothalami averaged 76.08% while controls averaged 79.20%, yielding to a 3.12% reduction (Supplemental Table 1). Non-CpG methylation levels in augmented and control hypothalami ranged from 0.86 to 1.46% of total cytosines (Supplemental Table 1) consistent with previous findings of normal neurodevelopment [24]. CpG hypomethylation was significant across all chromosomes with the exception of chromosome 13 and X, indicating a genome-wide impact of AMC on CpG methylation levels (Figure 1 C). Interestingly, windowing over different genomic features revealed hypomethylation specifically in the region 5 kb upstream of the transcriptional start site (TSS) of CpG island promoters (P = 0.0057) in AMC hypothalamus (Supplemental Figure 2). To identify potential hypothalamic genes and pathways contributing to AMC, large-scale differential DNA methylation patterns were identified from the WGBS data sets. Partially methylated domains (PMDs) were mapped from AMC and control WGBS data sets using a hidden Markov model (HMM) previously described for detecting partially methylated domains in neuronal cell lines [22]. AMC PMDs averaged 7.3% of the genome and control PMDs averaged 6.7% of the genome. Although consistent with genome-wide hypomethylation the increased PMD coverage with AMC was not significant (P = 0.158) (Supplemental Table 2). While PMDs were found throughout the genomes of both groups, only five PMDs were differential by treatment group (Table 1). Consistent with genome-wide hypomethylation shown in Figure 1 B, these five PMDs were found in AMC but not control hypothalami. AMC PMDs covered 674,924 bases and overlapped with seven genes. These PMD genes included Dbn1, Prr7, Foxo3, and Cnksr2, which function in neuronal development, and Pak6, which functions in the MAP kinase pathway to regulate cytoskeletal dynamics underlying multiple brain processes [25–29] (Table 1). AMC PMDs also contained AY383691 and Ankrd63, which are gene products without well-established functions.

Figure 1. Reduced DNA methylation in hypothalamus is associated with AMC. (A) Principal components 1 and 2 of the average DNA methylation in 20 kb windows with masking of CpG Islands separately groups AMC samples (red outline, n = 3) from control samples (blue outline, n = 3). Outlines represent 95% confidence intervals. (B) Average genome wide CpG methylation is significantly reduced (P = 0.015) by 3.12% in AMC samples (n = 3) compared with control samples (n = 3). (C) DNA hypomethylation in the AMC hypothalami at PND9 was significant (*P < 0.05) on every chromosome except chromosome 13 and X (P < 0.52 and 0.054, respectively). RM ANOVA AMC group X-chromosome interaction F (20.80) = 21.89, P < 0.0001; * Benjamini-Hochberg corrected t-tests for each chromosome all P < 0.05.

Table 1. Gene transcripts located within AMC PMDs.

size (bases)	genomic position	genes	PMD	function
24,869	chr17:9,703,631-9,678,703	Dbn1	AMC	Cell migration, extension of neuronal processes and plasticity of dendrites
		Prr7	AMC	Enriched in post-synaptic densities and dendritic spines
4,862	chr20:46,425,321-46,430,182	Foxo3	AMC	Transcriptional activator regulating neuronal apoptosis in the absence of survival factors upon oxidative stress
11,746	chr3:110,484,332-110,496,077	Ankrd63	AMC	Ankyrin repeat domain protein 63
		Pak6	AMC	Functions in cytoskeleton rearrangement, apoptosis, and the mitogen-activated protein (MAP) kinase signaling pathway
168,454	chrX:39,544,850-39,713,303	Cnksr2	AMC	Ternary complex assembly of synaptic proteins and coupling of signal transduction to membrane/cytoskeletal remodeling
464,993	chr5:35,189,670-35,654,662	AY383691	AMC	unknown

AMC DMRs are linked to neurodevelopmental signaling and neurotransmission

While PMDs highlight broad regions of differential methylation associated with AMC, few genes are located within these loci. Therefore, to identify small scale (<3 kb) methylation changes, differentially methylated regions (DMRs) were identified from WGBS data [23,30]. With this bioinformatics approach, a total of 9,439 DMRs between AMC and Control PND9 hypothalami were identified (Supplemental Table 3). Of these 9,439 AMC DMRs, 6,228 had reduced methylation levels in the AMC samples, which is consistent with the genome-wide hypomethylation shown in Figure 1. To investigate the relationship between large-scale PMDs and small-scale DMRs, the two genomic features were compared. Surprisingly, of the 9,439 DMRs identified, only one in Pak6 was found within the five AMC PMDs. Visualization of the DMRs using hierarchical clustering analysis (HCA) of significant (P < 0.05) DMRs from hypothalami clustered each rat by treatment (Figure 2 A; Supplemental Table 3). Both PCA of 20 kb windows (Figure 1 A) and HCA (Figure 2 A) independently partitioned the AMC and control hypothalami into separate groups, indicating that distinct methylation patterns distinguish the AMC from control hypothalamic samples at two different resolutions. Motif enrichment analysis of the DMRs revealed the top significant (q < 0.0001) enrichment was for CTCF binding sites, which were present in 4% of DMRs (Figure 2 B; Supplemental Table 4). CTCF binds these sites in a methylation sensitive manner and functions as a genomic insulator critical to chromatin architecture [31]. Of the 9,439 total DMRs identified, 5,284 were between 5 kb upstream and 1 kb to downstream of 4,023 coding and noncoding RefSeq genes in the rat genome (Rn6). Forty-four percent of the DMRs were intergenic and 35.3% were localized to introns. Fourteen percent of DMRs were immediately upstream of genes and only 0.4% and 1.1% were located in the 5’ and 3’ prime untranslated regions (UTRs) respectively (Supplemental Figure 3A). For the DMRs annotated to within genes, the majority (93.6%) were localized to coding genes and 4.1% were located in long noncoding RNAs (lncRNAs) (Supplemental Figure 3B).

Figure 2. Differences in DNA methylation between AMC (n = 3) and control (n = 3) P9 hypothalami. (A) Hierarchical clustering of significant (P < 0.05) DMRs. Individual percent methylation values are normalized to the mean of each DMR. (B) CTCF motif binding motif was significantly (q < 0.0001) enriched in DMRs with 4.24% of DMRs containing the binding motif compared to 3.15% of background sequences with the motif (q < 0.00001). (C) Significant (q < 0.05) enriched gene ontologies and pathways of genes mapped to significant (P < 0.05) DMRs. Overlap corresponds to genes observed compared with the total genes in the pathway.

To identify the functions associated with the 4,023 genes mapped to AMC DMRs, gene ontology (GO) and pathways analyses were performed (Figure 2 C). Gene ontology analysis revealed a significant (q < 0.05) enrichment for the biological processes of axon guidance, behavior, and protein auto-phosphorylation. Membrane raft, adherens junction, and synapse were the significantly (q < 0.05) enriched cellular components. Protein serine/threonine kinase activity, receptor signaling protein activity, and kinase binding were the most unique significant (q < 0.05) molecular functions. Pathway analyses revealed a significant (q < 0.05) enrichment for protein kinase A signaling, axonal guidance signaling, epithelial adherens junction signaling, Wnt/β-catenin signaling, and molecular mechanisms of cancer. Notably, the molecular mechanisms of cancer pathway represent a collection of developmental signaling pathways. Overall, all of the above GO terms and pathways are related to neurodevelopment as well as synaptic and cellular signaling (Figure 2 C).

To further examine the DMRs for relevance to AMC, the 9,439 DMRs were subject to permutation testing, which yielded seven high confidence or “gold” DMRs (family-wise error rate (FWER) <0.05) (Supplemental Table 5). These gold DMRs were annotated to LOC00911370, Tuba3b, Glb1l, Papolb, Grb10, Tenm3, and Zfp332a (Supplemental Table 5). Zfp322a and Tenm3 both exhibited hypomethylation and increased gene expression in AMC pups. Interestingly, Tenm3 plays a critical role in axon guidance and Zfp322a is a transcription factor that inhibits cellular differentiation [32,33]. The hypomethylation of the Tenm3 associated DMR was validated (P = 0.01) by an independent pyrosequencing assay in a validation cohort (Supplemental Figure 4).

Differential AMC gene expression influences translation and oxidative phosphorylation

To investigate the regulatory relationship between DMRs and gene expression, RNA-seq analysis was performed on postnatal day 9 (PND9) male hypothalami (Figures 3 A and B & Supplemental Table 6). Examination of the top 500 most variably expressed genes revealed clustering according to AMC condition (Figure 3 A). For differential analysis, a total of 2,464 significant (q < 0.05) differentially expressed (DE) transcripts were identified between AMC and control hypothalami (Figure 3 B, 3 C & Supplemental Table 7A & 7B). There were 1,411 genes with reduced expression and 1,053 with increased expression following AMC. Log₂ fold changes in gene transcripts ranged from 2.67 to -2.08 (Supplemental Table 7A & 7B) and the top 20 significant (q < 0.05) genes ranked by fold change are shown in Table 2. Consistent with previous studies that have established a key role for stress-induced neuropeptide Crh [5,11], we also found that Gal, Oxt, and Ucn3, transcripts encoding galanin, oxytocin, and the Crh homologue urocortin, respectively, were significantly reduced in handled pup hypothalamus as were transcripts for the Crh1 receptor (P < 0.05)(Table 2 & Figure 3 B).

Figure 3. Differences in gene expression between AMC (n = 3) and control (n = 3) PND9 hypothalami. (A) Principal components analysis (PCA) of the top 500 most variably expressed genes. Outlines represent 95% confidence intervals. (B) Scatter plot DE transcripts with the top DE transcripts labeled. Grey indicates no significant difference in expression, while red indicates a significant (q < 0.05) difference. (C) Hierarchal clustering analysis (HCA) of significant (q < 0.05) DE transcripts. (D) Significant (q < 0.05) enriched gene ontologies and pathways from significantly (q < 0.05) DE transcripts. (E) Most significant (q < 0.0001) transcription factor binding site identified by motif enrichment analysis.

Table 2. The top 20 elevated and repressed AMC RNA transcripts.

Gene Symbol	q-value	log2 FC	Gene Name
Gpr151	0.002	2.1	G Protein-Coupled Receptor 151
Ebf1	0.002	1.7	early B-cell factor 1
Rbm46	0.002	1.7	RNA Binding Motif Protein 46
Stk17b	0.002	1.7	Serine/Threonine Kinase 17b
Snhg4	0.006	1.6	Small Nucleolar RNA Host Gene 4
Dcun1d1	0.002	1.6	Defective In Cullin Neddylation 1 Domain Containing 1
Zfp347	0.002	1.6	Zinc finger protein 347
Slc5a7	0.002	1.6	solute carrier family 5 member 7
Rgs4	0.003	1.5	Regulator Of G-Protein Signaling 4
Pabpc2	0.008	1.5	poly(A) binding protein, cytoplasmic 2
Sh3bgrl	0.002	1.5	SH3 Domain Binding Glutamate Rich Protein Like
Lpar4	0.002	1.5	Lysophosphatidic Acid Receptor 4
Tlr4	0.005	1.5	Toll Like Receptor 4
Rybp	0.002	1.5	RING1 And YY1 Binding Protein
Thoc2	0.002	1.5	THO Complex 2
Mier1	0.002	1.5	Mesoderm Induction Early Response 1, Transcriptional Regulator
Wfdc1	0.009	1.5	WAP Four-Disulfide Core Domain 1
Stag2	0.002	1.5	Stromal Antigen 2
Arrdc3	0.002	1.4	Arrestin Domain Containing 3
Slc6a4	0.002	1.4	Solute carrier family 6 member 4 / SERT
Gal	0.002	−2.7	Galanin And GMAP Prepropeptide
Dgcr8	0.002	−2.5	DGCR8 Microprocessor Complex Subunit
Rsph1	0.002	−2.4	Radial Spoke Head 1 Homolog
Ucn3	0.005	−2.4	Urocortin 3
Upk2	0.023	−2.4	Uroplakin 2
Oxt	0.002	−2.1	Oxytocin/Neurophysin I Prepropeptide
Gstm1	0.002	−1.9	Glutathione S-Transferase Mu 1
Cox4i2	0.002	−1.9	Cytochrome c oxidase (COX)
Dad1	0.002	−1.8	Defender Against Cell Death 1
Sept1	0.005	−1.8	Septin 1
Gstm6	0.002	−1.8	Glutathione S-Transferase Mu 6
Ercc2	0.003	−1.7	Excision Repair Cross-Complementation Group 2
Tekt4	0.002	−1.7	Tektin 4
Pvalb	0.002	−1.7	Parvalbumin
Ocm2	0.003	−1.7	Oncomodulin 2 (aka Parvalbumin Beta)
Apoe	0.002	−1.6	Apolipoprotein E
Uqcrb	0.002	−1.6	Ubiquinol-Cytochrome C Reductase Binding Protein
Crip1	0.002	−1.6	Cysteine-rich intestinal protein
Ghrh	0.002	−1.6	Growth Hormone Releasing Hormone
Rax	0.039	−1.6	Retina And Anterior Neural Fold Homeobox

To identify the functions associated with the 2,464 transcripts, gene ontology (GO) and pathways analyses were performed (Figure 3 D). Gene ontology analysis revealed a significant (q < 0.05) enrichment for the biological processes of translation, protein localization to organelle, and response to oxidative stress. Mitochondrial membrane, focal adhesion, and synapse were the significantly (q < 0.05) enriched cellular components. Structural constituent of ribosome, peptidase activator activity, and oxidoreductase activity acting on NAD(P)H were the most unique significant (q < 0.05) molecular functions. Pathway analyses revealed a significant (q < 0.05) enrichment for eIF2 signaling, regulation of eIF4 and p760K signaling, oxidative phosphorylation, mitochondrial dysfunction, and mTOR signaling (Figure 3 D). Notably, most of the ontologies and pathways are related to the output of PI3K/AKT/mTOR/PTEN signaling [34], particularly translation by the ribosome (eIF2 and eIF4 signaling) and oxidative phosphorylation in the mitochondria. Motif enrichment analysis identified several transcription factor binding sites, with the top significant (q < 0.0001) enrichment representing binding sites for the transcription factor Elk1 that were present in 31% of transcripts (Supplemental Table 7). Elk1 is part of an insulin signaling network that cross-talks with PI3K/AKT/mTOR/PTEN signaling [34] and is also involved in synaptic plasticity, learning, and neurodegenerative disease [35]. Together, these results reveal differential regulation of a variety of cellular processes, and signaling pathways not previously associated with AMC, which are consistent with neurodevelopmental events occurring in PND9 hypothalamus during AMC [34].

To validate and expand on DE genes identified by genome-wide RNA-seq analysis, reverse transcriptase quantitative PCR (RT-qPCR) analysis was performed on select genes of interest in both males and females as well as a second, independent cohort of AMC offspring. RT-qPCR analysis confirmed that Crh (Supplemental Figure 5A) was significantly reduced in AMC hypothalamus (P = 0.0006). There was a trend (P = 0.07) for an increased expression of Ube3a ligase domain containing isoforms (Supplemental Figure 5B), which were discovered using RNA-seq. There was also a significant (P < 0.0001) reduction in transcripts containing the alternative 3’UTR of Ube3a1 (Supplemental Figure 5C), which has the noncoding function of a competing endogenous RNA (ceRNA) or “miRNA sponge” that is critical to neurodevelopment [36].

Reciprocal differences in miRNA and target gene expression

One of the top gene expression changes was decreased expression of the Dgcr8 microprocessor complex subunit gene (Table 2). Dgcr8, in combination with Drosha, is involved in the processing of primary miRNA (pri-miRNA) transcripts into precursor miRNAs (pre-miRNA) [37]. The decreased expression of Dgcr8 suggests alterations to miRNA processing. Also suggestive of alterations to miRNA regulation was the decrease in the Ube3a1 ceRNA. Since the RNA-seq data from PND9 hypothalamus did not accurately profile miRNA transcripts, we next evaluated the impact of AMC on miRNA expression using a multiplex miRNA expression assay on PND9 hypothalamus from handled and control pups. A total of 156 miRNAs were expressed at levels above background, with five showing significant (P < 0.05) differential expression between AMC and control samples. Of these five, transcript levels for rno-miR-488, rno-miR-144, and rno-miR-542-5p were increased while transcript levels for rno-miR-421 and rno-miR-376b-5p were reduced (Figure 4 A). The five significant (P < 0.05) DE miRNAs were then analyzed alongside the significant (q < 0.05) DE mRNAs for both a bioinformatically predicted match and reciprocal DE relationship in order to produce a list of 127 putative target genes (Supplemental Table 9). Among the 127 putative mRNA targets of DE miRNAs are the stress response genes Gal (rno-miR-144), Pomc (rno-miR-488), and Pvalb (rno-miR-488) [38–40]. Other miRNA targets include the sialidase Neu1 (rno-miR-542-5p) and the Amyloid beta precursor protein (APP) protease Adam10 (rno-miR-421), and transcripts for the neuronal regulators Pcdh9 (rno-miR-376b-5p) and Avpr1a (rno-miR-488) (Supplemental Table 9) [41–43]. Finally, rno-mir-542-5p was predicted to target the Ube3a1 ceRNA and the two showed a reciprocal relationship in expression (Figure 4 B).

Figure 4. Differences in miRNA expression between AMC (n = 7) and control (n = 5) PND9 hypothalami. (A) Hierarchal clustering analysis of significant (P < 0.05) differential miRNA expression between. The expression of each gene was normalized to a mean of 0 and a standard deviation of 1. The values range from -2.22 to 2.22. Downregulated genes with AMC have negative values (blue), genes with no change have a value of zero (grey), and upregulated genes with AMC have a positive value (red). (B) Base-pairing of the predicted relationship between rno-mir-542-5p and the Ube3a1 3’UTR.

Integrative epigenomic prioritization of differentially regulated genes

To identify a convergent set of genes regulated by AMC in hypothalamus we integrated genes identified by the DNA methylation, RNA-seq, and multiplex miRNA assay analyses. For this combined analysis the 127 unique genes potentially regulated by the five differential miRNAs from Supplemental Table 9 were compared with the 4,023 unique genes in proximity to DMRs shown in Figure 2. The overlap between these three gene lists resulted in 20 unique genes (Figure 5 A). The presence of factors from diverse processes suggests that many pathways contribute to the AMC response in hypothalamus (Table 3). For example, Avrp1a functions in the stress pathway and regulates social behaviors, Kcnh5 is an ion channel that regulates hormone and neurotransmitter release, and Ring1 is representative of several factors on this list that repress gene expression levels (Table 3) [44–46]. Ring1, a component of the Polycomb Repressive Complex 1 (PRC1) complex, together with Ube3a, another ubiquitin E3 ligase, has been recently shown to regulate genes associated with autism [23]. The integration of genome-wide RNA-seq, DMRs, and miRNA target data sets produces a refined set of gene transcripts (Figure 5 A), while also demonstrating the genome-wide nature of the differential expression and epigenetic modification that follows AMC (Figure 5 B).

Figure 5. AMC is associated with genome-wide differences in DNA methylation and gene expression. (A) Integration of DMR, RNA-seq, and miRNA analyses yields 20 prioritized genes (inset). (B) Circos plot of the different data sets. The outer heatmap represents the percent differences in methylation of DMRs, the inner heatmap represents the log₂ fold change of DE transcripts, and the links represent miRNAs and their predicted target genes.

Table 3. Integrated DMR and miRNA analysis.

Symbol	Name	FC	P value	miRNA	P value	FC	Region	Size	Distance to TSS	AreaStat
Avpr1a	Arginine vasopressin receptor 1A	−1.7	0.00680	miR-488	0.012	1.2	Intron	366	1949	23.9
Kcnh5	Potassium voltage-gated channel subfamily H member 5	2.0	0.00045	miR-376b-5p	0.036	−1.1	Intron	334	166874	−24.5
Aldh1b1	Aldehyde dehydrogenase 1 family, member B1	−2.3	0.00005	miR-542-5p	0.017	1.4	Intron	287	979	−10.6
Fbxo2	F-box protein 2	−1.6	0.00390	miR-542-5p	0.017	1.4	Intron	320	820	17.9
							Intron	166	2987	−21.2
Stambp	Stam binding protein	−1.6	0.00420	miR-144	0.039	2.0	Intron	407	5190	−47.0
							Intron	145	710	−14.6
Commd7	COMM domain containing 7	−1.6	0.00485	miR-542-5p	0.017	1.4	Upstream	47	−10	−23.8
Ercc2	ERCC excision repair 2, TFIIH core complex helicase subunit	−3.3	0.00010	miR-542-5p	0.017	1.4	Intron	108	997	−14.4
							Intron	183	1784	−11.1
Smpd1	Sphingomyelinphosphodiesterase1	−1.8	0.00060	miR-542-5p	0.017	1.4	Upstream	337	−1180	11.1
				miR-144	0.039	2.0
Fank1	Fibronectin type 3 and ankyrin repeat domains 1	−1.7	0.00685	miR-542-5p	0.017	1.4	Intron	441	13675	10.0
Gal	Galanin and GMAP prepropeptide	−6.4	0.00005	miR-144	0.039	2.0	Intron	454	1385	30.2
Mus81	MUS81 structure-specific endonuclease subunit	−1.6	0.00795	miR-488	0.012	1.2	Intron	447	923	16.2
Dcps	Decapping enzyme, scavenger	−1.7	0.00115	miR-488	0.012	1.2	Intron	532	38640	12.3
Nme6	NME/NM23 nucleoside diphosphate kinase 6	−1.7	0.00345	miR-542-5p	0.017	1.4	Intron	261	1364	−16.3
Abcc9	ATP binding cassette subfamily C member 9	1.7	0.00300	miR-376b-5p	0.036	−1.1	Intron	797	11185	−51.8
							Exon	86	32330	12.4
Olig1	Oligodendrocyte transcription factor 1	−1.6	0.00640	miR-542-5p	0.017	1.4	Upstream	536	−1021	21.7
Bnip1	BCL2/adenovirus E1B interacting protein 1	−1.7	0.00270	miR-542-5p	0.017	1.4	Upstream	201	−562	−11.7
Prpsap1	Phosphoribosyl pyrophosphate synthetase-associated protein 1	−2.0	0.00005	miR-542-5p	0.017	1.4	Intron	325	1445	−20.2
Gmpr	Guanosine monophosphate reductase	−1.7	0.00105	miR-488	0.012	1.2	Intron	199	9206	−20.9
Ring1	Ring finger protein 1	−1.6	0.00725	miR-144	0.039	2.0	Exon	363	1237	−60.9
Wfdc1	WAP four-disulfide core domain 1	2.7	0.00065	miR-376b-5p	0.036	−1.1	5′UTR	332	211	−12.4

Discussion

Early life experience has tremendous power to shape neurodevelopment with long-term consequences [1–3]. In humans, disjointed or reduced maternal care that occurs in cases of poverty, maternal depression, or addiction is associated with increased risk for developing cognitive and emotional deficits, neuropsychiatric disorders, and cognitive decline later in life [47,48]. On the opposite end of the spectrum, studies in humans suggest that a positive early-life environment can produce resilience to affective disorders [49–51]. In animal models, enhancements in maternal care can positively impact the developing brain leading to reduced anxiety, resilience to depression-like behaviors, and enhanced learning and memory [10,11,15]. In the AMC rat paradigm, while the brief removal of the pups is a stressful event for both the dam and pups, it triggers increased maternal licking and grooming as well as nursing of the pups. Together, these actions have been shown to re-program neonatal HPA circuits (14). The programming of the stress response in offspring receiving augmented maternal care can last a lifetime [52] and help offset age-related cognitive deficits [14]. Much of the previous work on AMC has focused on regulation of stress response genes such as Crh and the Nr3c1. While critical for the AMC phenotype [11] these two components of the stress system are likely not the only important factors for such large-scale and diverse reprogramming of the developing brain.

Together, our results provide the first integrative, genome-wide overview of the genes and pathways that contribute to the “reprogrammed” cognitive and emotional phenotype resulting from AMC. The results offer five novel observations. First, we identified a significant (P = 0.015) brain region specific, reduction of global DNA methylation in the hypothalamus, as well as five large-scale PMDs. Second, we identified 9,439 significant (P < 0.05) small-scale DMRs associated with 4,023 unique genes, which function in neurodevelopment and neurotransmission and may contribute to the establishment and maintenance of the behavioral phenotype. Third, we identified 2,464 significant (q < 0.05) DE transcripts by RNA-seq analysis that extends our understanding of the gene expression differences to neurotransmission and other crucial neurologic functions by highlighting a role for translation and oxidative phosphorylation. Fourth, we identified the predicted targets of 5 significantly (P < 0.05) altered miRNAs and also a putative interaction between mir-542-5p and the Ube3a1 3’UTR. Finally, we integrated these diverse data sets and identified a list of 20 genes with three lines of evidence (methylation, expression, and altered miRNA target) that are highly relevant to the stress resiliency phenotype resulting from AMC.

Consistent with prior descriptions for reduced methylation together with reduced CRH expression in the hypothalamus of AMC rats [21], the connection between genome-wide DNA methylation and gene expression was complex, and more apparent at the level of pathways than at individual genes. At the level of gene ontology, the cellular component term “synapse part” is significantly (q < 0.05) enriched in both datasets. Reactive oxygen species (ROS) generated by oxidative phosphorylation in the mitochondria are not only a key signaling molecule behind synaptic plasticity and memory formation [53], but are also critical to neurodevelopment since they regulate PI3K/AKT/mTOR/PTEN signaling [54]. Furthermore, the hypomethylated Grb10 gold DMR (FWER <0.05) belongs to a genomically-imprinted gene encoding for a receptor-binding protein that interacts with insulin and insulin-like growth factor receptors of the PI3K/AKT/mTOR/PTEN signaling pathway [55]. Although a CpG island in the promoter of Grb10 is hypomethylated in AMC, RNA-seq analysis does not show altered transcript abundance. This disconnect could reflect an ongoing regulatory process where transcriptional changes are not yet or no longer apparent or may be a limitation of our analysis. Translational control is achieved via eukaryotic initiation factors (eIFs) that interact with the ribosome to drive the protein synthesis needed for synaptic plasticity [56]. Additionally, PI3K/AKT/mTOR/PTEN signaling, Wnt/β-catenin signaling, and protein kinase A (PKA) signaling converge to shape neuroplasticity, neurogenesis, cell survival, and resilience [57]. Finally, protein kinase A signaling is activated by G-protein-coupled receptors (GPCRs) that have neurotransmitter and hormone ligands, which include CRH and URN3 [58]. The GPCRs of PKA signaling utilize cyclic adenosine monophosphate (cAMP) as a secondary messenger to activate PKA, which is also affected by handling to alter relevant gene expression profiles, potentially including the glucocorticoid receptor gene Nr3c1 [52]. Ultimately, the overlapping pathways between the methylation and gene expression suggest a broad and potentially long-term impact on numerous cellular and molecular functions important for hypothalamic function.

The genes and pathways identified from our unbiased analysis of AMC-induced hypothalamic changes may be used towards the development of biomarkers and therapeutics. Suderman et al. have shown that epigenetic programming of the HPA axis by maternal care is conserved in humans at a select locus [59]. Given the potential conservation, the suite of new candidate loci identified may be examined in humans towards the development of novel epigenetic biomarkers of resiliency [60]. Such a diagnostic may then be used to identify individuals that would benefit from corrective measures, which include environmental enrichment [61] and pharmacologic agents. Other studies have shown that the transcriptomic profile and epigenetic modifications of altered maternal care can be reprogrammed at adulthood by the global epigenetic modifiers, specifically the methyl donor L-methionine or the histone deacetylase inhibitor trichostatin A (TSA) [18,19,62,63]. More recently, it has been shown that pharmacologic agents can mimic some of the effects of AMC in vitro [12]. Therefore, some of the epigenomic profile established by maternal care is plastic and can be manipulated later in life. However, these previously examined modifiers result in attenuations that are global in nature and may have off-target effects. Our results could be used to overcome this limitation since the suite of therapeutic candidates identified can be utilized as prioritized targets for somatic epigenetic editing studies [64–67]. Ultimately, our results could form the molecular foundation for developing effective interventions for deficient early maternal care in humans or be developed as new treatments for psychiatric disorders [67].

Author contributions

A.V.C., B.I.L., D.S.S., T.Z.B., J.M.L., and D.H.Y. designed research; A.V.C., B.I.L., T.S.T., J.F., A.S-T, and D.H.Y. performed research; K.W.D. and C.E.M. contributed new reagents/analytic tools; A.V.C., B.I.L., K.W.D., R.L.C., and D.H.Y. analyzed data; and A.V.C., B.I.L., and D.H.Y wrote the paper.

Disclosure of Potential Conflicts of Interest

The authors declare that no financial or scientific conflicts of interest exist.

Acknowledgments

This work used the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley, supported by NIH [S10 OD018174] Instrumentation Grant.

Additional information

Funding

National Institute of Mental Health [grant numbers MH73136 and P50 096889 to T.Z.B and A.S-T]; National Institutes of Health [grant number T32MH073124-06 to A.V.C.]; National Institute of Health [grant number 1N01NS076263 to J.M.L.]; Canadian Institutes of Health Research (CIHR) [grant number MFE-146824 to B.I.L].