Spinal cord injury (SCI) is described as an acute traumatic lesion of the nervous tissue in the spinal canal that results in temporary or permanent sensory/motor impairments, bowel and bladder dysfunction, or paralysis (Thurman et al., 1995). SCI leads to a variety of neurological problems, including dysfunctional sensor and motor responses. Transcriptional programs controlled by many regulatory elements and global epigenome remodeling constitute the molecular mechanisms of SCI, showing remarkable temporal and geographical precision (Crunkhorn, 2019; Hutson et al., 2019; Song et al., 2021; Sugeno et al., 2021; Sun et al., 2022). The regulatory elements involved in SCI have been characterized by several epigenetic features, including open chromatin, histone modifications, noncoding RNA regulation, and local reduction in DNA methylation abundance. However, the predominant molecular functional processes after SCI in both ongoing primary and secondary injury remain unclear. Treatment of SCI is limited, and better understanding of the mechanisms may help identify novel treatment strategies.
Many studies have revealed genome-wide genetic or epigenetic alterations in neurodegenerative and neurotraumatic diseases by high-throughput sequencing. Our previous transcriptome research into SCI showed that gene expression patterns change post-injury over time (Li et al., 2022). These transcriptome assay results provide useful resources for investigating functional processes after SCI. In addition to transcriptional analysis, epigenetic analysis and microRNA high-throughput sequencing have also been explored and are attracting attention in SCI research (Liu et al., 2009). However, alterations in DNA methylation after SCI have only been explored in a few studies (Shi et al., 2018; Boni et al., 2020; Davaa et al., 2021).
Cytosine DNA methylation (mC) is a critical epigenetic mark that is highly dynamically regulated during development and disease onset (Suzuki and Bird, 2008; Deaton and Bird, 2011). Typically, hypomethylation of gene promoters is related to gene transcription, whereas hypermethylation is associated with gene repression (Fernandez et al., 2012). Furthermore, de novo methylation has been identified throughout the early stages of development (Li et al., 1993; Smith et al., 2012). The pattern of DNA methylation during brain development has been thoroughly investigated over recent years, with evidence demonstrating that DNA methylation in specific genomic regions (particularly those correlated with chromatin open access) is critical for transcription initiation, resulting in alteration of the transcriptome during development. Additionally, the presence of non-CpG (CHG and CHH) methylation is preferentially maintained in the fetal brain and embryonic stem cells, which comprise a large population of multipotent cells (Ramsahoye et al., 2000; Lister et al., 2009).
DNA methylation has been shown to be associated with the functional response after SCI. For example, in peripheral sciatic nerve lesion model experiments, DNA methylation was implicated in neuronal cell death (Chestnut et al., 2011). DNA methylation regulates axon regeneration and growth cone development in neurons (Mahar and Cavalli, 2018). In neural stem cell models, DNA methylation was shown to be involved in glial cell differentiation and fate determination (Wu et al., 2012). DNA methylation contributed to inflammatory responses and activation in a multiple sclerosis model (Celarain and Tomas-Roig, 2020). Moreover, DNA methylation was found to be closely related to angiogenesis (Sabbagh et al., 2018). Hence, alterations in DNA methylation status following SCI may play important regulatory roles and are worthy of investigation.
To explore the role of DNA methylation during SCI, we built a reduced-representation bisulfite sequencing (RRBS) library from samples of SCI model mice at various time points after injury. Our study had several aims. First, we investigated the dynamic alterations of DNA methylation during SCI. Second, we obtained key differential methylation regions (DMRs) that might be potential targets for future studies. Finally, we explored the epigenetic regulation of gene transcription driven by DNA methylation during SCI.
Twenty-four specific pathogen-free-grade adult C57BL/6 female mice (6–8-week-old, 20–30 g) were purchased from Shanghai Jiesijie Laboratory Animal Co., Ltd. (Shanghai, China; license No. SCXK (Hu) 2018-0004). Mice were fed a specific pathogen-free standard diet and housed following Division of Laboratory Animal Medicine of Shanghai Tongji Hospital guidelines (room temperature, 25°C; relative humidity, 60%). The mice were randomly divided into eight groups, including a sham group and seven SCI groups (days 0, 1, 3, 7, 14, 28 and 42 post-SCI) (n = 3/group). Surgical procedures and postoperative care were performed following protocols approved by the Institutional Animal Care and Use Committee of Shanghai Tongji Hospital and Tongji University (approval No. 2020-DW-005) on May 12, 2020. All surgeries were performed under aseptic conditions.
Crush SCI at T9
The surgical procedure for establishing a T9 spinal cord crush injury model was described in detail in our previous study (Li et al., 2020). Mice were anesthetized with inhaled 1.5% isoflurane (RWD Life Science, Shenzhen, China) mixed with carbon dioxide (95% O2/5% CO2) and were maintained under anesthesia during surgery. Under aseptic conditions and microscopy, T8–T10 spinous processes were exposed with microscissors. A sterile laminectomy was performed at the T9 vertebral level without disrupting the dura. We developed a forceps crush model with SCI in mice using No. 5 Dumont-type forceps (Fine Science Tools, North Vancouver, Canada) with force applied for 5 seconds. The sham group received identical treatment, including exposure, laminectomy, and forceps placement around the spinal cord, but crush injury was not induced. After the injury or treatment, the muscles and skin were sutured. This procedure resulted in hind limb paralysis in all experimental groups. Immediately following the surgery, the mice received Ringer’s solution for hydration (Procell, Wuhan, China) (1 mL, subcutaneous injection) and ibuprofen (Easton Biopharmaceuticals, Chengdu, China) (0.01 mg/kg) daily for 2 days. The bladder was expressed on SCI-exposed mice twice a day until spontaneous micturition. Sham mice underwent laminectomy and they were treated equally with SCI-exposed mice after surgery; however, they did not need bladder compression.
At the indicated time points after SCI (days 0, 1, 3, 7, 14, 28 and 42), mice were anesthetized with inhaled 1.5% isoflurane mixed with carbon dioxide and executed by cervical dislocation. Spinal cord tissue segments (1 cm in length and encompassing the lesion site) were dissected and used to produce an RRBS library (Li et al., 2022). The samples were prepared with Ovation RRBS Methyl-Seq System 1–16 (Tecan, Männedorf, Switzerland) following the manufacturer’s recommended protocol. Briefly, 100 ng genomic DNA from spinal cord tissue was used as input DNA. Bisulfite conversion was performed using the EpiTect Fast DNA Bisulfite Kit (Qiagen, Germantown, MD, USA). Post-library quality control was performed with the Qubit dsDNA High Sensitivity fluorometry assay (Invitrogen, Carlsbad, CA, USA). An equimolar pool of the prepared libraries was created at a concentration of 5 nM. Libraries were prepared by 100-bp paired-end sequencing using an Illumina HiSeq 2000 platform (Novogene Co. Ltd., Beijing, China).
Raw data processing
DNA reads from the RRBS library were aligned to the bowtie2-indexed reference genome GRCm38/mm10 by the Bismark tool (Krueger and Andrews, 2011) (Babraham Institute, Cambridge, UK) using a criterion of a maximum of two mismatches in a single direction. Only uniquely aligned reads were retained. Next, polymerase chain reaction duplicates were removed for inclusion in the RRBS library using default parameters. Methylated gene calling was performed using the Methylation Extractor module. A base-pair-level differential methylation analysis was performed using the methylKit R package (Akalin et al., 2012). Bases with very low (< 3×) and excessively high read coverage (> 10×) were discarded to prevent polymerase chain reaction bias and increase the power of the statistical analyses. Before calling the differentially methylated genes, each comparison was reorganized with methylKit, and the data were combined and subjected to differential methylation analysis. Pairwise comparisons were performed by Fisher’s exact test. With a minimum of five reads in each group, a differential methylation value of 20 (on a percentage scale) and P values less than 0.05 were considered threshold values for DMRs identification.
Data exploration and plot generation
Both the methylKit (https://github.com/al2na/methylKit) and RnBeads R packages (www.rnbeads.org) offer consistent pipelines for performing similarity, hierarchical clustering, and principal component analysis analyses. Similar studies can also be performed with SeqMonk (Babraham Institute), which executes these three analyses using custom probe/filtering, such as that for 10-kb tiling and promoters. For the association study between cytosine density and methylation levels, the RnBeads package was used to investigate CpG and non-CpG contexts independently. Using SeqMonk, bean plots and continuous genetic area plots were created to show the distribution of methylation levels over time, regions, and contexts. To identify genes with comparable DNA methylation patterns, we created temporal and spatial heatmaps as previously described (Lister et al., 2013). The raw methylation data matrix was obtained with SeqMonk and processed with the ComplexHeatmap R package (https://bioconductor.org/packages/ComplexHeatmap/) after log10 transformation. Every upstream and downstream 100-kb area of a gene body was separated into 100 1-kb bins. The probe generator function in SeqMonk was used to divide the gene body region into 100 equal bins, irrespective of the length of the gene body.
DMR, Gene Ontolog, and Kyoto Encyclopedia of Genes and Genomes analyses
To identify DMRs, we performed both count-based and methylation level-based calculations. Furthermore, because our experimental design included three biological replications at each time point, DMR analysis data were processed via both logistic regression (SeqMonk) and Fisher’s exact test (methylKit). A slight difference was observed between the results obtained from these two statistical methods. From the DMR results, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis data were parsed with the clusterProfiler R package. To amplify the dynamic functional mechanisms during SCI, the Mfuzz R package was used to describe temporal changes in methylation. Locus overlap analysis was performed with the RnBeads R package. The listing of genes with DMRs (CpG, CHG, and CHH regions) within promoters was updated and is presented in Additional Table 1.
All data were analyzed using GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA, www.graphpad.com) and are expressed as mean ± standard deviation (SD) between independent controls and experimental subjects. Statistical evaluations were performed with Student’s t-test, one-way analysis of variance. Post hoc analysis of variance comparisons was made using the Bonferroni (one-way analysis of variance) correction. A P value < 0.05 was considered as statistically significant. To obtain the differential methylation regions/genes, Kruskal-Wallis test (P < 0.05 after Benjamini-Hochberg correction) was applied among the multiple groups. Pearson correlation analysis was used to perform the correlation of DNA methylation values of gene promoters and mRNA transcription levels obtained from RNA-Seq.
Reduction in non-CG DNA methylation during SCI
We generated SCI model mice as described in Methods; we used female mice because of the different DNA methylation in sex chromosomes, which may generate unavoidable batch effects. The SCI mice showed an immediate and complete locomotor loss in the lower limbs, and spontaneous recovery was observed from day 3 post-injury.
Spinal cord tissue was collected from three biological repeats at various time points after SCI (days 0, 1, 3, 7, 14, 28, and 42) for RRBS library establishment. Approximately 34–47 million reads with high bisulfite conversion rates (> 99.4%) were obtained in all 24 samples (Additional Table 2), and the mean mapping efficiency was 68.6% (mm10) (Figure 1A). The global cytosine fraction (including the CHG and CHH fractions, called the CH fraction) was predominated by non-CpG methylation marks (> 86%). However, because of the low non-CpG methylation level (0.6–0.8%), non-CpG methylation marks accounted for only 1% of the 7% genome-wide DNA methylation marks. CpG dinucleotides, which carried approximately 85% of the mC of the 14% total genome-wide cytosine, presented an mC level approximately sixfold higher than other gene regions (Additional Figure 1). To minimize the bias in the RRBS process, cytosine bases with greater than 10× coverage and in the 99.9 percentile were discarded. The remaining bases with CpG and non-CpG DNA methylation marks were included for subsequent study. Determination of the percentage of methylation marks per base in each sample revealed that the majority of loci were either hypo- or hypermethylated, with bimodal distribution predicted in all spinal cord tissues (Additional Figure 2).
The global level of DNA methylation decreased modestly after SCI (P = 0.0287, Figure 1B). This decline was largely caused by a decrease in CH methylation (P = 0.0194) and not CpG methylation (P = 0.5473, Figure 1B). To obtain a higher temporal resolution of DNA methylation throughout the entire SCI and recovery processes, DNA methylation was examined from 15 minutes to 42 days after injury. Total DNA methylation levels steadily declined until day 14 after SCI and then recovered to baseline levels. The DNA methylation levels on days 3, 7, and 14 were reduced compared with levels in the sham group (Figure 1C). CpG methylation level did not show significant variation across time points (Figure 1C). CH methylation levels on days 1, 7, and 14 were decreased compared with levels in the sham group (Figure 1C).
Functional stages of SCI as determined by DNA methylation levels
The post-SCI period was classified into three stages: the early stage (days 0, 1, and 3), intermediate stage (days 7 and 14), and late stage (days 28 and 42). RNA-Seq and single-cell RNA-Seq have identified expression patterns that distinguish these stages (Li et al., 2022).
CpG methylation showed a clear stage-specific clustering because of its dominant global methylation status during SCI, and CpG methylation (> 0.95) were similar among the stages post-SCI (Figure 2A). In contrast, non-CpG methylation mark abundance was profoundly decreased after injury, with distinct similarities across stages (Additional Figure 3A); however, these marks, particularly CHG methylation marks, were not markedly reduced in a temporal or functional pattern. Notably, analyzing the cluster differences between CpG and non-CpG methylation, the fraction of mC was considered. The SCI periods were grouped into three stages using hierarchical clustering methods. The first cluster contained the day 0, 1, and 3 groups, the second cluster included the day 7 and 14 groups, and the third cluster included the day 28 and 42 groups. These three groupings thus represented early, intermediate, and late stages (Figure 2B). Moreover, the distances between groups in a stage cluster in a principal component analysis revealed a clear temporal orientation. In the early stage, the day 0 group was closest to the sham group, and the day 3 group was the farthest from the sham group. In the intermediate cluster, the day 7 and 14 groups were not separated. The day 42 group showed more similarities with the sham group than with the day 28 group (Figure 2C). Our CpG methylation analysis revealed a significant difference between the day 3 and day 7 groups, indicating that the molecular features on the third day after SCI were most likely to be associated with the early stage of the disease. The days 7 and 14 groups (in the intermediate stage) were interspersed between the early and late stages. In accordance with the restoration of total mC levels on days 28 and 42 (Figure 1C), the distance between the SCI and sham groups in the late stage was short, suggesting that alterations to DNA methylation had been attenuated after 28 days and before initiation of the late stage. Both principal component analysis and hierarchical cluster analysis performed effectively in separate time points from day 0 to day 42. However, some biological variability among replicates and some degree of overlap among periods were also observed, indicating the changes in DNA were mild and variable after SCI (Figure 2A).
We also investigated similarities among non-CpG methylation conditions (Additional Figure 3). In contrast to the CpG methylation pattern, the CHG methylation pattern did not show a clearly characterized temporal-injury-stage-related clustering pattern (Additional Figure 3B and C). The day 28 samples showed a CHG methylation pattern that differed from that of the other samples. The clustering of CHH methylation marks was similar to that of CpG methylation marks, in which early and intermediate stage groups diverged in opposite directions (Additional Figure 3B and C).
Next, we reevaluated the DNA methylation level at the three injury stages. There was a significant decrease in the total mC levels at the intermediate stage compared with the levels in the Sham group (P = 0.0045, Figure 2D). Non-CpG methylation accounted for a significant percentage of this methylation decline (P = 0.0008, Figure 2D).
Global distribution of DNA methylation marks on genetic element variants
Gene expression is mediated by DNA methylation in promoter regions and gene bodies (Lister et al., 2013). We therefore used an unbiased approach to classify patterns of mCG and mCH within each annotated gene and in flanking regions extending 100 kb up- and downstream of genes. After normalizing the methylation pattern around each gene based on the local baseline mCG or mCH level in each sample, we combined the features into a large data matrix comprising data from 7200 individual DNA methylation measurements for each gene (three contexts, eight time points, 300 bins of each gene). Using the k-means method, 26,816 genes were identified and allocated to five clusters based on DNA methylation distribution patterns. Although the absolute methylation values were markedly different, the distribution patterns in CpG, CHG, and CHH methylation were strikingly comparable. During the early period, the methylation baseline of CpG and CHH regions was elevated, and hypomethylation at transcription start sites (TSSs) was frequently observed in these gene regions. We then examined the changes in DNA methylation in gene clusters. Neither CpG nor non-CpG methylation alterations were significant, indicating that SCI exerted a limited influence on the overall methylation pattern, at least compared with its effect on the developmental stage of the CNS. Moreover, genes in Cluster 1 showed a low DNA methylation level of the whole gene body. Clusters 2, 3, and 4 exhibited the typical DNA methylation distribution in coding genes. Cluster 5 did not show demethylation around TSSs (Figure 3A).
In cells in the spinal cord, DNA methylation marks are widely spread throughout the genome, showing universal traits as well as unique patterns. From the annotations in the USCS database (mm10), probes in at least 10 observations at five depths were found to be clustered based on genomic region. CpG dinucleotides were found in high concentrations surrounding TSSs in areas of protein-coding genes (Additional Figure 4A). A fivefold increase in the peak density in the CpG context was observed at TSSs compared with that at the baseline in the same context. In promoter regions (–2000 bp upstream of TSS), CpG dinucleotides assembled close to TSSs (Additional Figure 4A). Furthermore, the density of CpG dinucleotides surrounding CpG islands was approximately four to fivefold greater than that in shore regions (Additional Figure 4A). In accordance with its biological function in transcription, DNA methylation marks around TSSs were largely absent (Additional Figure 4B), and a negative correlation was observed between methylation level and CpG density following SCI (Additional Figure 4C). The overall distribution of the non-CpG methylation marks on promoters, 5′ untranslated regions (UTRs), TSSs, CpG islands, and other regions was comparable to the distribution of CpG methylation marks and was relatively flat because of the low density in the non-CpG context. Although the baseline non-CpG density (0.6) around gene areas was similar to the CpG context density (0.5), the peak site densities (0.75) of the former were much lower than those in the latter.
To better understand the dynamic feature of DNA methylation on gene elements following SCI, temporal CpG and non-CpG methylation data on promoters, 5′UTRs, exons, introns, and 3′UTRs were integrated (Figure 3B). CHG and CHH were shown to be particularly altered by SCI in terms of their DNA methylation status after intergenic region methylation data were excluded. Overall, both CpG and non-CpG methylation marks revealed low methylation levels in TSS regions, and the absolute methylation value was consistent with our above findings (Figure 1A and Additional Figure 1A). The CpG methylation levels in the five areas did not change considerably after SCI; in contrast, the CHG and CHH regions were particularly hypomethylated after SCI, with the most significant hypomethylation evident in the intermediate stage, and the DNA methylation levels in the early and late stages were comparable. Specifically, CHG and CHH exhibited considerably greater levels of DNA methylation in intergenic areas than in other regions. Next, we looked at the distribution of methylation marks in each location based on bean plots (Additional Figure 4D). A typical spindle shape was seen in CpG methylation when the entire genome was subjected to a consecutive tiling analysis. Intron methylation contributed the most to the high levels of methylation (methylation value = 60–80), comprising the majority of the methylation marks. In promoter areas, most CpG methylation levels were low (methylation value = 10), with only a negligible percentage of the CpG methylation levels in the promoter region found to be high (methylation value = 75–100). A bipolar distribution of high and low CpG methylation in exons was found; however, the number of exons was substantially smaller than the number of introns. An examination of CHG and CHH regions revealed negligible hypermethylation with the majority of methylation found to be less than 5. CHG and CHH are prevalent in gene areas; for example, the number of CHG and CHH regions in introns was fourfold and twofold higher than that of CpG regions, respectively. Notably, the minor hypomethylation change significantly after SCI.
DMRs after SCI
We used SeqMonk and its annotation databases to investigate DMRs in various genomic elements. We performed pairwise comparisons of promoters, exons, introns, 5′UTRs, and 3′UTRs, as well as in noncoding regions, such as noncoding RNAs, short interspersed nuclear elements (SINEs), long interspersed nuclear elements, repetitive DNA, and long terminal repeats (LTRs) between the SCI groups and the sham group. Differentially hypermethylated and hypomethylated regions in the CpG, CHG, and CHH contexts were identified by Kruskal-Wallis test (P < 0.05 after correction; the detailed statistical information of promoter methylation is provided in Additional Table 3). Global and temporal hypomethylation was detected in most of the examined regions, particularly in the non-CpG contexts (Figure 4A and Additional Figure 5A). Because of the different contexts, CH methylation marks are more than CpG methylation marks in mapping to a larger number of genes. Notably, a substantial proportion of DMRs was found in introns. More DMRs were found in annotated LTRs, SINE, and long interspersed nuclear element regions compared with other regions. Non-CpG methylation was associated with more DMRs in LTRs and SINEs than CpG methylation, indicating that DNA methylation was significantly altered in these areas following SCI. The biological significance of these findings requires additional investigation. Notably, we did not detect DMRs in annotated areas such as low-complexity regions, DNA repeats, simple repeats, signal recognition particle RNA, satellites, or transfer RNA.
Non-CpG but not CpG DMRs were significantly altered after SCI, in agreement with the previous findings, and the absolute methylation levels were much lower in non-CpG regions than in CpG regions. Both CpG and non-CpG DMRs showed hypomethylation status in the early and intermediate stages after SCI until day 14. After day 14, the methylation levels were gradually restored, although there was still a considerable degree of hypomethylation compared with that in the sham group (Figure 4A).
To elucidate the functional process associated with DNA methylation at various times post-SCI, we conducted GO analysis on promoter DMR clusters at the seven time points (Figure 4B). The results indicated that CpG methylation was primarily involved in neurogenesis-related regulatory mechanisms in both the early and late stages; specifically, synaptic reorganization, neuronal regeneration, axonal regeneration, spinal cord development, and neuronal projection formation terms were enriched. Notably, in the intermediate stage, CpG methylation was related not only to functional processes of the nervous system and nerve injuries, such as synapse and forebrain development and apoptosis but also to a variety of previously undiscovered mechanisms, such as the development of reproductive and gland systems and skeletal development. In the cases of CHG (Additional Figure 5B) and CHH (Additional Figure 6) promoters, the DMRs did not show the same characteristic functional mechanism cluster as observed for CpG promoters, but there were still some pathological mechanisms worth noting, such as transforming growth factor, insulin-like growth factor, and WNT signaling pathways, which showed more significant clustering than CpG promoter regions. The continuous methylation of CpG, CHG, and CHH on day 28 was associated with functional processes related to neuron, synapse, and axon development. Surprisingly, DNA methylation was not associated with functional processes such as glial cells, inflammatory cells, and angiogenesis; only CHG demonstrated B-cell differentiation–related clustering, which was identified on day 28.
To gain detailed information on differentially methylated genes, we investigated the distribution of promoter DMRs using volcano plots with comparison data (false discovery rate < 0.01, change difference > 10% applied to CpG regions; change difference > 3 applied to CHG and CHH regions, Additional Figure 7). The regulator factor X5 (Rfx5), zinc finger protein 451 (Zfp451), and makorin ring finger protein 3 (Mkrn3) genes were hypomethylated in the early stage. Moreover, protease-related genes, such as proteasome 26S subunit ATPase 4 (Psmc4) and glycine cleavage system protein H (Gcsh), were linked to the early stage. Genes involved in various signaling pathways, including a protein kinase (protein kinase C and casein kinase substrate in neurons 1 (Pacsin1)), voltage-gated channel (sodium voltage-gated channel alpha subunit 5 (Scn5a)), transcription factor (CCAAT enhancer-binding protein epsilon (Cebpe)), extracellular matrix retinoschisin 1 (Rs1) and extracellular matrix protein 1 (Ecm1), were identified in the intermediate stage. In the late stage, the number of genes with promoter DMRs was reduced. Furthermore, numerous genes, including Gcsh, immunoglobulin heavy chain variable 6-6 (Ighv6-6), and Rs1, maintained a differentially methylated state throughout the whole injury response. In contrast to DMRs in the CpG context, DMRs were substantially less abundant in the CHG and CHH contexts. Only a few genes were found in both the early and late stages in CHG and CHH regions. Furthermore, the non-CpG DMRs were significantly hypomethylated during the intermediate stage, indicating that genes were substantially regulated during this period.
Locus overlap analysis in RnBeads software was performed on the DMRs to enrich and annotate the regions that showed variations in genomic elements. Among the hypermethylated genes, genes encoding chromatin-binding proteins (SUZ12 polycomb repressive complex 2 subunit (Suz12); PHD finger protein 19 (Phf19); metal response element binding transcription factor 2 (Mtf2)), histone modification proteins (ring finger protein 2 (Rnf2) and Jumonji and AT-rich interaction domain-containing 2 (Jarid2)), and neurogenesis-related proteins (Notch1) were found in the early stage. Polycomb or polycomb-like/related genes constitute most of the list (seven of eight genes). These results were found throughout the stages, except the late stage, demonstrating that the activity of polycomb group (PcG) proteins was regulated following SCI. The promoter of the stemness-specific transcription factor Nanog was also hypermethylated. In contrast, the histone acetyltransferase p300 gene was hypomethylated after SCI (Additional Figure 8A and B). Notably, all hypomethylated DMRs identified in the intermediate stage of SCI were transcription factors, including CBFA2/RUNX1 partner transcriptional corepressor 2 (Cbfa2t2), LYL1 basic helix-loop-helix family member (Lyl1), RUNX family transcription factor 1 (Runx1), LIM domain-binding 1 (Ldb1), signal transducer and activator of transcription 6 (Stat6), signal transducer and activator of transcription 5b (Stat5b), E1A-binding protein p300 (Ep300), and CCAAT enhancer-binding protein alpha (Cebpa), indicating that DNA methylation of transcription factors after SCI was mediated via the epigenetic modification of PcG proteins and transcription factors (Additional Figure 8C).
Temporal characteristics of DNA methylation after SCI
Next, we identified the temporal changes in DNA methylation following SCI. The corresponding genes were classified into 16 categories from the methylation patterns of the DMRs in the promoter region (Mfuzz R package), and a functional annotation assay (clusterProfiler R package) was performed (Figure 5).
Among the CpG hypermethylated genes, persistent hypermethylation of Cluster 8, corresponding to the appendage development signaling pathway and pluripotent stem cell terms, was observed, suggesting that the expression of stem cell-related genes was suppressed following SCI. The clusters with temporary hypermethylation (Clusters 3, 4, 6, and 12) were very weakly connected with biological function of the nervous system, except for Cluster 3, which was associated with axon guidance only in the intermediate stage. Some pathways, such as cAMP, RAS, and Alzheimer’s disease signaling pathways, warrant additional investigation. Cluster 1 was hypermethylated in the early and intermediate stages but recovered to near-normal levels in the late stage; it was connected with WNT signaling and protein processing in the endoplasmic reticulum, suggesting that protein maturation might be impeded in SCI.
Among the persistent hypomethylated genes, those in Clusters 2, 5, and 13 were associated with axon guidance and axon regeneration, implicating CpG methylation may play regulatory roles in neuron regeneration after SCI. Moreover, genes in signaling pathways such as epithelial tube morphogenesis, the Hippo pathway, and PI3K-AKT signaling were hypomethylated. In a transient hypomethylated cluster, Clusters 9 and 10 showed regulation of genes related to postsynaptic specialization as well as to Notch and calcium signaling pathways. Additional findings included hypomethylation of genes in histones and endocytosis during the early stage, hypomethylation of genes in retrograde endocannabinoid signaling and epithelial cell proliferation during the intermediate stage, and hypomethylation of genes in histone modification and endocytosis during the late stage. The methylation levels in Clusters 11 and 16 fluctuated after SCI, but no pattern in neural-related processes was detected. The results of a temporal clustering analysis with non-CpG methylation were similar to those of gene annotation. Methylated CpG regions were seldom associated with the nervous system but were compatible with alterations in CpG methylation (Additional Figures 9 and 10).
DNA methylation is closely related to gene expression after SCI
DNA methylation of genes exerts a profound influence on the development, maturation, and disease incidence of the mammalian central nervous system because it is directly linked to transcription (Moore et al., 2013). To better understand the mechanism of DNA methylation and gene expression regulation following SCI, we conducted a Pearson correlation analysis on the levels of DNA methylation and gene expression in contexts and regions of the spinal cord after injury. Gene expression data were collected from a previous study (Li et al., 2022). With regard to CpG methylation, promoter and intron methylation levels were negatively connected with gene expression levels, but exon CpG methylation levels were not correlated with gene transcription levels. The methylation levels of CHG and CHH in promoter regions, exons, and introns were inversely linked with gene transcription levels (Figure 6A). Additionally, our results demonstrated a negative connection between DNA methylation in a promoter region and gene transcription at various time intervals, as shown in a heatmap in Figure 6B. CpG and CH Clusters 1 and 2 were expressed at modest levels, and their methylation levels were greater than those of the other clusters. Clusters 5 and 6 demonstrated the opposite pattern, with greater gene expression levels but lower DNA methylation levels. While the accumulated alterations in gene transcription within 15 minutes of spinal cord damage (day 0) were not significant, substantial variations in DNA methylation, particularly CHG and CHH methylation, were seen compared with those in the sham group. By examining the early and intermediate stages, we discovered that DNA methylation levels were markedly reduced in transcriptionally upregulated genes; however, DNA methylation levels did not increase as predicted in transcriptionally repressed genes.
SCI is a serious medical issue and adequate treatment strategies are lacking; the underlying mechanisms have not been fully elucidated. Preliminary transcriptome and single-cell profiling investigations have revealed that the functional response to SCI can be generally categorized into early, intermediate, and late stages on the basis of changes in molecular events (Li et al., 2022). Many epigenetic mechanisms in SCI have been extensively described in the literature (York et al., 2013; Hutson et al., 2019; Danilov et al., 2020; Zhang et al., 2020). DNA methylation is a key regulatory mechanism of gene transcription, and its role in SCI has not been fully determined.
In this study, we used an RRBS library of mouse T9 SCI to obtain several notable findings. First, the level of global DNA methylation was reduced after SCI, with most of the decline in non-CpG methylated regions. Using hierarchical clustering and principal component analysis, we identified temporal windows of methylation-related diseases, which revealed that days 0–3 represented the early stage, days 7–14 represented the intermediate stage, and days 28–42 represented the late stage. Methylation loss was the most significant in the intermediate stage; hypomethylation levels were similar in the early and late stages. Analyses of different genomic areas showed that DNA methylation was concentrated mostly in intergenic regions. CpG methylation at TSSs was lost after SCI, as expected. The methylation levels of CHG and CHH in introns were substantially higher than those in other areas. The baseline CpG and CHH methylation levels increased after SCI, whereas CHG methylation levels did not. CpG methylation levels in multiple gene regions remained virtually unaltered in the SCI groups compared with the sham group. In contrast, the levels of CHG and CHH methylation in promoter regions, 5′UTRs, exons, introns, and 3′UTRs were markedly reduced. We focused on alterations in DNA methylation in the promoter region of genes. Through statistical analysis, we identified genes that were hypermethylated or hypomethylated post-SCI and conducted clustering and functional annotation. Generally, hypomethylated genes in promoters were associated with neural regeneration, whereas hypermethylated genes were associated with functional processes such as stemness maintenance, oxidative stress, and histone modification. Unexpectedly, DNA methylation was not related to many classic molecular events, such as glial cell activation, the inflammatory response, and angiogenesis.
The function of DNA methylation in traumatic brain injury has been previously reported (Zhang et al., 2007; Duan et al., 2021). Hypomethylation of a subset of reactive microglia/macrophages is thought to be the predominant contributor to the global methylation loss in traumatic brain injury. Identifying the cell type critical for SCI-induced hypomethylation was a challenge in our investigation; furthermore, GO analysis did not provide solid evidence of a link between DNA methylation and glia-/inflammatory-related events throughout the course of SCI. Differential cytosine methylation between neurons and glia has been demonstrated in a variety of states, including mammalian CNS development (Yao and Jin, 2014), neurodegenerative disease (Gasparoni et al., 2018), and traumatic injury (Haghighi et al., 2015). Future studies are required to determine DNA methylation alterations in distinct cell types following SCI.
SCI causes a series of severe neurotraumatic pathologies, such as neuronal death and subsequent apoptosis and necrosis, glial scarring, neuroinflammation, and microenvironmental damage. While CpG is more highly conserved than CH in the mammalian postnatal central nervous system, whether hypomethylation in the CH context exerts a greater impact on post-injury spinal cord transcription should be examined in future studies. Although non-CpG methylation is critical in the adult CNS (Lister et al., 2013), the level of non-CpG methylation is reduced during development and is negligible in numerous somatic cells (Ziller et al., 2011). Furthermore, CH methylation marks accumulate in neurons but not in glia and eventually become the dominant sources of methylation-related change in the neuronal genome (Lister et al., 2013). In this study, we hypothesized that DNA methylation, particularly CH methylation, was critical for neuronal activity following SCI. However, neither CHG nor CHG were identified in genes related to neuron-related activity. More investigation is required to determine the alteration of DNA methylation in distinct cell types, or at least the differences, to distinguish neuronal cells from other cell types.
Methylation changes are particularly likely in areas with high methylation levels, such as promoters (Shearstone et al., 2011). We examined the methylation status of imprinted regions in the injured spinal cord. Methylation at satellites, DNA repeats, and low-complexity repeats, as well as noncoding RNAs, such as lncRNAs and miRNAs, was conserved. However, long interspersed nuclear element, SINE, and LTR regions were significantly hypomethylated in the intermediate stage. Future studies should investigate the functional connection between methylation of these regions with recovery after SCI.
Many profound questions could not be answered by RRBS technology of mass tissue. RRBS, which is based on a methylation-sensitive restriction enzyme (Msp1), enriched the cytosine methylation of the CCGG content around in the genome. Despite these outcomes, low genome coverage was evident because the number of CpG-containing recognition sites was too small (Hsu et al., 2017). As an alternative to the present approach, whole genome bisulfite sequencing (GBS) offers advantages for researching areas of interest outside CpG islands (Yong et al., 2016). RRBS also has another limitation: it cannot be used to discriminate between 5mC and 5hmC, which are mediated by TET protein family members (Hahn et al., 2015). Thus, it will be interesting to apply additional oxBS-Seq should be applied to reveal the precise pattern of the 5hmC modification (Booth et al., 2012). Additionally, ChIP-Seq of Tet1 should be carried out. Moreover, we would like to investigate the detailed information about how DNA methylation regulates gene transcription in different cell types. In other words, single-cell RRBS or RRBS with flow cytometry is fascinating to be performed.
DNA methylation has been implicated in numerous mechanisms of SCI, although the results conflict with certain expectations. The most important results indicated that DNA methylation is linked to neuron/axon regeneration, synapse formation, ion channels, vesicle transport, microtubule formation, and neural circuit projection construction from the early to late stages of injury. Additionally, many neurodevelopmental signaling pathways, including the WNT, transforming growth factor, and insulin-like growth factor pathways, depend on DNA methylation. DNA methylation shows clear temporal characteristics, which are strikingly comparable to those that have been identified in the development of the central nervous system. DNA methylation is also critical for downstream processes mediated by protein modification systems, such as phosphorylation, ubiquitination, glycosylation and other modifications, and the Ras, MAPK, PI3K-AKT, and Hippo pathways, indicating that it plays a vital role in protein maturation, activation, inactivation, and degradation. DNA methylation marks also interact with other epigenetic alteration marks, primarily histone modifications, indicating that the control of gene transcription by DNA methylation after SCI involves coordination with other epigenetic processes. Finally, DNA methylation is involved in regulatory transcriptional complex-related processes, particularly those involving PcG family transcription factors, which are critical for gene transcription. Taken together, our findings provide a detailed characterization of the functional role played by DNA methylation and the dynamic characteristics of this modification during SCI. These results indicate that DNA methylation plays an essential function as a transcription regulator by targeting regeneration events after SCI.
Author contributions:Study design and supervision, and SCI model establishment: LC; RRBS library establishment and data analysis: ZW; other experiments implementation: YC, RZ; manuscript draft: ZW, TCC, LC. All authors have read and approved the final version of the manuscript.
Conflicts of interest:The authors declare that they have no conflict of interest.
Data availability statement:The data that support the findings of this study are openly available in Figshare database at http://doi.org/10.6084/m9.figshare.21405477 (CpG_methylation_annotations_stages), http://doi.org/ 10.6084/m9.figshare.21405483 (CHG_methylation_annotations_stages), http://doi.org/.10.6084/m9.figshare.21405486 (CHH_methylation_annotations_stages), http://doi.org/10.6084/m9.figshare.21401466 (CpG_methylation_annotations_all_24), http://doi.org/10.6084/m9.figshare.21401523 (CHG_methylation_annotation_all_24), http://doi.org/., http://doi.org/10.6084/m9.figshare.21401529 (CHH_methylation_annotations_all_24).
Open peer reviewers:Ilaria Palmisano, Imperial College London, UK; Elena Giusto, San Camillo Hospital, Italy.
Additional Table 1: Gene list of differentially methylated regions within promoter.
Additional Table 2: The amount, distribution and context of methylated cytosine through the whole genome after spinal cord injury.
Additional Table 3: Detailed statistical information of CHH, CPG, and CHH promoter methylation with eight times annotated.
Additional Figure 1: The expression levels of DNA methyltransferases, hydroxymethyltransferases and methylation binding-proteins were associated with CpG and CH methylation post-spinal cord injury.
Additional Figure 2: Percentage of methylation marks per base in each sample.
Additional Figure 3: The similarity and clustering of CHG and CHH contexts with methylation patterns.
Additional Figure 4: The relationship of cytosine density/position and methylation levels in the whole genome.
Additional Figure 5: The DMRs changes at non-coding region post-SCI and the GO analysis of CHG promoter DMRs.
Additional Figure 6: Promoter DMRs of CHH in comparison between each time group and sham group.
Additional Figure 7: The volcano plots demonstrated pairwise comparisons among the promoter DMRs at different stages.
Additional Figure 8: The Locus overlap for genomic region sets and regulatory elements.
Additional Figure 9: The GO and KEGG analysis for temporal changes of promoter DMGs in CHG context.
Additional Figure 10: The GO and KEGG analysis for temporal changes of promoter DMGs in CHH context.
Additional file 1: Open peer review reports 1 and 2.
P-Reviewers: Palmisano I, Giusto E; C-Editor: Zhao M; S-Editors: Yu J, Li CH; L-Editors: Yu J, Song LP; T-Editor: Jia Y
Additional Table 1
Gene list of differentially methylated regions within promoter.
Additional Table 3
Detailed statistical information of CHH, CPG, and CHH promoter methylation with eight times annotated.
The authors are thankful to Dr. Yu Tao at Department of Molecular, Cell and Developmental Biology, University of California, Los Angeles for her kind assistance in establishing of RRBS library; and Dr. Simon Andrews at Babraham Institute for his help in analysis with the SeqMonk pipeline.
1. Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, Mason CE 2012. methylKit:a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol 13:R87.
2. Boni JL, Kahanovitch U, Nwaobi SE, Floyd CL, Olsen ML 2020. DNA methylation:a mechanism for sustained alteration of KIR4.1 expression following central nervous system insult. Glia 68:1495–1512.
3. Booth MJ, Branco MR, Ficz G, Oxley D, Krueger F, Reik W, Balasubramanian S. 2012 Quantitative sequencing of 5-methylcytosine and 5-hydroxymethylcytosine at single-base resolution. Science 336:934–937.
4. Celarain N, Tomas-Roig J 2020. Aberrant DNA methylation profile exacerbates inflammation and neurodegeneration in multiple sclerosis patients. J Neuroinflammation 17:21.
5. Chestnut BA, Chang Q, Price A, Lesuisse C, Wong M, Martin LJ 2011. Epigenetic regulation of motor neuron cell death through DNA methylation. J Neurosci 31:16619–16636.
6. Crunkhorn S 2019. Epigenetic route to recovery for spinal cord injury. Nat Rev Drug Discov 18:420.
7. Danilov CA, Gu Y, Punj V, Wu Z, Steward O, Schönthal AH, Tahara SM, Hofman FM, Chen TC 2020. Intravenous delivery of microRNA-133b along with Argonaute-2 enhances spinal cord recovery following cervical contusion in mice. Spine J 20:1138–1151.
8. Davaa G, Hong JY, Kim TU, Lee SJ, Kim SY, Hong K, Hyun JK 2021. Exercise ameliorates spinal cord injury by changing DNA methylation. Cells 10:143.
9. Deaton AM, Bird A 2011. CpG islands and the regulation of transcription. Genes Dev 25:1010–1022.
10. Du H, Shi J, Wang M, An S, Guo X, Wang Z 2018. Analyses of gene expression profiles in the rat dorsal horn of the spinal cord using RNA sequencing in chronic constriction injury rats. J Neuroinflammation 15:280.
11. Duan K, Mayer AR, Shaff NA, Chen J, Lin D, Calhoun VD, Jensen DM, Liu J 2021. DNA methylation under the major depression pathway predicts pediatric quality of life four-month post-pediatric mild traumatic brain injury. Clin Epigenetics 13:140.
12. Fernandez AF, Assenov Y, Martin-Subero JI, Balint B, Siebert R, Taniguchi H, Yamamoto H, Hidalgo M, Tan AC, Galm O, Ferrer I, Sanchez-Cespedes M, Villanueva A, Carmona J, Sanchez-Mut JV, Berdasco M, Moreno V, Capella G, Monk D, Ballestar E, et al. 2012 A DNA methylation fingerprint of 1628 human samples. Genome Res 22:407–419.
13. Gasparoni G, Bultmann S, Lutsik P, Kraus TFJ, Sordon S, Vlcek J, Dietinger V, Steinmaurer M, Haider M, Mulholland CB, Arzberger T, Roeber S, Riemenschneider M, Kretzschmar HA, Giese A, Leonhardt H, Walter J 2018. DNA methylation analysis on purified neurons and glia dissects age and Alzheimer's disease-specific changes in the human cortex. Epigenetics Chromatin 11:41.
14. Haghighi F, Ge Y, Chen S, Xin Y, Umali MU, De Gasperi R, Gama Sosa MA, Ahlers ST, Elder GA 2015. Neuronal DNA methylation profiling of blast-related traumatic brain injury. J Neurotrauma 32:1200–1209.
15. Hahn MA, Li AX, Wu X, Pfeifer GP. 2015 Single base resolution analysis of 5-methylcytosine and 5-hydroxymethylcytosine by RRBS and TAB-RRBS. Methods Mol Biol 1238:273–287.
16. Hsu FM, Yen MR, Wang CT, Lin CY, Wang CR, Chen PY 2017. Optimized reduced representation bisulfite sequencing reveals tissue-specific mCHH islands in maize. Epigenetics Chromatin 10:42.
17. Hutson TH, Kathe C, Palmisano I, Bartholdi K, Hervera A, De Virgiliis F, McLachlan E, Zhou L, Kong G, Barraud Q, Danzi MC, Medrano-Fernandez A, Lopez-Atalaya JP, Boutillier AL, Sinha SH, Singh AK, Chaturbedy P, Moon LDF, Kundu TK, Bixby JL, et al. 2019. Cbp-dependent histone acetylation mediates axon regeneration induced by environmental enrichment in rodent spinal cord injury models. Sci Transl Med 11:eaaw2064.
18. Krueger F, Andrews SR. 2011 Bismark:a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27:1571–1572.
19. Li C, Zhu X, Lee CM, Wu Z, Cheng L 2020. A mouse model of complete-crush transection spinal cord injury made by two operations. Ann Transl Med 8:210.
20. Li C, Wu Z, Zhou L, Shao J, Hu X, Xu W, Ren Y, Zhu X, Ge W, Zhang K, Liu J, Huang R, Yu J, Luo D, Yang X, Zhu W, Zhu R, Zheng C, Sun YE, Cheng L 2022. Temporal and spatial cellular and molecular pathological alterations with single-cell resolution in the adult spinal cord after injury. Signal Transduct Target Ther 7:65.
21. Li E, Beard C, Jaenisch R. 1993 Role for DNA methylation in genomic imprinting. Nature 366:362–365.
22. Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery JR, Lee L, Ye Z, Ngo QM, Edsall L, Antosiewicz-Bourget J, Stewart R, Ruotti V, Millar AH, Thomson JA, Ren B, Ecker JR. 2009 Human DNA methylomes at base resolution show widespread epigenomic differences. Nature 462:315–322.
23. Lister R, Mukamel EA, Nery JR, Urich M, Puddifoot CA, Johnson ND, Lucero J, Huang Y, Dwork AJ, Schultz MD, Yu M, Tonti-Filippini J, Heyn H, Hu S, Wu JC, Rao A, Esteller M, He C, Haghighi FG, Sejnowski TJ, et al. 2013. Global epigenomic reconfiguration during mammalian brain development. Science 341:1237905.
24. Liu NK, Wang XF, Lu QB, Xu XM. 2009 Altered microRNA expression following traumatic spinal cord injury. Exp Neurol 219:424–429.
25. Mahar M, Cavalli V. 2018 Intrinsic mechanisms of neuronal axon regeneration. Nat Rev Neurosci 19:323–337.
26. Moore LD, Le T, Fan G. 2013 DNA methylation and its basic function. Neuropsychopharmacology 38:23–38.
27. Ramsahoye BH, Biniszkiewicz D, Lyko F, Clark V, Bird AP, Jaenisch R. 2000 Non-CpG methylation is prevalent in embryonic stem cells and may be mediated by DNA methyltransferase 3a. Proc Natl Acad Sci U S A 97:5237–5242.
28. Sabbagh MF, Heng JS, Luo C, Castanon RG, Nery JR, Rattner A, Goff LA, Ecker JR, Nathans J 2018. Transcriptional and epigenomic landscapes of CNS and non-CNS vascular endothelial cells. Elife 7:e36187.
29. Shearstone JR, Pop R, Bock C, Boyle P, Meissner A, Socolovsky M 2011. Global DNA demethylation during mouse erythropoiesis in vivo. Science 334:799–802.
30. Shi GD, Zhang XL, Cheng X, Wang X, Fan BY, Liu S, Hao Y, Wei ZJ, Zhou XH, Feng SQ. 2018 Abnormal DNA methylation in thoracic spinal cord tissue following transection injury. Med Sci Monit 24:8878–8890.
31. Smith ZD, Chan MM, Mikkelsen TS, Gu H, Gnirke A, Regev A, Meissner A. 2012 A unique regulatory phase of DNA methylation in the early mammalian embryo. Nature 484:339–344.
32. Song HH, Song TC, Yang T, Sun CS, He BQ, Li H, Wang YJ, Li Y, Wu H, Hu YM, Wang YJ. 2021 High mobility group box 1 mediates inflammatory response of astrocytes via cyclooxygenase 2/prostaglandin E2 signaling following spinal cord injury. Neural Regen Res 16:1848–1855.
33. Sugeno A, Piao W, Yamazaki M, Takahashi K, Arikawa K, Matsunaga H, Hosokawa M, Tominaga D, Goshima Y, Takeyama H, Ohshima T. 2021 Cortical transcriptome analysis after spinal cord injury reveals the regenerative mechanism of central nervous system in CRMP2 knock-in mice. Neural Regen Res 16:1258–1265.
34. Sun ZC, Liang F, Yang J, Hai Y, Su QJ, Liu XH. 2022 The mechanism by which hyperbaric oxygen treatment alleviates spinal cord injury:genome-wide transcriptome analysis. Neural Regen Res 17:2737–2742.
35. Suzuki MM, Bird A. 2008 DNA methylation landscapes:provocative insights from epigenomics. Nat Rev Genet 9:465–476.
36. Thurman D, Sniezek J, Johnson D, Greenspan A, Smith S 1995. Guidelines for surveillance of central nervous system injury. Atlanta, GA: US Department of Health and Human Services, Centers for Disease Control and Prevention.
37. Wu H, Zhang Y 2011. Mechanisms and functions of Tet protein-mediated 5-methylcytosine oxidation. Genes Dev 25:2436–2452.
38. Wu Z, Huang K, Yu J, Le T, Namihira M, Liu Y, Zhang J, Xue Z, Cheng L, Fan G. 2012 Dnmt3a regulates both proliferation and differentiation of mouse neural stem cells. J Neurosci Res 90:1883–1891.
39. Yao B, Jin P. 2014 Cytosine modifications in neurodevelopment and diseases. Cell Mol Life Sci 71:405–418.
40. Yong WS, Hsu FM, Chen PY 2016. Profiling genome-wide DNA methylation. Epigenetics Chromatin 9:26.
41. York EM, Petit A, Roskams AJ 2013. Epigenetics of neural repair following spinal cord injury. Neurotherapeutics 10:757–770.
42. Zhang BY, Chang PY, Zhu QS, Zhu YH. 2020 Decoding epigenetic codes:new frontiers in exploring recovery from spinal cord injury. Neural Regen Res 15:1613–1622.
43. Zhang ZY, Zhang Z, Fauser U, Schluesener HJ. 2007 Global hypomethylation defines a sub-population of reactive microglia/macrophages in experimental traumatic brain injury. Neurosci Lett 429:1–6.
44. Ziller MJ, Müller F, Liao J, Zhang Y, Gu H, Bock C, Boyle P, Epstein CB, Bernstein BE, Lengauer T, Gnirke A, Meissner A 2011. Genomic distribution and inter-sample variation of non-CpG methylation across human cell types. PLoS Genet 7:e1002389.