Perinatal protein malnutrition results in genome-wide disruptions of 5-hydroxymethylcytosine at regions that can be restored to control levels by an enriched environment

ABSTRACT

Maternal malnutrition remains one of the major adversities affecting brain development and long-term mental health outcomes, increasing the risk to develop anxiety and depressive disorders. We have previously shown that malnutrition-induced anxiety-like behaviours can be rescued by a social and sensory stimulation (enriched environment) in male mice. Here, we expand these findings to adult female mice and profiled genome-wide ventral hippocampal 5hmC levels related to malnutrition-induced anxiety-like behaviours and their rescue by an enriched environment. This approach revealed 508 differentially hydroxymethylated genes associated with protein malnutrition and that several genes (N = 34) exhibited a restored 5hmC abundance to control levels following exposure to an enriched environment, including genes involved in neuronal functions like dendrite outgrowth, axon guidance, and maintenance of neuronal circuits (e.g. Fltr3, Itsn1, Lman1, Lsamp, Nav, and Ror1) and epigenetic mechanisms (e.g. Hdac9 and Dicer1). Sequence motif predictions indicated that 5hmC may be modulating the binding of transcription factors for several of these transcripts, suggesting a regulatory role for 5hmC in response to perinatal malnutrition and exposure to an enriched environment. Together, these findings establish a role for 5hmC in early-life malnutrition and reveal genes linked to malnutrition-induced anxious behaviours that are mitigated by an enriched environment.

KEYWORDS:

• 5-hydroxymethylcytosine
• perinatal Malnutrition
• enriched Environment
• anxiety
• ventral Hippocampus
• epigenetics
• neuronal Development

Introduction

Early life adversities such as malnutrition, maternal neglect and abuse, prenatal alcohol exposure, and traumatic stress exposure can have long-lasting effects on physical and emotional health, resulting in a major burden to the families of affected individuals and the healthcare system. The World Health Organization estimates that 45% of deaths among children under 5 years of age are related to malnutrition. Studies in rodents have shown that early-life malnutrition leads to cognitive and learning deficits [1], affects social interactions [2], and increases anxiety and depression-like behaviours [3,4]. The negative impact of perinatal malnutrition is associated with structural, functional, and molecular changes in the brain development, resulting in altered dendritic arborization, myelination, and neuronal plasticity [5], and increasing the likelihood to develop cognitive and behavioural-like disorders [4].

Epigenetic mechanisms such as DNA methylation have emerged as environmentally sensitive biochemical modifications that modulate gene expression [6]. Alterations in DNA methylation levels are linked to a number of neurodevelopmental disorders, including learning disabilities, attention-deficit/hyperactivity disorder, autism, as well as anxiety and depressive disorders [7–11]. Studies in rodents have shown that DNA methylation is nutrition-sensitive, resulting in altered DNA methylation levels on neuropsychiatric-related genes, which suggests epigenetic regulatory mechanisms have a key role in the development of early life adversity-related behaviours [12,13].

While studies in DNA methylation have focused heavily on 5-methylcytosine (5mC), oxidation of 5mC by the ten-eleven translocation (TET) family of enzymes has found stable molecular derivatives, including 5-hydroxymethylcytosine (5hmC) [14]. 5hmC may have central nervous system (CNS)-specific functions, as it is enriched more than 10-fold in the CNS compared to peripheral tissues, associated with the regulation of neuronal activity, and accumulates in the brain during neuronal development and maturation [15,16]. Moreover, perturbations in 5hmC levels are linked to several developmental disorders (e.g. Rett syndrome, autism, and mood disorders) and neurodegenerative diseases (e.g. Alzheimer’s, Parkinson’s, and Huntington diseases) [17,18], suggesting a role for 5hmC in the establishment, maintenance, and survival of the CNS [19–21]. Thus, disruptions in 5hmC levels associated with the development of cognitive and behavioural disorders are of high interest.

The molecular impact of a social and sensory stimulation (enriched environment) on the development of the brain has been deeply studied in rodents: enriched environment improves spatial learning [22] and memory [23], and ameliorates anxiety-like behaviour[24]. An enriched environment also enhances neuronal plasticity and attenuates cognitive deficits caused by brain lesions, such as cortical or hippocampal injuries [25], degeneration of cerebellar cortex [26], prenatal alcohol exposure [27], and anoxia [28]. Recently it was shown that an enriched environment modulates 5hmC levels in the hippocampus and improves learning and memory in aged rodents [29], suggesting a role for 5hmC in neuronal plasticity induced by an enriched environment. Since the ventral region of the hippocampus modulates neuronal processes associated with anxiety-like behaviour[30], the present study examined the genome-wide 5hmC distribution in the ventral hippocampus of adult female mice that were perinatally exposed to a low protein diet that resulted in anxiety-like behaviours and found disruptions of 5hmC across the genome. In addition, post-weaning social and sensory stimulation of the malnourished mice alleviated the anxiety-like behaviours and restored 5hmC levels to control levels in 34 genes, revealing key genes that might be playing a role in the anxiety-like phenotype. Together, these findings provide new insight into the dynamic life-long sensitivity of 5hmC and its influence on neuronal circuits that modulate adult behaviours.

Results

Enriched environment attenuates the anxiety-like behaviour induced by perinatal protein malnutrition

To determine the impact of early-life malnutrition on adult female behaviour, mice were subjected to a low protein diet (8% casein) or standard protein diet (20% casein) during gestation and lactation periods (Supplementary Fig. 1). Offspring were continuously monitored for weight and as we previously published, at post-natal day (P) 20 offspring on the low protein diet weighed significantly less than offspring on the standard diet (Figure 1a; One-way ANOVA with repeated measure showed effects of diet (F = 18.62, P-value = 0.0002), time (F = 574.0, P-value <0.0001) and interaction between diet and time (F = 19.54, P-value = 0.0001) [3,31]. At weaning, offspring from both diets were randomly assigned to either a standard or enriched environment. Subsequent analyses were conducted on four distinct groups: effect of perinatal low protein: standard protein diet and standard environment (CTL) versus low protein diet and standard environment (LP) and effect of enriched environment on perinatal low protein: perinatal low protein diet and enriched environment (LP+EE) versus low protein diet and standard environment (LP). Mice were weighed every 3–4 days and by P42 mice in the enriched environment groups (EE and LP+EE) weighed significantly less than controls (Control and LP, respectively), suggesting that increased stimulation results in weight loss. Indeed, malnourished mice (LP and LP+EE) showed decreased body weight throughout the experiment compared to controls (Control and EE, respectively; Figure 1b; Two-way nested ANOVA with a repeated measure identified significant effect of diet (F = 54.49, P-value <0.0001), environment (F = 15.39, P-value = 0.0002), time (F = 1388.73, P-value <0.0001), the interactions between diet and time (F = 8.15, P-value <0.0001), environment and time (F = 11.88, P-value <0.0001), and diet and environment (F = 32.53, P-value <0.0001)).

Figure 1. Results of perinatal malnutrition and enriched environment in anxiety-like behaviour in females. (a) Body weights of Control (CTL, white) and perinatal malnutrition (LP, black) female mice during lactation at postnatal day 2 and 20 (x-axis; N_CTL = 12, N_LP = 17). (b) Weights of CTL (white dots), enriched environment (EE, light grey dots), LP (black triangles), and perinatal low protein diet and enriched environment (LP+EE, dark grey triangles) female groups measured at several timepoints (x-axis). Significant P-value < 0.05 is depicted as * (LP vs. CTL and LP+EE vs. EE), # (EE vs. CTL), and + (LP+EE vs. LP). (c) Time spent in open arms (seconds; y-axis) in the elevated plus maze (EPM) test. (d) Time spent at the far end of open arms in the EPM test. (e) Square root of time spent in the centre in the open field (OF) test. (f) Square root of latency to enter the centre of the OF. Orthogonal contrast performed to compare treatments (N_CTL = 11–12, N_EE = 12, N_LP = 8–10, N_LP+EE = 7–8). All values are expressed as averages ± SEM and significance is indicated as *P-value < 0.05, **P-value < 0.01, and ***P-value < 0.001, unless otherwise noted above

At P56, mice were randomly assigned to behavioural tests (N_CTL = 12, N_EE = 12, N_LP = 10, N_LP+EE = 8) or sacrificed for molecular tests (N = 6 per experimental group) and female offspring were examined for anxiety-like behaviours using the elevated plus maze and open field tests. Elevated Plus Maze: The LP group spent significantly less time in the open arms than the Control group (Two-way ANOVA for diet (F = 1.44, P-value = 0.2377), environment (F = 24.44, P-value <0.0001), and interaction between diet and environment (F = 3.31, P-value = 0.0769); LP = 40.2 ± 8.8; Control = 75.0 ± 12.5; P-value = 0.035; Figure 1c), suggesting perinatal malnutrition results in anxiety-like behaviour in female mice. In contrast, the EE group as well as the LP+EE group spent significantly more time in the open arm compared to the Control and LP groups, respectively (EE = 109.9 ± 11.6, P-value = 0.028 and LP+EE = 117.1 ± 11.6, P-value = 0.0001; Figure 1c), suggesting the enriched environment mitigates the anxiety-like behaviour. Additionally, the LP+EE group spent significantly more time at the far end of the open arms compared to LP group (Two-way ANOVA for diet (F = 0.07, P-value = 0.7902), environment (F = 6.75, P-value = 0.0133), and interaction between diet and environment (F = 1.98, P-value = 0.1672); LP = 6.5 ± 3.8; LP+EE = 30.0 ± 7.3; P-value = 0.012), unlike compared to the Control and EE groups (Control = 13.2 ± 5.2; EE = 20.2 ± 6.5; P-value >0.05, Figure 1d).

Open Field: The LP group spent significantly less time in the centre of the open field compared to the Control group (Two-way ANOVA for diet (F = 3.35, P-value = 0.0758), environment (F = 13.21, P-value = 0.0009), and interaction between diet and environment (F = 2.20, P-value = 0.1472); LP = 8.6 ± 0.9; Control = 11.4 ± 0.7; P-value = 0.022; Figure 1e), suggesting perinatal malnutrition results in anxiety-like behaviour in female mice. In contrast, LP+EE female mice spent significantly more time in the centre of the open field compared to the LP mice (LP+EE = 12.9 ± 0.8; P-value = 0.0025; Figure 1e). In addition, the time spent to enter the inner centre of the open field was significantly greater in LP mice compared to Control and LP+EE mice (Two-way ANOVA for diet (F = 3.06, P-value = 0.08), environment (F = 10.25, P-value = 0.002), and interaction with diet and environment (F = 1.61, P-value = 0.2); LP = 12.9 ± 2.5; Control = 7.8 ± 1.7; P-value = 0.039; LP+EE = 5.3 ± 1.2; P-value = 0.0066; figure 1f), suggesting the enriched environment rescued the anxiety-like behaviour. The total distance travelled among each group was not significantly different (Control = 88.5 ± 5.5; EE = 99.9 ± 5.5; LP = 79.0 ± 8.7; LP+EE = 93.9 ± 8.5; P -values >0.05; Supplementary Fig. 2). Together, these results showed perinatal protein malnutrition results in anxiety-like behaviour in adult female mice that can be attenuated by post-weaning exposure to an enriched environment.

Long-lasting disruptions of 5hmC levels following perinatal protein malnutrition and enriched environment

To determine the effect of perinatal protein malnutrition and an enriched environment on the genome-wide distribution of 5hmC, we employed an established chemical labelling and affinity purification method, coupled with high-throughput sequencing technology [15,32,33]. A pool of two age-matched Control, EE, LP, and LP+EE female offspring (N = 3 pools for each group) were sacrificed at P56 as independent biological replicates. 5hmC containing DNA sequences were enriched from ventral hippocampal total DNA and high-throughput sequencing resulted in a range of ~20–30 million uniquely mapped reads from each biological replicate (Supplementary Table 1). Read density mapping had high correlations between all mice and no significant differences in the correlations among the mice from the same group and between mice from different groups, suggesting that differences between groups are small and only in a short regions of the genome. Moreover, these data provided no visible differences among the chromosomes, except depletion on chromosome X, which is consistent with previously reports (Supplementary Fig. 3A) [16,33]. Annotation of 5hmC reads to standard overlapping genomic features obtained from the UCSC Table for NCBI37v1/mm9, showed that 5hmC mapped reads were found in all regions with enrichment in intergenic regions compared with genic structures (Supplementary Fig. 3B). Nonetheless, these distributions demonstrated an overall equal distribution of 5hmC levels in all groups on these defined genomic features.

To identify distinct 5hmC distribution patterns throughout the genome of female mice following perinatal malnutrition with and without exposure to an enriched environment, we characterized differentially hydroxymethylated regions (DhMRs) between: 1) LP and Control groups, to determine the long-lasting effects of perinatal protein malnutrition on 5hmC levels; as well as 2) LP and LP+EE groups, to investigate the dynamic effects of an enriched environment on perinatal protein malnutrition-associated 5hmC levels. A total of 274 significant LP increases in hydroxymethylation (LP-specific 5hmC) and 561 significant LP decreases in hydroxymethylation (Control-specific 5hmC) were found and these loci were distributed across all chromosomes (P-value <0.05 and Fold-change >1.5; Figure 2a; Dataset 1). We proceeded to define the genomic features associated with perinatal low protein diet by annotating the DhMRs to overlapping gene structures or to the intergenic region if more than 10 kilobases away from any gene structure (Figure 2b). To determine the relative abundance of DhMRs in each annotated region we compared, via binomial testing, the proportion of DhMRs to the total number of 5hmC peaks detected for each region (Materials and Methods). The distribution of all DhMRs indicated a significant under-representation of DhMRs on the central promoter (1–2 kilobases upstream of start) and 5ʹ UTR regions and a significant over-representation of DhMRs on other intron (introns 2+) and distal intergenic regions, suggesting a role for 5hmC in gene regulation. The LP-specific DhMRs were significantly over-represented within exonic regions and significantly under-represented within intronic regions, while Control-specific DhMRs were significantly over-represented within intronic and distal intergenic regions and significantly under-represented within the central and proximal promoter, 3ʹUTR, and exonic regions (Figure 2b). Together, these data suggest that perinatal protein malnutrition alone results in stable non-random changes in 5hmC levels that can last for a lifetime.

Figure 2. Characterization of 5hmC profiles across the genome. (a and c) Modified Manhattan plots depict 5hmC peaks from female ventral hippocampus are distributed across the entire genome that were assessed for effects of (a) perinatal protein malnutrition- (c) perinatal protein malnutrition+enriched environment mice. Chromosomes across the genome alternate between black and grey. The significance (y-axis, – log10(P-value)) is shown for each dot line, which represents a region examined for differential hydroxymethylation. Solid lines illustrate a P-value < 0.05. LP- and LP+EE-specific 5hmC (N = 274 and 712, respectively) were found above the upper solid line while Control- and LP-specific 5hmC (N = 561 and 349, respectively) are displayed below the lower solid line. (b and d) Bar plots showing the percentage (y-axis) of all DhMRs (grey), Control- (white) and LP- (black) specific DhMRs and LP- (white) and LP+EE- (black) specific DhMRs compared to the genomic background of all 5hmC peaks (hatched) on various genomic structures (x-axis) in comparisons of (b) perinatal protein malnutrition- (d) perinatal protein malnutrition+enriched environment mice. Significant over- and under-representation of DhMRs compared to all 5hmC peaks were generated by permutation testing and are indicated with an asterisk. *P-value < 0.05

An enriched environment on protein malnutrition-associated 5hmC levels resulted in a total of 712 LP+EE increases in hydroxymethylation (LP+EE-specific 5hmC) and 349 LP+EE decreases in hydroxymethylation (LP’-specific 5hmc) were found and these loci were distributed across all chromosomes (P-value <0.05 and Fold-change >1.5; Figure 2c; Dataset 2). The distribution across genomic features of all DhMRs (Figure 2d) indicated a redistribution of DhMRs compared to those found by LP alone (Figure 2b), by exhibiting a significant under-representation of DhMRs on the distal intergenic and distal promoter regions and a significant over-representation on the first intron and 3ʹ UTR regions (P-value <0.05). The LP+EE-specific DhMRs were significantly over-represented on the proximal promoter (<1 kilobases upstream of start), first intron, other exon (exons 2+), and 3ʹ UTR regions and significantly under-represented on the distal promoter (2–3 kilobases upstream of start) and distal intergenic regions. The LP’-specific DhMRs were significantly over-represented on distal intergenic regions and significantly under-represented in the central promoter and 3ʹUTR (P-value <0.05; Figure 2d). Together, these data indicate that malnutrition+enriched environment-associated DhMRs are unique to malnutrition alone DhMRs and are not randomly distributed throughout the genome.

Annotation of malnutrition- and malnutrition+enriched environment-associated DhMRs to genes reveals genes related to neuronal development, anxiety, and stress response

To first assess the long-lasting molecular effects of perinatal low protein diet on female mice, we annotated the malnutrition-associated DhMRs to genes and found 170 and 338 genes contained either LP- or Control-specific 5hmC, respectively (Dataset 1). To examine the biological relevance of the malnutrition-related genes, we performed gene ontological (GO) analysis of the LP- and Control-specific 5hmC-associated genes (N = 508) and found a significant enrichment of terms linked to central nervous system development (90% of the first 20 terms), modification of postsynaptic actin cytoskeleton (P-value <0.05), as well as cortisol synthesis and secretion (P-value <0.05; Figure 3a; Dataset 3). These findings suggest that a perinatal low protein diet leads to long-term modulation of 5hmC in genes related to neuronal development and stress response.

Figure 3. The top 20 significant results from gene ontological analyses. (a and b) Heatmaps of the top 20 significant gene ontology terms from perinatal protein malnutrition-induced DhMR-associated genes (a) and perinatal protein malnutrition+ enriched environment-induced DhMR-associated genes (b) (y-axis). Heatmap colour gradient indicates the percent of DhMR-associated genes contributing to the input list size relating to each term, and range from 0 genes (yellow) to >80 genes (red)

To examine the dynamic molecular effects of perinatal malnutrition followed by an enriched environment on female mice, we annotated the malnutrition+enriched environment-associated DhMRs to genes and found 540 and 203 genes contained either LP+EE- or LP’-specific 5hmC, respectively (Dataset 2). Annotation of the LP+EE- and LP’-specific 5hmC-associated genes (N= 743) to GO pathways found a significant enrichment of terms linked to central nervous system development such as neuron projection development (P-value = 1.20E-7), neuron differentiation (P-value = 9.94E-6) and axonogenesis (P-value = 2.87E-4; Figure 3b; Dataset 4), suggesting that perinatal malnutrition followed by an enriched post-weaning environment affects neuronal plasticity, which may contribute to the attenuation of the anxiety-like behavioural phenotype. Finally, we independently investigated the overlap of both the malnutrition- and malnutrition+enriched environment DhMR-associated genes to genes know to be involved in anxiety and neuronal development and found a significant overlap with the malnutrition+enriched environment-DhMR-associated genes with anxiety (N = 114 of 704; Figure 4a) and neuronal development (N = 286 of 704; Figure 4b). This finding was not observed with the perinatal malnutrition-DhMR-associated genes, suggesting the enriched environment has a distinct effect on neuronal developmentally important genes and may target them to ameliorate the anxiety-like phenotype.

Figure 4. Comparison of DhMR-associated genes to genes related with anxiety and CNS development functions. (a) Venn diagrams display the overlap of DhMRs-associated genes with known anxiety-related genes, as identified in the LP vs. Control (left panel) and LP+EE vs. LP (right panel) comparisons. (b) Venn diagrams displaying the overlap of DhMRs-associated genes with known CNS development genes, as identified in the LP vs. Control (left panel) and LP+EE vs. LP (right panel) comparisons. Significant overlaps were identified by Chi-square test and are indicated by asterisks (*p < 0.05 **p < 0.01)

Neurodevelopmental disorders can involve the disruption of transcription factor function [34,35]; thus, we investigated the potential for the perinatal malnutrition to alter DNA-binding of transcription factors by testing for enrichments of known transcription factor sequence motifs among the DhMRs (Material and Methods). Several transcription factor sequence motifs were enriched in the regions containing CTL-, LP-, LP+EE- and LP’-specific 5hmC. Notably, many of the transcription factors that bind these motifs have links to neurodevelopmental disorders, such as Lrf, Hif-1a, and Arnt (Supplementary Fig. 4) [36–38]. Together, these findings suggest that 5hmC influences transcription factor binding following perinatal malnutrition and enriched environment.

Enriched environment restores the aberrant 5hmC levels of 34 genes in the perinatal malnutrition mice to control levels

To identify genes contributing to the rescued anxiety-like behaviour in the LP+EE mice, we examined the DhMRs found in the same gene of both LP versus Control and LP+EE versus LP comparisons but exhibited an opposite effect on 5hmC abundance between groups. This approach revealed 34 DhMR-associated genes that exhibited a non-significant difference in 5hmC levels between Control versus LP+EE groups (i.e. the LP+EE 5hmC levels were restored to control levels; Table 1). The relationship between these molecular data and the rescue of the anxiety-like behaviour can be summarized through the known synaptic proteins, receptors, and signalling pathways associated with anxiety in humans (Figure 5). The enriched environment restores altered 5hmC levels to control levels on neurodevelopmental-related genes and pathways related to anxiety, including genes with a known role in synapse organization, neurite growth, and axon guidance (e.g. Fltr3, Itsn1, Lman1, Lsamp, Nav, and Ror1) and epigenetic processes (e.g. Hdac9 and Dicer1). Together, these data indicate that the perinatal malnutrition causes dynamic epigenetic disruptions in orthologous genes with common roles in anxiety-related pathways.

Table 1. Genes with 5hmC levels restored by environmental enrichment

Decreased Armc9, Cherp, Dicer1, Fabp6, Itsn1, Nav1, Pitpnb, Zc3h12a
Increased1810013L24Rik, Adgrb3, Clip4, E030024N20Rik, Flrt3, Gm3985, Hdac9, Heatr3, Lman1, Lsamp, Me3, Ppp2r2b, Psd3, Rab30, Rheb, Rnf170, Rngtt, Ror1, Rpp25l, Sfi1, Sox6, Stard4,

Figure 5. A working model of molecular pathways effected by perinatal protein malnutrition and the genes exhibit a restoration of 5hmC levels to control level in the LP+EE mice. Schematic representation of genes/pathways involved in mechanisms that intervene in neuronal growth mediated by the dynamic stability of microtubules, maintenance of neural circuits by cell-adhesion proteins, trafficking and transport of molecules, and epigenetics regulations. Genes marked with a red circle indicate LP-specific 5hmC, while green circles indicate LP+EE-specific 5hmC

Discussion

Here, we show that perinatal malnutrition-induced anxiety-like behaviours in adult female mice and that post-weaning exposure to a social and sensory stimulation (enriched environment) rescued this behaviour. Several studies in both humans and rodents have demonstrated that maternal malnutrition affects neurocognitive and neuropsychological development [39] and increases the risk to develop behavioural disorders [24,33,40]. In contrast, the combination of enriched environment improves both behavioural and molecular components: reduces anxiety-like responses [41], enhances dendritic branching and spine density in the hippocampus [42], increases the cortical and whole-brain weight [43], and stimulates hippocampal neurogenesis [44]. Thus, since neurogenesis in ventral dentate gyrus specifically plays a role in emotional and affective disorders [45], these previous findings may suggest that an enriched environment promotes neurogenesis in the ventral hippocampus to mitigate the protein malnutrition-induced anxiety-like behaviour observed here in adult female mice.

Early-life stressors that cause anxiety-like behaviours are linked to stable changes in 5hmC levels throughout the genome [33]. In addition, exposure to an enriched environment alters epigenetic mechanisms, including 5hmC abundance [29]. Together, these data led to the hypothesis that exposure to an enriched environment may reverse 5hmC modifications caused by perinatal malnutrition. Consistent with this hypothesis, we found that perinatal malnutrition results in more decreases in 5hmC levels, while the enriched environment results in more increases in 5hmC levels (Figure 2). Notably, it was previously shown that an enriched environment decreases genome-wide 5hmC levels in the hippocampus [29]; however, this study examined the whole hippocampus and utilized a different enrichment model, based on aged mice (18-month-old) and limited enrichment exposure (3 hours every day during 4 week period). Thus, it is difficult to compare the findings of these two studies, other than to acknowledge that they both find dynamic changes in 5hmC levels that are specific to brain regions and enriched environmental paradigms.

Bioinformatic analysis revealed that a significant number of known anxiety- and CNS development-related genes harbour alterations in 5hmC levels following exposure to an enriched environment. Since these significant findings were not observed with perinatal malnutrition alone, they likely indicate genes targeted for the rescue of the anxiety-like behaviour. These results are consistent with previous findings of 5hmC linked to neuronal development, particularly on genes involved in axon guidance [29] and neuronal functions, such as synaptic transmission and neuronal activity in young adult mice [16]. Mechanistically, it is intriguing that the binding motifs of several transcription factors with known roles in neuronal development were found among the DhMRs, including roles in cellular homoeostasis and hypoxia. It is notable that the enriched transcription factor sequence motifs were distinct between malnutrition- and malnutrition+enriched environment-associated DhMRs. Together, these findings suggest that 5hmC could be involved in the regulation of genes linked to early-life malnutrition-induced anxious behaviours and their rescue following exposure to an enriched environment. Further studies are needed to definitively determine the role of differential hydroxymethylation within transcription factor binding sites in stress-induced developmental brain disorders.

Enriched environment restoration of aberrant 5hmC levels is found on genes related to microtubule, cell adhesion, and Golgi complex functions

It was compelling to find 5hmC levels restored to control levels on 34 genes in the anxiety-like behaviour rescued LP+EE mice. Since most genes related to neurodevelopmental phenotype are not expressed at P56 (when these mice were sacrificed) it will be imperative to examine gene expression at multiple timepoints throughout development to find correlated changes in 5hmC, gene expression, and the observed anxiety-like phenotype, such as at birth, P21, and P56. Many of these 34 genes have known neuronal roles related to microtubule, cell adhesion, and Golgi complex functions, and represent neuroplasticity genes that could be targeted to form the basis of early interventions and therapies. Below is a summary of the possible relevance of these genes to the observed anxiety-like behaviours, which also are depicted in Figure 5.

Neuronal outgrowth and survival require the participation of microtubules [46]. Several genes affecting microtubules exhibited restored 5hmC levels in the LP+EE mice. For example, Ror1 encodes a receptor tyrosine kinase which modulates neuronal growth and branching patterns in hippocampal neurons and impacts the expression levels of microtubules-associate proteins like MAP2 and MAP1B [47]. The gene encoding the signalling scaffold intersectin-1 protein (Itsn1) is a component of the Reelin pathway. ITSN1 is associated with very-low-density lipoprotein receptor (VLDLR) and adaptor protein Disabled 1 (DAB1) to facilitate Reelin signalling[48] that regulates microtubule dynamics during neural migration and synaptic plasticity in the hippocampus [49]. Finally, Nav1 encodes neuron navigator 1, a microtubule-associated protein that plays an important role in neuronal migration [50]. Nav1 is a known stress-related gene and its DNA methylation levels also are altered in human anxiety [51]. Together, these data suggest that malnutrition could impact microtubular stability and the dendritic plasticity of neurons, and that exposure to an enriched environment might result in an elongation of dendrites, which likely contributes to the rescue of the impaired behaviour.

Another group of genes that exhibited restored 5hmC levels in the LP+EE mice encoding cell-adhesion proteins, which are essential for processes related to neuronal tissue organization such as cell migration, neurite outgrowth, and maintenance of the neuronal circuits [52]. For example, Fltr3 encodes a member of the fibronectin leucine-rich transmembrane protein family, an adhesion molecule that promotes axon outgrowth by contacting neurons following axonal injury [53]. The FLTR3 protein plays a transcellular adhesion role via its interaction with ADGRL3 [54]. Adhesion G protein-coupled receptor is expressed at post-synaptic terminals, promoting glutamatergic synapse development, neuronal migration and is linked with multiple behavioural domains related to attention-deficit/hyperactivity disorder [55]. Lsamp encodes a neural cell adhesion molecules that regulates neural projection development and neurite formation, and is linked to anxiety-like behaviours [56,57]. LSAMP also plays an important role in the ability to adapt to novel stimuli and environment [58], which is relevant for the model system studied here. Adgrb3 encodes brain-specific angiogenesis inhibitor 3 in neurons and is enriched in the postsynaptic density with a role as a regulator of synaptic plasticity and synaptogenesis [59]. Together, alterations of 5hmC levels on these genes support that perinatal malnutrition could be affecting the formation of neuronal circuits and impacts neurotransmitter systems [60] through neural adhesion proteins.

The restoration of 5hmC levels also was found on several genes involved in the integrity of Golgi complex; such as Piptnb, Rab30 and Lman1. PIPTNB is essential for correct signal trafficking during neurogenesis [61], RAB30 maintains the Golgi complex [62], and LMAN1 serves as trafficking factor for transmembrane and soluble proteins, including γ-aminobutyric acid type A receptors ([GABA_ARs] the major mediators of synaptic inhibition in the brain) and serotonin receptors (5HT), respectively [63]. In neurons, the Golgi complex is crucial for proper trafficking of ion channels, receptors, and other signalling molecules and also mediates transport of exogenous molecules by retrograde and trans-synaptic pathways [64]. The restoration of 5hmC levels on these Golgi complex-related genes suggests that the enriched environment could modulate the Golgi complex functionality in the LP+EE mice.

Restoration of 5hmC levels was found on genes that function in histone deacetylation (Hdac9) and miRNA biogenesis (Dicer1). Histones deacetylases are known to play a role in the prevalence of anxiety disorders [65] and early-life stress [66]. DICER1 is a target of β-catenin, when coupled this complex mediates anxiolytic effects and increases behavioural flexibility and resilience [67]. These results suggest that 5hmC may modulate histone deacetylation and miRNA production to impact anxiety-like behaviours.

Finally, we previously showed a similar rescue of anxiety-like behaviour in perinatal malnourished male mice that was linked to the expression of three miRNAs (miR-187-3p, miR-369-3p, and miR-132-3p) [24]. None of these miRNAs are closer than 100,000 base pairs from any of the DhMRs, indicating that it is unlikely that DNA hydroxymethylation is directly regulating their expression following perinatal malnutrition. It was demonstrated that miR-369-3p targets LSAMP [68], which suggests a link between miR-369-3p and 5hmC mediated regulatory mechanisms to influence LSAMP expression levels. Moreover, miR-187-3p exhibits a gender-specific expression pattern [31], suggesting that the behavioural rescue observed in both sexes may involve different molecular mechanisms. Further studies are needed to determine the relationship between 5hmC and miRNA in gene regulation underlying malnourished-related behavioural outcomes.

Conclusion

Here, we establish dynamic fluctuations of 5hmC levels on genes related to early malnutrition and its rescue following exposure to an enriched social and sensory environment. In particular, we have identified 34 genes that exhibited a restoration of aberrant 5hmC levels and future studies focused on the function and expression levels of these neuroplasticity genes could form the basis of early interventions with the potential of preventing malnourished individuals to have a heightened risk to develop anxiety disorder.

Materials and Methods

Animals and Experimental Design

CF-1 mice were obtained from the Bioterio Central, Facultad de Ciencias Exactas y Naturales (University of Buenos Aires, Argentina) and housed under uniform conditions in a pathogen-free mouse facility with a 12-h light/dark cycles. Experiments were performed in accordance with local regulations and the National Institutes of Health (NIH) guide of the care and use of laboratory animals (NIH publication 80–23/96) and approved by the Ethical Committee (CICUAL – Protocol N° 0024) of the Facultad de Ciencias Exactas y Naturales (University of Buenos Aires). The special diets employed in this study during mating and gestation were prepared according to the AIN-93 Final Report [69] and purchased from Research Diets Inc. (New Brunswick, NJ). Nutrient composition of standard protein diet (SP) and low-protein (LP) diet were used previously [24,31,70] and SP and LP diets contain 20% and 8% of casein, as source of protein, respectively.

The experimental design was conducted as previously published [24]. Adult male and female mice (F₀) were provided with standard protein diet and their weight were obtained every 2 days. Upon weight stabilization, (F₀) male and female mice were randomly assigned (computer-based method) to either SP or LP groups and dams were maintained on the specific diet until weaning (postnatal day (P) 21). During mating, one male was housed with two nulliparous females for 5 consecutive days and sires in the same cage with dams assigned to the low protein group were fed a low protein diet during the 5 days of mating. Following mating, the two pregnant females remained together until they were individually housed one week before they gave birth. After offspring birth, efforts were made to balance the litter size within a range of 10 to 12 pups per dam on P2 with a 1:1 female: male ratio when possible. Litter with less than 8 pups were removed from the experiment. At weaning (P21), SP- and LP-fed offspring were provided with standard laboratory chow (balanced meal Cooperación, containing 23% of protein) and housed either in a standard environment (SE, 3 or 4 mice in a standard mouse cage (30 × 18.5 × 12 cm) or an enriched environment (EE, 8 to 10 mice in a bigger size cage (64 × 45 × 25 cm) with novel objects (i.e. ramps, boxes, tunnels, toys, plastic houses) that were randomized every 3 days for 5 weeks (P56) (Supplementary Fig. 1). At P56, mice were randomly assigned to behavioural tests or sacrificed for molecular test, and brains were flash frozen for tissue collection and molecular analysis. Subsequent analysis was conducted on four distinct groups: effect of perinatal low protein: standard protein diet and standard environment (CTL) vs low protein diet and standard environment (LP) and effect of enriched environment on perinatal low protein: perinatal low protein diet and enriched environment (LP+EE) vs low protein diet and standard environment (LP).

Behavioural testing

At P56, animals were tested for anxiety-related behaviours as previously published [3,24]. Tests were performed during the light cycle (9 am to 12 pm, lights on at 8 am). With one week between tests, the same group of mice were subjected to the following sequence of behavioural tests: elevated plus maze and open field. An experimenter who was blind to the animal’s group scored video recording for each test.

Elevated plus maze (EPM)

The apparatus consisted of a dark floor maze with a central platform (5 × 5 cm) from which four arms radiated (30 cm in length). Two arms consisted of dark polystyrene walls (15 cm in height, closed arms), and other two arms consisted of a small ledge (0.5 cm in height, open arms). The maze was elevated 50 cm above the floor during testing. Each mouse was individually placed in the centre of the apparatus facing one of the open arms and allowed to explore for 5 min with lights set at 100 luxs. A video camera was placed above the maze, and ANY-Maze™ was used to record and track the animal’s movements. The total distance travelled, the total time spent in each compartment of the maze, and the time spent at the end of the open arms was analysed.

Open Field test (OF)

The apparatus consisted of a dark square arena (45 × 45 cm) with dark walls (30 cm in height). Each mouse was individually placed next to one of the walls and allowed to explore the apparatus for 20 min with lights set at 100 luxs. Using ANY-Maze™ software a central zone (23x23 cm) was virtually drawn. A video camera was placed above the arena, and ANY-Maze™ was used to record and track the animal’s movements. Total distance travelled, the total time spent in the central zone of the arena, and the latency to enter to the centre zone of the arena was analysed.

Statistics for weight analysis and behavioural tests

Mice were weighed every 3 or 4 days during the environmental housing. For analysis of weights during lactation, we used 12 CTL and 17 LP mice, for analysis of weights during environmental exposure, we used 21–24 CTL, 20–21 LP, 18–24 EE, and 18–20 LP+EE mice. For behavioural tests, we used 11–12 CTL, 8–10 LP, 12 EE, and 7–8 LP+EE mice. A one-way ANOVA (from postnatal day 2 to postnatal day 20) and two-way ANOVA (from P21 to P56) with repeated measures was performed for analysing weight parameter in the same individual during a series of days. Tukey post-test was used to compare treatments between them. For behavioural parameters, a two-way ANOVA and orthogonal contrast performed in order to compare treatments. The data which showed a studentized residual below −2.5 or above 2.5 were considered outliers and they were removed from the analysis. When data failed to pass normality testing a square root transformation was applied. In the cases data did not pass the normality test after the transformation, Kruskal Wallis non-parametric test was employed in the raw data. Levene’s test was used to control the homoscedasticity of the variance.

Tissue acquisition and DNA/RNA extractions

The effects of perinatal protein malnutrition and the enriched environment during development on 5-hydroxymethylation was obtained on normal or low protein fed mice exposed to normal or enriched environment. At P56 mice (a pool of two age-matched female offspring, N = 3 pools for each group) were sacrificed and a 3D printed brain matrix [71] was used to slice the brain and extract the ventral hippocampus (from −2 mm to −3.5 mm of Bregma). DNA and RNA were extracted using AllPrep DNA/RNA mini kit.

5hmC enrichment of genomic DNA

Chemical labelling-based 5hmC enrichment was described previously [15,21,32,33]. Briefly, DNA extracted from ventral hippocampus tissue was sonicated to 300–500 bp and incubated for 1 h at 37°C in the following labelling reaction: 1.5 µl of N3-Uridine diphosphate Glucose (2 mM); 1.5 µl β-Glucosyltransferase (β-GT; 60 µM); and 3 µl of 10X β-GT buffer, in a final volume of 30 µl. Next, 5 µl Biotin (1 mM) was added and the reaction was incubated at 37°C for 2 h before capturing on streptavidin-coupled Dynabeads (Invitrogen, 650.01). Enriched DNA was released from the beads during a 2-h incubation at room temperature with 75 µl Dithiothreitol (100 mM) (DTT; Invitrogen, 15,508,013), which was removed using a Bio-Rad column (Bio-Rad, 732–6227). Capture efficiency was approximately 5–7% for each sample.

Library preparation and high-throughput sequencing of genomic DNA

5hmC-enriched libraries were generated using the NEBNext Ultra II DNA Library Prep Kit for Illumina sequencing, according to the manufacturer’s protocol. Briefly, the 5hmC-enriched DNA fragments were purified after the adaptor ligation step using NEBNext Sample Purification Beads AMPure XP beads (Agencourt A63880). An Agilent 2100 BioAnalyzer was used to quantify the amplified library DNA and 20-pM of diluted libraries were used for sequencing. A 50-cycle single-end sequencing reaction, which was performed by the University of Wisconsin Biotechnology Centre. Image processing and sequence extraction were done using the standard Illumina Pipeline.

Sequence alignment of 5hmC data and peak calling

Raw FASTQ files were first inspected for quality using FastQC. Each paired-end mate file per sample was individually inspected. As each sample was run on two lanes for sequencing, upstream-mates were concatenated to one file and downstream-mates were concatenated to a separate file. An average of 30 million paired-end reads was sequenced per sample. The alignment of paired-end mates to the NCBI37v1/mm9 genome was performed using Bowtie v1.2.2 [72]. During alignment, only uniquely mapping reads with no more than two mismatches throughout the entire read and having the least amount of mismatches (-a – best – strata), with a maximum insertion length of 600bp, were reported and used for downstream analysis. An average of 25.1 million (~83%) paired-end reads were successfully aligned to the genome.

The Model-based Analysis of ChIP-Seq 2 (MACS2) algorithm v2.1.2 [73] was used to estimate fragment size, call peaks, and identify peak summits from aligned paired-end reads using the following parameters: paired-end format, effective genome size of 1.87e9, bandwidth of 300bp, an FDR cut-off of 0.05, auto pair model process enabled, local bias computed in a surrounding 1kb window, and a maximum of one duplicate fragment to avoid PCR bias. Summits were extracted for each peak for each sample and extended ±500bp for downstream analysis.

DhMR identification and genomic annotations

Peaks from control and treatment samples were imported to R environment and converted to Granges objects. Peaks were redefined for each group for differential hydroxymethylation analysis as follows: peaks from each sample were pooled together, overlapping peaks were merged to form one region, peaks were extracted for further analysis if at least two samples had such an overlap. Significant accumulation of unique mapped reads represents hydroxymethylated regions throughout the genome. This produced 48,564 candidate regions for the LP vs CTL comparison, 44,407 candidate regions for the LP+EE vs LP comparison, 44,496 candidate regions for the EE vs CTL comparison and 39,005 candidate regions for the LP+EE vs EE comparison. To form candidate regions for differentially hydroxymethylated region (DhMR) analysis, peaks from control and treatment groups were similarly pooled and merged. DhMRs were identified using R package edgeR [74]. Reads from each sample were extended by the estimated fragment size determined by MACS2. Normalization factors, common dispersion, and tagwise dispersion were calculated for each comparison, aimed at ameliorating differences in library size between samples. Regions were termed DhMRs using a P-value cut-off of 0.05 and a log₂ fold-change of 0.58, which produced 835 DhMRs for the LP vs CTL comparison, 1061 DhMRs for the LP+EE vs LP, 411 for the LP+EE vs EE comparison and 1448 DhMRs for EE vs CTL comparison. DhMRs were annotated to the nearest genomic structure and gene using R package ChIPseeker and Bioconductor packages TxDb.Mmusculus.UCSC.mm9.knownGene and org.Mm.eg.db [75]. A± 3,000bp window surrounding the transcription start site was used termed the promoter region. For the analysis of genomic distribution of DhMRs, the gen was fragmented including Distal Promoter (2–3 kb), Central Promoter (1–2 kb), Proximal Promoter (≤1kb), 5ʹ UTR, First exon, First Intron, Other exons, Other introns, 3ʹ UTR, Downstream (≤3 kb) and Intergenic (>3kb). DhMRs were annotated to the nearest genomic structures and permutation test (p-value <0.05) was used to determine the under/over-representation of DhMRs in these genomic features. Note, DhMRs falling greater that 3,000bp upstream of the transcription start site of the nearest gene or greater than 300bp downstream of the nearest gene were deemed ‘Distal Intergenic’ and not associated with any gene.

Heatmap and Gene Ontology analysis

Genes associated with DhMRs of LP vs CTL and LP+EE vs LP comparisons were studied for ontology enrichment through GO Biological process and KEGG pathway network analysis using Cytoscape 3.7.2 (ClueGO plugin). GO terms were defined as statistically different if the p-value corrected by the Bonferroni method was <0.05. The top 20 most significantly enriched terms for LP vs CTL and LP+EE vs LP are presented in Figure 3a and 3b, respectively.

Heatmaps with enriched terms were generated using the percentage of genes corresponding to CTL, LP and LP+EE-specific 5hmC adjusted to the input list size and were visualized using the pheatmap package in R [76]. The top 20 most significative-enriched terms are presented from the lowest P-value to the highest P-value.

Overlap between DhMR-associated genes and Anxiety- CNS development-related genes

In order to analyse if the overlap between DhMR-associated to a unique gene of LP vs CTL (N= 477) and LP+EE vs LP (N = 704) comparison and genes known to play a role in anxiety (LP vs CTL N = 1,571/LP+EE vs LP N = 1,527) or CNS development (LP vs CTL N = 3,944/LP+EE vs LP N = 3,867), a Chi-square test was used. Known-related genes were extracted from Gene Cards database, using the following terms for CNS development: synaptic plasticity; synapse formation; neurite outgrowth and dendritic arborization. The gene universe employed for the Chi-square test were all the genes associated with peaks from sequenced data (LP vs CTL N = 11,404 and LP+EE vs LP N = 11,224). A Chi-square test was also employed when DhMRs-associated genes of LP vs CTL (N= 508) and LP+EE vs LP (N = 743) were compared to each other.

Sequence motif analysis

For motif discovery analysis, the Hypergeometric Optimization of Motif EnRichment (HOMER) suite of tools was utilized [77]. DNA sequences corresponding to DhMR coordinates were obtained from the mm9 genome and compared against background sequences, using the given size of the region. Enriched known motifs of vertebrate transcription factors (N = 428) were determined using binomial testing and an FDR cut-off of 0.05.

Acknowledgments

This work was supported by the Agencia Nacional de Promoción Científica y

Tecnológica, Argentina (Grants PICT-2013-0653 EC and PICT 2016-0218 MC), Universidad de Buenos Aires (Grant UBACYT20020130100011) (EC), CONICET (PIP 2015-2017 N°11220150100175CO) (MC), University of Wisconsin-Madison department of Neurological Surgery (RSA), NARSAD Young Investigator Grant from the Brain & Behavioral Research Foundation #22669 (LP), and a Ruth L. Kirschstein National Research Service Award (MH113351-02) (AM). The authors would like to thank the University of Buenos Aires for granting a fellowship for short-stay abroad (CA), which facilitated the collaboration.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the Agencia Nacional de Promoción Científica y Tecnológica [PICT-2013-0653]; Agencia Nacional de Promoción Científica y Tecnológica [PICT 2016-0218]; Brain and Behavior Research Foundation [22669]; National Institute of Mental Health [MH113351-02]; Consejo Nacional de Investigaciones Cientifi y Tecnicas (argentina) [PIP 2015-2017 N°11220150100175CO]; School of Medicine and Public Health, University of Wisconsin-Madison; Universidad de Buenos Aires [UBACYT20020130100011].