Comprehensive mapping of 5-hydroxymethylcytosine epigenetic dynamics in axon regeneration

ABSTRACT

In contrast to central nervous system neurons, dorsal root ganglia (DRG) neurons can switch to a regenerative state after peripheral axotomy. In a screen for chromatin regulators of the regenerative responses in this conditioning lesion paradigm, we identified Tet methylcytosine dioxygenase 3 (Tet3) as upregulated in DRG neurons, along with increased 5-hydroxymethylcytosine (5hmC). We generated genome-wide 5hmC maps in adult DRG, which revealed that peripheral and central axotomy (leading to no regenerative effect) triggered differential 5hmC changes that are associated with distinct signaling pathways. 5hmC was altered in a large set of regeneration-associated genes (RAGs), including well-known RAGs, such as Atf3, Bdnf, and Smad1, that regulate axon growth potential of DRG neurons, thus supporting its role for RAG regulation. Our analyses also predicted HIF-1, STAT, and IRF as potential transcription factors that may collaborate with Tet3 for 5hmC modifications. Intriguingly, central axotomy resulted in widespread 5hmC modifications that had little overlap with those of peripheral axotomy, thus potentially constituting a roadblock for regeneration. Our study revealed 5hmC dynamics as a previously unrecognized epigenetic mechanism underlying the divergent responses after axonal injury.

KEYWORDS:

• 5hmC
• axon injury
• axon regeneration
• conditioning lesion
• DNA methylation
• DRG
• epigenetics
• regeneration-associated genes (RAG)
• Tet

Introduction

DNA methylation regulates diverse sets of genes that control cellular identity and differentiation state.¹ Initially thought to be a stable and heritable modification, DNA methylation is in fact a highly dynamic process, serving as a basic epigenetic mechanism to regulate neurodevelopment² and neural plasticity in learning and memory.³ DNA methylation is also frequently altered in CNS disorders.^4,5 The identification of enzymes responsible for DNA (de)methylation further highlights its dynamic nature. The addition of methyl group to cytosine (5mC) is catalyzed by DNA methyltransferases (Dnmts). DNA demethylation, on the other hand, occurs through an intermediate base, 5-hydroxymethylcytosine (5hmC), and is catalyzed by the Tet methylcytosine dioxygenase family (Tet1–3).⁶

Besides serving as an intermediate in DNA demethylation, 5hmC may have regulatory roles in its own right.⁷ The 5hmC modification is relatively abundant in the central nervous system (CNS), —more than 10-fold higher than in embryonic stem (ES) cells—, and accounts for ∼40% of modified cytosines in CNS^8–10 Comprehensive genome-wide mapping of 5hmC has been performed in ES cells,¹¹ and in CNS tissues at different developmental stages^7,12 or under pathological conditions, such as cocaine abuse.¹³ Collectively, they form the foundation to understand the influences of this base modification in neurodevelopment and CNS disorders. In contrast, virtually nothing is known about the 5hmC dynamics in nerve injury and axon regeneration.^14,15 Hence, insights into the 5hmC state would provide a new perspective for understanding epigenetic regulation of regenerative injury responses.

Compared to their embryonic counterparts, adult CNS neurons exhibit a diminished capacity for axon regrowth.^16,17 This is partly attributed to epigenetic repression of pro-growth genes or stable expression of growth-inhibiting genes after completion of axonal wiring.¹⁸ To assess the role of 5hmC in the regenerative responses after axon injury, we took advantage of the well-established conditioning lesion paradigm in sensory neurons in dorsal root ganglia (DRG). DRG neurons share unique features of both CNS and peripheral nervous system (PNS) neurons, as they project both a central and a peripheral axon branch. Whereas central axotomy leaves DRG neurons in a state refractory to regeneration, peripheral axotomy switches DRG neurons into a regenerative state that promotes regrowth of both peripheral and central axon branches.^19,20 The conditioning lesion effect is transcription dependent, highlighting the importance of gene regulatory mechanisms for axon regeneration.²¹

Earlier effort has been focused on identifying regeneration-associated genes (RAGs),^22-25 yet little is known of the transcriptional mechanisms underlying the induction and sustained expression of the pro-growth genes. Dynamic changes of chromatin landscapes are proposed to set the stage for coordinated regulation for entire classes of RAGs.¹⁸ Previously, our laboratory has identified histone acetylation as a critical epigenetic mechanism in the conditioning lesion paradigm.²⁶ We showed that increased histone acetylation promoted access of the pro-regenerative transcription factor Smad1 to target genes. Increasing histone acetylation levels by specific HDAC inhibitors resulted in orchestrated transcriptional changes of a large repertoire of RAGs, as well as enhanced axon growth potential and sensory axon regeneration in a spinal cord injury (SCI) model.²⁶ However, HDAC inhibitors induced only a fraction of the RAGs tested, indicating engagement of additional regulatory mechanisms for the regenerative responses after the conditioning lesion.

To identify novel chromatin regulators of regenerative injury responses, we screened for epigenetic factors that are differentially regulated in adult DRG after peripheral axotomy as compared with naive state with no injury. We identified Tet3 as specifically upregulated in conditioned DRG, along with elevated 5hmC levels. To understand the influence of 5hmC reconfigurations in the regenerative responses, we generated comprehensive epigenomic maps of 5hmC distributions and dynamics under different regenerative states of adult DRG. While a large number of genomic loci displayed stable 5hmC configurations, axonal injury initiated widespread 5hmC alterations. A concept emerging from our analyses is that contrary to previous notion, central axotomy of DRG neurons is not a static event, but instead a dynamic process that results in genome-wide 5hmC reconfigurations with little overlap with peripheral axotomy. Hence, central axotomy appears to engage divergent 5hmC-regulated signaling pathways, which may constitute a roadblock for regeneration. Furthermore, we demonstrated that conditioning lesion resulted in unique acquisition or deletion of 5hmC modifications in about half of RAGs, many of which play important roles in neurite outgrowth. Our study also predicted a set of transcription factor families, including hypoxia-inducible factor 1 (HIF-1), signal transducer and activator of transcription (STAT), and interferon regulatory factor (IRF), that may assist Tet3 for 5hmC modifications. Together, our results suggest a novel role for 5hmC epigenetics in regulating regenerative injury responses.

Results

Tet3 is upregulated after conditioning lesion

In order to identify novel epigenetic regulators during axon regeneration, we screened by qRT-PCR array transcriptional changes of over 80 chromatin regulators in adult DRG 24 h after sciatic nerve transection (peripheral lesion, PL) as compared with contralateral uninjured control DRG (ctrl). We identified Tet3 as the most upregulated gene after the conditioning lesion (Fig. 1A). We then compared transcriptional changes of Tet family members (Tet1–3) after PL or T8 dorsal column transection (central lesion, CL). Tet3 was the only member of the Tet family that was upregulated at 24 h in DRG after PL relative to CL (Fig. 1B). Time-course analysis revealed that Tet3 was significantly induced in PL-DRG at 12 h post-lesion, remained elevated at 24 h, but started to decline after 3 d (Fig. 1C).

Figure 1. Induction of Tet3 in conditioned DRG (A) Screen for transcriptional changes of epigenetic regulators in conditioned DRG vs. contralateral naive DRG by qRT-PCR array at 24 h after lesion. Tet3 showed highest induction among genes tested (blue box). (B) qRT-PCR showed specific induction of Tet3, but not other Tet family members in conditioned DRG with peripheral lesion (PL) as compared with DRG with central lesion (CL) at 24 h after injury. qRT-PCR of conditioned DRG also showed upregulation of Gadd45a, but not of other enzymes involved in DNA (de)methylation. Results were normalized to housekeeping gene Rpl13a. (C) Time-course analysis of Tet3 induction by qRT-PCR comparing conditioned PL-DRG to contralateral control DRG. (D) Diagram of DNA (de)methylation pathways. (E) Immunohistochemistry showed upregulation of TET3 (green, arrows) in the nuclei of conditioned PL-DRG neurons (NeuN, red) at 24 h after peripheral lesion as compared with control with no injury. DAPI (blue) was used for nuclear counterstaining. (F-G) Immunocytochemistry of cultured DRG neurons showed levels of TET3 (F) and 5hmC (G) in nuclei of DRG neurons at 2h, 24 h and 96 h in vitro. Scale bars: 50 μm (E), 25 μm (F-G). Representative images from three independent sets of experiments are shown.

As TET enzymatic activity plays a key role in DNA hydroxymethylation/demethylation (Fig. 1D), we performed additional qRT-PCR analysis on genes involved in regulation of DNA methylation. We detected a marked induction of growth arrest and DNA damage-inducible gene 45 α (Gadd45a) in conditioned PL-DRG at 24 h as compared with CL-DRG, but no significant changes in other genes tested (Fig. 1B). Consistently, Gadd45a, a regulator of active DNA demethylation ,²⁷ has been shown previously to be a highly upregulated RAG in the conditioning lesion paradigm.²⁸ Tet3 upregulation was confirmed by immunohistochemistry (IHC) in conditioned DRG, with predominantly nuclear localization in sensory neurons (Fig. 1E).

Cultured DRG neurons undergo a process that is analogous to the in vivo conditioning lesion, with similar regenerative responses, as the cell dissociation step severs both axon branches.^29–31 We therefore conducted immunocytochemistry (ICC) on cultured DRG neurons, which allows more detailed cellular analysis. We found that Tet3 protein levels were low in freshly dissociated DRG neurons, increased in nuclei at 24 h, and waned by 96 h after plating, mirroring the in vivo results (Fig. 1F). Tet3 induction hence occurs early after axonal injury, supporting the importance of early modulation of chromatin states to facilitate RAG induction.

Upregulated levels of Tet3 have been shown to lead to accumulation of 5hmC,^32–34 we therefore examined 5hmC levels by ICC, which confirmed elevated 5hmC immunosignals in the nuclei of cultured DRG neurons (Fig. 1G). Specifically, freshly dissociated DRG neurons showed a moderate baseline level of 5hmC at 2 h after plating, but a marked increase at 24 h in culture, consistent with initiation of an in vitro conditioning program. Notably, the 5hmC increases were sustained at 96 h in vitro, even though Tet3 upregulation had waned (Fig. 1F, G). These data confirmed a correlation between Tet3 induction and 5hmC accumulations in conditioned DRG neurons and supported the model that while Tet3 expression changes may be relatively transient, 5hmC may serve as a lasting epigenetic mark.^7,35

Dynamic genomic 5hmC patterns in adult DRG in different regenerative states

To identify genome-wide 5hmC reconfigurations in adult DRG under different regenerative states, we collected lumbar 4–6 DRGs at 24 h after either PL or CL in independent triplicate experiments. Contralateral, uninjured DRGs from mice subjected to PL were used as controls. Extracted chromatin was subjected to 5hmC capture using the β-glucosyltransferase-based covalent method (Hydroxymethyl Collector),⁹ followed by genome-wide short-read sequencing. Around 60 million uniquely mapped, non-duplicated reads per sample were obtained (Table S2). As input controls, DNA samples without 5hmC enrichment were sequenced for each condition.

We first assessed global distributions of 5hmC signals using input-normalized 10 kb binned read count data. Principal component analysis revealed that the triplicates of each DRG condition clustered together (Fig. 2A), indicating both reproducibility of each replicate and highly similar patterns of genomic 5hmC configurations within each DRG condition, but divergent 5hmC patterns among the three DRG conditions.

Figure 2. Genomic mapping revealed distinct 5 hmC distributions in adult DRG under different injury conditions (A) Principal component analysis (PCA) of input-normalized 5hmC binned data (10 kb bins) showed clustering of three independent biological replicates for each DRG injury condition. (B) Hierarchically-clustered heatmap of input-normalized 5hmC signals in genome-wide binned data (10 kb bins, filtered for RPKM>1) in adult DRGs under different regenerative conditions. Blue box highlights cluster of bins with conditioning lesion-specific 5hmC gain, and orange box highlights bins with conditioning-lesion-specific 5hmC loss. (C) Global visualization of 5hmC signals in and around gene bodies by ngs.plot. Read densities at all RefSeq transcripts were summed for each DRG injury condition and distributed from 3 kb upstream of transcription start site (TSS) to 3 kb downstream of transcription end site (TES). (D) Bar graph showing numbers of all 5hmC peaks in different DRG conditions as determined by MACS (top), and numbers of gene-associated 5hmC peaks (from TSS-3kb to TES+1kb; bottom). (E) Distributions of 5hmC MACS peaks per chromosome in different DRG conditions, as compared with values expected according to chromosome size or number of encoded genes. See also Table S2.

When the binned 5hmC data were compared and visualized as a hierarchically-clustered heatmap, we could distinguish genomic loci with differential 5hmC configurations among the control, CL, and PL DRG conditions (Fig. 2B). Notably, one cluster displayed 5hmC gains in PL relative to ctrl or CL-DRG, while another cluster showed 5hmC depletions in PL-DRG. It was also evident that a substantial number of loci displayed 5hmC alterations in CL relative to ctrl or PL-DRG (Fig. 2B). This suggests that even though CL induced neither Tet3 upregulation nor enhanced regenerative potential, it nevertheless resulted in widespread 5hmC reconfigurations throughout the genome. Together, our data lend evidence for dynamic 5hmC modifications in adult DRG induced after either peripheral or central axotomy.

Earlier studies revealed a predominant association of 5hmC with gene bodies,¹¹ we therefore generated global visualization of 5hmC read densities in and near gene bodies, which showed that the distributions of 5hmC densities were not random; instead, the highest 5hmC densities were located in gene bodies and in regions ∼1 kb upstream of transcription start site (TSS) (Fig. 2C). In contrast, there was a gradual trailing along 3 kb downstream of transcription end site (TES). Moreover, 5hmC densities were markedly lower at the TSS sites, similar to that in earlier 5hmC mappings in CNS tissues.¹² Of note, the overall genic 5hmC densities patterns were comparable among the three DRG conditions (Fig. 2C).

Peak-calling analysis revealed that the genome-wide 5hmC peak number was highest in PL-DRG (314,203) as compared with CL-DRG (274,544) or ctrl-DRG (195,662) (Fig. 2D). PL-DRG also had the highest number of gene-associated 5hmC peaks (18,520) as compared with CL-DRG (18,032) or ctrl-DRG, which had far less (16,206) (Fig. 2D).

Chromosomal distribution analysis showed a relatively even distribution of 5hmC peaks with a noticeable depletion on the X chromosome (Fig. 2E). Substantial depletion of 5hmC on chromosome X has also been observed in CNS tissues,¹² supporting conserved 5hmC-associated regulatory mechanisms on chromosome X in both CNS and PNS. In comparison, the relative distributions of 5hmC peaks on autosomes were largely proportional to chromosome size and gene numbers encoded in each chromosome in all three DRG conditions with a few exceptions such as chromosome 1 and 7 (Fig. 2E). Together, our results demonstrate that genomic features of 5hmC distributions in PNS are similar to those reported in CNS.¹² Our data also highlight specific loci with 5hmC changes induced by the conditioning lesion, thus raising the tantalizing possibility of a link between 5hmC modifications and the regenerative responses/RAG regulation in conditioned DRG.

Central and peripheral axotomy lead to distinct differentially hydroxymethylated regions

To comprehensively localize regeneration-associated differentially 5-hydroxymethylated regions (DhMRs) on a genome-wide scale, we first conducted pairwise comparisons of 5hmC signals between different DRG states using the diffReps program.³⁶ We identified 12,244 DhMRs in PL-vs-ctrl comparison, and 22,944 DhMRs in PL-vs-CL comparison (nominal P-value < 10⁻⁴; Supplemental Data Sets 1 and 3). There was a relatively even split of 5hmC gain and loss for PL-vs-ctrl DhMRs, and a bias toward 5hmC loss in PL-vs-CL DhMRs (3:1 ratio, Fig. 3A). There was nearly three times the number of CL-vs-ctrl DhMRs (35,242), with a strong bias toward 5hmC gain (5:1 ratio; Supplemental Data Set 2). The fact that more DhMRs display 5hmC gain in CL than in PL may suggest that increased levels of Tet3 in PL-DRG could contribute to full cytosine demethylation. Of note, the number of DhMRs may also be influenced by technical variations among the triplicates.

Figure 3. Distinct 5hmC reconfigurations after peripheral and central lesion of DRG (A) Bar graph showing the number of differentially hydroxymethylated regions (DhMRs) among pairwise comparisons between different DRG injury conditions. The numbers of DhMRs with 5hmC gain or loss are listed on the right. (B) Pie charts showing association of DhMRs with genomic features. (C) Pie charts showing association of DhMRs with CpG islands, shores, shelves, and open sea. (D) Bar graph showing the number of genes containing DhMRs for each comparison. The proportion of genes containing single or multiple DhMRs are indicated by colored boxes. (E) Venn diagram revealing small overlap between the two sets of DhMRs induced by CL or PL relative to control (left), and between the two sets of DhMRs that are associated with gene bodies (right). (F) Alignment of genome-wide DhMRs according to the extent of 5hmC gain or loss highlighting little overlap between different sets of pairwise comparisons. (G) Hierarchical clustering of enriched pathways associated with either CL-vs-ctrl or PL-vs-ctrl DhMRs by Ingenuity Comparative Analysis shows the top differential categories for Diseases and Functions (left) and Canonical Pathways (right).

Genomic distribution analyses showed that ∼55% of DhMRs were in gene bodies, ∼5–7% at 1–3 kb promoter regions of TSS, ∼2% within 1 kb promoter regions of TSS, and ∼30% in intergenic regions (Fig. 3B). In relation to CpG islands, the overwhelming majority of DhMRs (∼90%) was in “open sea,” while ∼5% were detected at CpG shelves, ∼6% at CpG shores, and less than 0.5% on CpG islands (Fig. 3C).

We next determined the number of genes containing the DhMRs. There were nearly twice as many genes containing CL-vs-ctrl DhMRs (9101) than PL-vs-ctrl DhMRs (4736), while 7229 genes contained PL-vs-CL DhMRs (Fig. 3D). Notably, a significant number of genes contained multiple DhMRs, ranging from 2 to as many as 55 (Fig. 3D).

To determine whether PL and CL resulted in similar or distinct 5hmC modifications, we intersected the two sets of DhMRs (PL-vs-ctrl and CL-vs-ctrl), and surprisingly found only ∼3% overlapping DhMR clusters (1051). Similarly, there were only 3% overlapping gene body-associated DhMRs between the two sets (Fig. 3E). Of note, on rare occasions, neighboring DhMRs in one set may share the same overlapping DhMR from the other set, in which cases, they were consolidated into one DhMR cluster, thus rendering a slightly decreased overall number of DhMR clusters than that of DhMRs (compare Fig. 3A and 3E). To further compare genomic distribution of DhMRs, we aligned all the 70,430 DhMRs, sorted according to the extent of 5hmC gain or loss for each pair-wise comparison, and it was evident that PL and CL resulted in markedly different 5hmC modification patterns with little overlap between the two (Fig. 3F).

For functional annotation of DhMR-associated genes, we performed comparative Ingenuity Pathway Analysis (IPA) of gene sets associated with PL-vs-ctrl or CL-vs-ctrl DhMRs. Remarkably, the most divergent Diseases and Functions categories include growth of axons, outgrowth of neurites, and regeneration of neurites according to their activation z-score (Fig. 3G). Likewise, the top differential DhMR-associated Canonical Pathways included many of the signaling pathways known to regulate axonal growth and regeneration, such as PTEN³⁷ and BMP signaling,³⁸ as well as pathways important for cytoskeleton dynamics, such as RhoGDI and actin cytoskeleton signaling (Fig. 3G). Of note, the directionality of the z-score is dictated by the direction of 5hmC changes in DhMRs, and thus does not necessarily indicate increased or decreased functional activation or signaling strength. Together, our results indicate that PL and CL elicit remarkably different genomic patterns of 5hmC modifications, which exert divergent influences on signaling pathways important for neurite outgrowth. These results also provide a fresh perspective on the underlying epigenetic mechanism of the conditioning effect induced by PL in contrast to a repressed regenerative state induced by CL.

5hmC dynamics in RAG regulation after conditioning lesion

Previously, we have performed Affymetrix array studies to identify differentially expressed genes in adult DRGs after conditioning lesion compared with central axotomy.³¹ This comparison has the benefit of excluding genes triggered by injury alone. A total of 365 RAGs were identified in PL-vs-CL comparison after 24 h post-injury, with 231 genes being upregulated and 134 downregulated (Fig. 4A). To explore the relationship of 5hmC alternations and transcriptional regulation, we intersected this set of RAGs with the PL-vs-CL DhMRs. We found that among the 365 RAGs, 161 genes (44%) contained DhMRs, among which 99 were upregulated and 62 downregulated RAGs (Fig. 4A). Hence, nearly half of the PL-induced RAGs displayed 5hmC alterations, suggesting an active role of 5hmC-associated mechanisms in the transcriptional regulation of RAGs.

Figure 4. 5hmC dynamics and RAG regulation in conditioned DRG (A) Table listing the number of differentially regulated regeneration-associated genes (RAGs) in adult DRG at 24 h after PL compare with CL, and the number of RAGs that contained PL-vs-CL DhMRs. (B) Pie charts showing distribution of DhMRs associated with RAGs regarding genomic features (left) and CpG densities (right). (C) Diagram showing distribution of DhMRs associated with RAGs according to directionality of transcriptional and 5hmC changes after PL compare with CL at 24 h post-injury. DE: differentially expression genes. DR: DhMR. (D) Tables listing the top up- (left) or downregulated RAGs (right) in PL-DRG as compared with CL-DRG and the characteristics of the associated DhMRs. (E-F) Examples of genomic views of 5hmC changes of the three different DRG conditions in two highly regulated RAGs (E), and genomic views of PL-vs-CL DhMRs in the indicated RAGs (F), with green boxes indicating gain in 5hmC, and red boxes loss in 5hmC. (G) Pie chart showing the proportion of RAGs containing one or multiple DhMRs. See also Fig. S1.

The genomic distribution analysis of the 370 DhMRs associated with the 161 RAGs showed that they were mostly located in gene bodies (88%), with only a small fraction located at promoter regions (10% at 1–3 kb of TSS and 2% within 1 kb of TSS) (Fig. 4B). Even for these actively regulated genes in the conditioning lesion paradigm, there was no apparent association of DhMRs with CpG islands in that 5hmC changes occurred mostly in “open sea” (Fig. 4B). The 161 RAGs were distributed relatively evenly throughout the genome with noticeable paucity on chromosome X (Fig. S1).

We next analyzed correlation of transcriptional changes of these 161 RAGs with 5hmC alterations (Fig. 4C). Activating transcription factor 3 (Atf3) was the top upregulated RAG that also displayed largest 5hmC gain (Fig. 4D). Atf3 is a well-known RAG that enhances the intrinsic growth state of DRG neurons.³⁹ It contained one PL-vs-CL DhMR with 5hmC gain in exon 2 (Fig. 4E). Brain-derived neurotrophic factor (Bdnf), an important neurotrophin regulating axon growth,⁴⁰ is another top ranked upregulated RAG with concomitant 5hmC gain in the last exon (Fig. 4E). Sema6a, ranked as the second upregulated RAG, is a transmembrane axon guidance molecule with axon repulsion activity through Plexin-A2.^41,42 It encompasses two PL-vs-CL DhMRs located in introns, with one showing 5hmC gain and one 5hmC loss (Fig. 4F). Smad1, a regeneration-associated transcription factor that has been extensively characterized in our earlier studies, is upregulated specifically after PL, and its activation enhances axon growth potential of DRG neurons.^31,38 It contained two DhMRs in introns, both showing 5hmC loss (Fig. 4F). Among the most downregulated RAGs, many contained multiple DhMRs (Fig. 4G). For instance, Ntrk2 (also known as TrkB, a BDNF receptor) contained 5 DhMRs with 5hmC gain and 6 DhMRs with 5hmC loss, all of which are located in introns except for one located in the last exon (Fig. 4F).

We next tested the model that gene upregulation may be correlated with 5hmC gain, and, conversely, gene downregulation with 5hmC loss. However, when the 370 RAG-associated DhMRs were classified into four groups according to the directionality of 5hmC and transcriptional changes, there was no clear correlation, with RAGs scattered in all four quadrants (Fig. 4C). This may be attributed to the fact that close to half of the RAGs contained multiple DhMRs (Fig. 4D and 4G), often with divergent 5hmC changes in a single gene. We then selected one DhMR with the largest absolute fold-change to represent that gene, but still found no obvious correlation as suggested by the model: only around half of the RAGs showed concordant 5hmC changes in relation to the direction of transcriptional changes (Fig. S1B). We next focused on promoter region DhMRs, which represent only 12% of RAG-associated DhMRs. Even when separating promoter from gene body DhMRs, there was still no clear evidence that the directionality of 5hmC changes was correlated with the direction of transcriptional changes (Fig. S1C). We thus conclude that complex regulatory mechanisms must be at play to convey the 5hmC modifications to genomic output. It is noteworthy that the functional annotation of multiple loci within one gene displaying 5hmC modifications remains unclear at this point.

Unique 5hmC modifications specifically regulated in the conditioning lesion paradigm

Next, we reasoned that comparisons of 5hmC distributions across the three DRG conditions would allow us to identify DhMRs that are unique to the conditioning lesion or commonly shared by both PL and CL, the latter of which reflects injury-induced 5hmC changes. A Venn diagram of overlapping DhMRs from the three pair-wise comparisons delineated seven distinct DhMR groups (Fig. 5A and S2–3). Group a, or PL-unique DhMRs, consisted of 1,036 DhMR clusters that exhibited altered 5hmC signals specifically in PL relative to both ctrl and CL-DRG (Fig. 5A, Supplemental Data Set 4). Among them, 552 loci showed 5hmC gain and 471 showed 5hmC loss (Fig. 5B). Hence, the three way comparisons allowed us to significantly narrow down the DhMRs to the ones uniquely associated with the regenerative state triggered by PL.

Figure 5. Conditioning lesion-specific DhMRs (A) Venn diagram of three sets of DhMRs from pairwise comparisons identifying DhMR clusters that were unique to PL or CL, as well as injury-associated DhMRs that were shared by both. (B) PL-specific (group a) DhMRs separated into the group with specific 5hmC gain (left) and the group with specific 5hmC loss (right). Normalized 5hmC signals (reads per million divided by bin size) were plotted continuously upstream, within, and downstream of each PL-specific DhMRs and averaged across these regions. Total number of DhMRs is indicated at the top. (C) Pie chart showing the association of group a DhMRs with genomic features (left). The numbers of genes associated with group a DhMRs (5hmC gain or loss) are shown on the right. (D) Volcano plot of PL-specific DhMRs according to 5hmC fold changes and P-values. The location of Bdnf is highlighted in red. (E) Top ranked group a DhMRs with the largest fold change (FC) of either 5hmC gain or loss. (F) Genomic views of top ranked PL-specific DhMRs, exhibiting unique gain (P2rx3 and Defb28) or loss (Ap1g1) of 5hmC after PL compared with CL or control. (G) Ingenuity annotations of gene functions of group a DhMR-associated genes. (H) Genomic views of Stat3 and Camk2d showing locations of DhMRs. See also Figs. S2-S5.

Around half of the group a PL-specific DhMRs were located in gene bodies (Fig. 5C). Among a total of 614 genes containing group a DhMRs, 284 showed 5hmC gain and 310 genes showed 5hmC loss in PL, while 20 genes showed both 5hmC gain and loss (Fig. 5C). Volcano plot identified genes with the most dynamically changed DhMRs (Fig. 5D-E). Genome browser views of the top ranked genes, e.g., P2rx3 or Defb28, clearly demonstrated intragenic 5hmC gain at introns or 5’-UTR, respectively (Fig. 5F). Ap1g1, on the other hand, showed markedly decreased 5hmC signals in an intron (Fig. 5F). Annotation of biological function of group a DhMR-associated genes by IPA identified enzymes (100 genes), transcriptional regulators (55), kinases (43), as well as transporters (36), G protein coupled receptors (8) and microRNA (5) (Fig. 5G). Notably, Stat3, a transcription factor shown to promote axon growth in DRG neurons,⁴³ displayed PL-specific 5hmC loss; calcium/calmodulin-dependent protein kinase II delta (Camk2d), a central regulator of plasticity,⁴⁴ also showed PL-specific 5hmC loss (Fig. 5H).

Interestingly, the number of loci displaying CL-specific 5hmC alterations (group d) outnumbered group a DhMRs by 10-fold (10,882 DhMRs; Fig. S2A-C; Supplemental Data Set 7). This further supports an active chromatin remodeling process initiated by CL. Notably, an overwhelming majority of group d DhMRs (10,396) showed 5hmC gain (Fig. S2B).

Group c consisted of 881 axotomy-associated DhMR clusters that were shared by both PL and CL relative to ctrl (Fig. S2D, Supplemental data set 6). Among them, 540 loci showed 5hmC gain and 309 loci showed 5hmC loss after axotomy (Fig. S2D-G). Group b had 170 DhMR clusters that were overlapped by all three pair-wise comparisons (Fig. S3A; Supplemental Data Set 5). Among them, 69 DhMR showed 5hmC gain in both PL- and CL-relative to control DRGs, but the extent of 5hmC gain was more pronounced in CL than PL (Fig. S3A). Sixty-two DhMR showed divergent 5hmC changes, with 5hmC gain in CL and loss in PL compared with control (Fig. S3A). A complete list of the other scenarios of 5hmC changes is shown in Fig. S4. Group e-g DhMRs included a large number of loci displaying 5hmC changes present in only one pairwise comparison, without overlap with other pairwise comparisons (Figs. S3B-D, Supplemental Data Sets 8–10). As stated earlier, intersecting pair-wise comparisons of 5hmC maps allowed us to narrow down to specific DhMRs reflecting regeneration- or injury-associated 5hmC signatures. For instance, PL-vs-ctrl comparison yielded 12,244 DhMRs, the majority of which belong to group e DhMRs (10,144); and only 1036 belong to group a, PL-specific DhMRs.

We further analyzed the 161 DhMR-containing RAGs induced by PL and asked whether they contain group a, PL-specific DhMRs. We found that 18 of them exhibited PL-specific 5hmC changes (Fig. S5A), including growth factor signaling genes, such as Bdnf and Ntrk2 (TrkB), as well as cytoskeletal regulators such as Nedd4l, RhoC and Arhgap24 (Fig. S5B).

Lastly, we applied Ingenuity Pathway Analysis (IPA) for enrichment of Canonical Pathways of DhMR-associated genes. In genes associated with group a DhMRs, enriched pathways included neuropathic pain signaling in dorsal horn neurons, synaptic long-term potentiation or depression (Fig. 6A). Other top-ranked pathways included NFAT, CREB, and Neurotrophin/TRK signaling, all of which have been linked to axonal growth.^45,46 (Fig. 6A). In comparison, the top enriched pathways for group d genes included Protein Kinase A, Hippo, Integrin, and Ephrin receptor signaling (Fig. 6A). It also included PI3K/AKT signaling, axon guidance signaling, ERK/MAPK signaling and PTEN signaling. In contrast, injury-associated, group c genes were enriched for Wnt, Cdc42, TGF-β, and STAT3 pathways (Fig. 6A). Together, these data implicate distinct signaling pathways as regulated by 5hmC-associated mechanisms in different axonal injury settings.

Figure 6. Canonical pathways and predicted TFs associated with DhMRs (A) Top ranked Canonical Pathways of gene sets that are associated with different groups of DhMRs, as determined by Ingenuity Pathway Analysis (IPA). (B) Top ranked over-represented TF motifs identified by MEME-AME (with adjusted P-value <0.001 as cutoff) in group a and group d DhMRs (from all genomic regions, not restricted to only promoter regions). Shared motifs are highlighted in blue. (C) Table listing the adjusted P-values showing significant enrichment of known Tet3 DNA binding motifs (shown at the bottom) in all the three groups of DhMRs. (D) Venn diagram showing overlap of TFs whose binding motifs were enriched in group a, PL-specific DhMRs and TFs identified in a study of axon regeneration associated signaling networks.²⁵ See also Fig. S6.

Enriched transcription factor binding motifs in DhMRs

To uncover potential transcription regulators of 5hmC dynamics after axonal injury, we performed transcription factor motif (TF motif) enrichment analysis among the DhMRs in group a and group d (using DhMRs of all regions, and not limited to promoter regions). One hundred overrepresented TF motifs were identified in group a DhMRs (adjusted P-value < 0.001), including Smad1, a known pro-growth transcription factor^26,47 For group d DhMRs, 361 TF motifs were overrepresented. Comparing the top 10 most significantly enriched TF motifs, we found that seven TF motifs were unique to each group, notably, IRF and STAT motifs were highly overrepresented in group a DhMRs, while forkhead box (FOX) motifs were enriched in group d DhMRs (Fig. 6B). Three TF motifs were shared, i.e., hypoxia-inducible factor 1-α (HIF1A), aryl hydrocarbon receptor nuclear translocator (ARNT), and interferon regulatory factor 4 (IRF4) (Fig. 6B). Interestingly, HIF-1 is a heterodimeric TF consisting of HIF1A, the α subunit, and ARNT, the β subunit. The fact that both α and β subunits of HIF-1 were identified as top ranked enriched TF motifs in our unbiased approach supports its central role in regulating 5hmC dynamics after axonal injury.

The CXXC DNA-binding domain of TET3 has three recognized binding motifs.³⁴ We found that CXXC motif 2 was overrepresented in Group a DhMRs (adjusted P-value < 1.3 × 10⁻⁴), while all the three CXXC motifs were enriched in group d DhMRs (Fig. 6C). We also analyzed group f DhMRs, which excluded PL-associated DhMRs, and found that all three Tet3 CXXC motifs were also enriched (Fig. 6C). These results suggest that even though Tet3 is upregulated only in PL-DRG, the baseline Tet3 activity may contribute to 5hmC reconfigurations in all three DRG conditions.

Next, we investigated whether TFs identified through our unbiased motif analysis of DhMRs correlated with signaling network-associated transcription factors that had been described to regulate RAGs after conditioning lesion.²⁵ We found that among the 39 signaling network-associated transcription factors, 14 overlapped with the ones identified in our motif search, including HIF1A, STAT family, FOXO (Fig. 6D).

Finally, we conducted Ingenuity upstream regulator analyses of the genes associated with PL- or CL-specific DhMRs. Among the top 10 most significant predicted upstream regulators, four (FOS, ERG, TGFB1, and ESR1) were shared by PL and CL, but the others did not overlap (Fig. S6), suggesting engagement of both shared and distinct regulators by PL and CL to regulate DhMR-associated genes.

Discussion

Cytosine hydroxymethylation is an epigenetic mark that modifies DNA-protein interactions, thereby influencing transcriptional output and cellular state. To set up a platform to understand the function of 5hmC dynamics in axon regeneration, we constructed detailed epigenomic maps of 5hmC modifications in adult DRG under three different regenerative conditions: peripheral lesion, which switches DRG neurons into a growth state; central lesion, which leaves DRG neurons in a state refractory to regeneration; and naive state with no injury. Functional annotations and comparative analyses of 5hmC maps provide a new perspective of a previously unrecognized epigenetic mechanism underlying the divergent regenerative responses after axonal injury.

Common genomic features of 5hmC distributions in PNS and CNS

The enhanced regenerative potential of PNS neurons after conditioning lesion provides an entry point to understand epigenetic mechanisms underlying the regenerative injury responses and coordinated regulation of RAGs that are absent in CNS neurons. Previous mapping studies have revealed distinct genomic distributions of 5hmC in ES cells and brain tissues.^9,11,12 To our knowledge, the current study represents the first genome-wide mapping of 5hmC in the PNS tissue. We revealed common genomic features of 5hmC configurations in CNS and PNS. For instance, there is a striking depletion of 5hmC marks on chromosome X in adult DRG, similar to that in CNS tissues, such as hippocampus and cerebellum.¹² This suggests engagement of similar 5hmC regulatory mechanisms in both CNS and PNS. Another conserved feature of 5hmC distributions between CNS and PNS is a marked low density of 5hmC near TSS sites and enriched 5hmC at gene bodies.^7,12 It is noteworthy that in ES cells, 5hmC has been found to accumulate at certain TSS sites, where it is associated with repressor complexes.^11,48

Previous studies on CNS tissues indicated that 5hmC signals remain relatively stable in certain genomic loci, but exhibit localized changes in others.¹² Our results similarly showed relative stability of 5hmC signals at many loci in adult DRG under injury conditions, echoing the findings in CNS studies. However, dynamically regulated loci are also abundant, with some exhibiting 5hmC modifications specific to peripheral lesion, others to central lesion or both. The intersection of three pairwise comparisons allowed us to narrow down the list of genomic loci that show PL-specific or CL-specific 5hmC modifications, as well as injury-associated 5hmC changes, thus providing additional understanding of this base modification in different injury settings.

5hmC-mediated regulation of RAGs after conditioning lesion

Overexpression of a single RAG, such as ATF3 or STAT3, is usually not sufficient to activate the axon regeneration program.^39,49 Instead, coordinated activation of many RAGs involved in multiple pathways is likely needed. Epigenetic changes hold potential for regulating entire classes of growth supporting or growth inhibiting genes. We found that nearly half of the RAGs induced by PL display 5hmC modifications at 24 h after lesion, including well-know RAGs, such as Atf3 and Bdnf, thus supporting an active role for 5hmC-associated mechanism in their regulation. Our findings on 5hmC-associated regulation of Bdnf in conditioned DRG also extends earlier reports on an active role of DNA methylation for activity-dependent gene regulation of Bdnf during neurodevelopment⁵⁰ and learning and memory.⁵¹

The mechanism by which 5hmC modifications regulate transcription remains to be determined. Oxidation of 5mC reduces the affinity of 5mC-specific binding proteins, including HDAC-associated repressive complexes, thereby switching transcriptional repression to transcriptional activation or a permissive state.⁵² In this regard, 5hmC may function primarily to offset the repressive effect of 5mC by simply reducing 5mC levels. On the other hand, 5hmC may interact independently with specific DNA binding proteins to regulate transcription. It has been reported that methyl-CpG binding protein, MeCP2, can bind to 5hmC in gene bodies.⁵³ Of note, we find concordant changes in the direction of 5hmC modification and transcriptional alteration only in a subset of RAGs. Hence, for other RAGs, 5hmC may serve as an intermediate for cytosine demethylation. Our analyses of separating out gene body and promoter DhMRs also did not find evidence that the directionality of 5hmC changes predicts the direction of transcription changes. Another confounding factor is that around half of DhMR-associated RAGs contain multiple DhMRs, often displaying both 5hmC gains and losses. The proportion of the RAGs containing multiple DhMRs appeared similar to that of all DhMR-associated genes, which argues against a model whereby multiple DhMRs in a gene may denote more active regulation. The functional significance of multiple DhMRs in one gene requires future clarification and complex regulatory mechanisms may be at play to convey the 5hmC configurations into genomic output.

Our data did not distinguish 5hmC at CpG vs. non-CpG cytosine methylation sites (CpH; H stands for A, C, or T). DNA methylation at CpH sites has been identified for nearly 25% of all methylated cytosines in mouse ES cells, and is enriched in gene bodies and implicated in regulation of developmentally activated genes.^11,54 Higher resolution maps will be needed to resolve 5hmC modification on a single base level. A previous DNA methylation microarray analysis found only a modest number of genes differentially methylated after PL or CL, but none of the RAGs examined displayed significant alterations in DNA methylation at promoters or CpG islands.¹⁵ In fact, most gene promoters were not significantly methylated or showed any changes in either axon injury condition. Of note, this microarray analysis focused on gene promoters and CpG islands, thus could not directly assess quantitative methylation levels at other genomic loci. Nevertheless, it echoes our findings that the overwhelming majority of 5hmC modifications, hence DNA methylation reconfigurations, occur at gene bodies and at “open sea,” and not on gene promoters or CpG islands.

The 5hmC modification by itself is insufficient to predict transcriptional changes, since only a subset of genes containing PL-vs-CL DhMRs showed significant transcriptional changes in our Affymetrix array screen for RAGs. It suggests that a changing chromatin landscape merely sets the stage for coordinated gene regulation. Our current study focuses on one time point, i.e., 24 h after axotomy, and it is likely that widespread 5hmC changes at 24 h may lead to sustained transcriptional changes of more genes at later time points. Additionally, Affymetrix array studies might underestimate transcriptional changes as compared with next-generation sequencing, a task worthy of future endeavor. Future studies conducted at later time points will address the timeframe under which specific 5hmC changes are maintained after axonal injury. The refinement of genome-wide distributions of 5hmC modifications in specific cell types will also shed further light on the influences of 5hmC in regulating cell type specific injury responses in adult DRG.

Central axotomy may switch DRG neurons to a repressed growth state

In the conditioning lesion paradigm, following axotomy of the peripheral branch of DRG neurons, retrograde transport of injury signals to the cell body of the neuron results in induction of RAGs.⁵⁵ Peculiarly, injury to the central branch of DRG axons does not lead to such activation of RAGs, thereby resulting in regenerative failure after spinal cord injury.⁵⁶ Our study revealed distinct 5hmC dynamics as a potential epigenetic mechanism underlying major differences in intrinsic transcriptional networks and, correspondingly, distinct regenerative states of DRG induced by PL or CL.

Long thought to be a static event, our study nonetheless highlights that the central axotomy results in widespread 5hmC modifications, despite the absence of Tet3 upregulation. Clearly, baseline enzymatic activities of TET3 or other TETs may induce and modify 5hmC throughout the genome in CL-DRG. In fact, we detected nearly 3-fold the numbers of DhMRs triggered by central axotomy as compare with peripheral axotomy, with little overlap between the two. Thus, the difference in global levels of 5hmC is not as relevant as the distinct genomic patterns of 5hmC modifications in response to peripheral vs. central axotomy. Aside from changes of expression level, Tet3 activity and its recruitment to specific genomic loci may be also otherwise regulated, affecting target gene changes in axon regeneration. Strikingly, IPA comparison analysis showed divergent Diseases and Functions and engagement of differential Canonical Pathways of the gene sets associated with PL- and CL-induced DhMRs. An intriguing perspective raised by our study is that CL may initiate a detrimental injury response that actively sets DRG neurons to a growth-repressed state. Hence, regenerative strategy would need not only to mimic favorable changes of the injury responses incurred by the PL, but also to mitigate those unfavorable changes incurred by the CL.

The underlying mechanisms by which PL and CL lead to distinct 5hmC modifications remain unclear, although our IPA analysis suggests that different axotomy paradigms may engage distinct upstream regulators to collaborate with TET3 or other Tets for converting 5mC to 5hmC, or further oxidization toward unmethylated cytosine. Indeed, motif enrichment analyses revealed enriched TF binding motifs at the DhMRs, with some specific for PL or CL, and others shared by both. For instance, IRF or STAT family TFs may participate in PL-specific 5hmC modifications, while FOX family TFs may contribute to CL-specific 5hmC changes. Our earlier studies demonstrated a collaborative mechanism of Smad1 and HDAC1 in RAG regulation,²⁶ hence different regeneration-associated TFs may interact with different epigenetic factors for gene regulation in the regenerative responses.

One notable shared TF motif among DhMRs is HIF-1, with both α and β units ranked as the top overrepresented TF motifs in our unbiased approach, suggesting a central role of HIF-1 in regulating injury responses after axotomy. A recent study demonstrated a role of HIF-1 in regulating multiple injury-induced genes in DRG neurons and a contribution to the conditioning lesion effect.⁵⁷ Our study thus extends further in suggesting 5hmC dynamics as a potential mechanism by which HIF-1 regulates transcriptional changes after axotomy. It is worth discussing, however, that our studies found enrichment of HIF-1 motifs in both PL- and CL-specific DhMRs, thus additional co-regulators must be at play to impart the specificity of HIF-1 target genes after PL as compared with CL.

A recent computational analysis identified 39 hub TFs that provide cross talk between 400 redundant axonal signaling networks in response to axonal injury in DRGs.²⁵ Encouragingly, there was a substantial overlap between our list of TFs implicated in 5hmC modifications and the 39 hub TFs implicated in DRG injury response, including HIF1A and STAT. Furthermore, recent multi-level bioinformatics analyses have revealed a specific regeneration-associated gene network that is coordinately regulated by RAGs, signaling pathways, and TFs, many of which also showed overlap with our 5hmC regulated Canonical Pathways and upstream TFs such as EGR1, STAT, and SMAD.⁵⁸ We also revealed potential 5hmC-associated mechanisms in regulating the hub TFs, such as ATF3 and STAT3, both of which contain PL-induced DhMRs. Hence, our data constructed new connections between signaling pathways, epigenetic regulator Tet, and transcription co-regulators, such as HIF-1, STAT, and IRF, for coordinated gene regulation in the regenerative responses after axonal injury.

In summary, genome-wide 5hmC mapping points to major influences of 5hmC in axon injury responses and regenerative potential through coordinated regulation of entire classes of genes in collaboration with specific transcription factors. The 5hmC reference epigenome in adult DRG provides a framework for future exploration of 5hmC-mediated mechanisms in the regenerative responses and RAG regulation. Altogether, our data suggest that changing the transcriptional state of adult neurons through epigenetic factors is a worthy strategy for enhancing regenerative potential.

Materials and methods

Animal Procedures

All surgeries were performed on 6- to 8-week-old CD1 female mice (Charles River Laboratories) in accordance with the guidelines and protocols approved by the Institutional Animal Care and Use Committee at Icahn School of Medicine at Mount Sinai. For dorsal column lesion, a T8 laminectomy was performed and ascending sensory fibers were transected using iris micro scissors (Fine Science Tools) to a depth of 0.8 mm. For the peripheral conditioning lesion, the right sciatic nerve was transected at mid-thigh level. Ipsilateral lumbar 4, 5, and 6 (L4–L6) DRGs and contralateral naive L4–L6 DRGs with no injury were collected. Of note, sciatic nerve lesion affects largely DRG neurons in lumbar L4–6 (all subpopulations), whereas T8 dorsal column lesion mainly affects medium/large-sized DRG neurons with low-threshold mechanoreceptors (Aβ fibers).

Immunohistochemistry

Mice were perfused intracardially with 4% paraformaldehyde (PFA) in 0.1 M phosphate buffer, and L4, L5, and L6 DRGs were dissected. DRGs were post-fixed in 4% PFA for 5 min, cryoprotected in 30% sucrose, embedded in OCT, and cryosectioned at 10 μm thickness. Sections were permeabilized with 0.1% Triton X-100 (Sigma), treated with blocking buffer (0.1 M phosphate buffer, 5% normal goat serum), and incubated overnight at 4°C with antibodies rabbit anti-TET3 (EMD Millipore, 1:500), rabbit anti-5hmC (Active Motif, 1:500), or mouse anti-β-III-Tubulin (Sigma, 1:1,000). After incubation with fluorescent secondary antibodies (Invitrogen), the sections were counterstained with DAPI (1:1,000; Invitrogen) and mounted with Fluoromount (EMS). Images were captured with an Axio Imager.A2 (Zeiss) microscope equipped with an AxioCam MRc camera.

Genomic DNA isolation

L4, L5, and L6 DRGs were dissected from animals subjected to dorsal column lesion or peripheral lesion 24 h prior. Contralateral DRGs from animals subjected to the peripheral lesion were collected as uninjured controls. DRGs from three separate animals were pooled for each biological replicate. Tissue was digested overnight at 65°C in TNES with 0.5 mg/ml Proteinase K. RNA was removed by digestion with RNase A (Thermo Scientific) for 1 h at 37°C. After RNase treatment, DNA was precipitated using high salt and ethanol.

5hmC capture and library construction

DNA was sheared into fragments of an average size of 250 bp using the Biorupter Twin (Diagenode) device, using the “low” setting and 30 cycles of 30 s “on,” followed by 30 s “off.” Fragment size was confirmed by gel electrophoresis, and DNA concentrations were confirmed before enrichment. Enrichment for 5hmC was performed with the Hydroxymethyl Collecter-Seq Kit (Active Motif) according to manufacturer's instruction. Briefly, 4.3 μg of sheared genomic DNA was used for each sample. The DNA was treated with a β-glucosyltransferase enzyme to chemically modify 5hmC nucleotides with a sugar moiety, which was then conjugated to biotin. The biotinylated 5hmC containing DNA fragments were captured with streptavidin conjugated to magnetic beads, washed, and eluted as 5hmC-enriched DNA. The 5hmC enriched DNA samples, plus one un-enriched sample from each experimental condition, designated “input,” were prepared for Illumina sequencing using the NEBNext ChIP-seq kit (New England Biolabs). Briefly, 10 ng DNA of each 5hmC enriched or input sample were subjected to end repair, A-tailing, adaptor ligation, size selection, and PCR enrichment according to manufacturer's directions. DNA concentrations and library fragment sizes were assessed by Qubit and Bioanalyzer. Libraries were submitted to the New York Genome Center for Illumina 50 bp paired-end sequencing on an Illumina HiSeq 2,500 device, with 6 barcoded samples in each lane.

qRT-PCR

RNAs from freshly collected DRGs or cultured DRG neurons were isolated using RNeasy Micro kit (QIAGEN). A total of 500 ng of total RNA was used for reverse transcription reaction (SuperScript III RT, Invitrogen). qRT-PCR was performed using SYBR Green PCR master mix (QuantaBiosciences) in ABI PCR System. The melting curve of each sample was measured to ensure the specificity of the products. Data were normalized to the housekeeping gene RPL13a and analyzed using the ΔΔC_t method. Primers used in quantitative PCR are listed in Table S1. At least three independent experiments were performed for each condition.

Bioinformatics analyses

The 5hmC sequencing data was first controlled for quality using metrics generated by FastQC (v0.11.2).⁵⁹ Raw sequencing reads were aligned to the mouse mm9 genome using default settings of Bowtie (v2.2.0).⁶⁰ Only uniquely mapped reads were retained and the alignments were subsequently filtered using the SAMtools package (v0.1.19)⁶¹ to remove duplicate reads. For the analysis of global distributions of 5hmC signals, the mouse genome was divided into 10 kb bins and the number of aligned reads for each bin were counted. Bins with RPKM > 1 for at least one sample were kept and the rest discarded. Read counts were variance stabilized, then normalized by input to yield an enrichment value. Finally, bins with enrichment > 1 for at least one sample were retained for use in the analysis, resulting in ∼20,000 bins. A hierarchical-clustered heatmap was generated using the default R heatmap.2 function, which defaults to clustering Euclidean distance using the “complete” tree-building method.

Peak-calling was performed using MACS (v2.1.1)⁶² with the “--broad -q 0.01” setting. Differential analysis between samples was performed using diffReps (v1.55.4)³⁶ with settings of 200 bp window size, 20 bp moving step size, and a nominal P-value of 0.0001 to generate candidate regions. Annotation of called peaks and differential regions to their genomic features (promoters, gene bodies, intergenic, etc.) was performed using the Region Analysis program of diffReps (v1.55.4),³⁶ and read alignment profile plots were generated using ngs.plot (v2.47).⁶³

The transcription factor motif enrichment analysis was performed using the AME program of the MEME Suite (v4.11.1)⁶⁴ with “--method ranksum --scoring avg --rsmethod quick” settings.

Accession numbers

The 5hmC-enriched DNA sequencing data reported in this study has been deposited at the NCBI Gene Expression Omnibus (GEO) data archive and can be accessed with accession code GSE85972.

Statistical analysis

Prism GraphPad software was used for Student's t test analysis of qRT-PCR data.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Author contributions

E. L and L. S. performed all bioinformatics analyses, data interpretation, and graph preparation. A. K.-C. performed experiments and participated in all aspects of the study and data quantification. R. F. participated in experimental design and data interpretation. H.Z. designed the study, analyzed the data, and wrote the paper. All authors contributed to preparation of the manuscript.

Acknowledgments

We thank all members of Zou laboratory for discussions. This work was supported by NIH/NINDS (R01 NS073596), New York State Spinal Cord Injury Research program (DOH01-C30603GG-3450000, DOH01-C30832GG-3450000), and Irma T. Hirschl / Monique Weill-Caulier Scholar Award to H.Z.