DNA methylation confers a cerebellum-specific identity in non-human primates

INTRODUCTION

The mammalian brain is an exceedingly complex organ that performs many unique neuronal activities in distinct brain regions, such as cognition in the prefrontal cortex, memory consolidation in the hippocampus, and motor coordination in the cerebellum. Tremendous efforts using molecular genetic tools have identified multiple transcription factors and regulatory networks that govern these region-specific functions. In parallel, epigenetic modifications are increasingly recognized as crucial regulators of cellular identity and adaptive responses to environmental stimuli. Despite their fundamental role in modulating gene expression, the relationship between epigenetic modifications and gene expression in the context of brain subdivision identity remains largely underexplored.

DNA methylation serves as a repressive epigenetic mechanism at the DNA level. Distinct DNA methylation patterns are established across tissues, shaping transcriptional landscapes that define tissue-specific functions (Ghosh et al., 2010; Shiota et al., 2002; Suzuki et al., 2007). Increasing evidence suggests that DNA methylation plays a pivotal role in brain development and pathology by influencing gene regulation and cellular differentiation (Altuna et al., 2019; De Souza et al., 2016; Feng & Fan, 2009; Gopalakrishnan et al., 2008; Robertson & Wolffe, 2000), dynamically regulated by various catalytic enzymes and reader proteins. In mammalian genomes, DNA methylation occurs at the 5-position of cytosine (5-methycytosine, 5mC) and is catalyzed by DNA methyltransferases (DNMTs), which enzymatically transfer a methyl group onto cytosine residues (Lyko et al., 1999; Okano et al., 1998). CpG sites often cluster to form CpG islands, which are frequently located within gene promoters. Unmethylated CpG islands are generally associated with transcriptional activation (Bird, 2002; Li, 2002), whereas extensive 5mC enrichment at promoters, CpG island shores, and enhancers is linked to transcriptional repression (Deaton & Bird, 2011). DNA methylation is thought to contribute to the maintenance of tissue identity by regulating genes essential for tissue-specific functions (Blake et al., 2020; Cui et al., 2020; Wan et al., 2015).

The ten-eleven translocation (TET) enzyme family, comprising TET1, TET2, and TET3, catalyzes the conversion of 5mC to 5-hydroxymethycytosine (5hmC) (Cimmino et al., 2011; Koh et al., 2011). Beyond serving as an intermediate in DNA demethylation, 5hmC is increasingly recognized as a stable component of DNA with a regulatory role in gene expression (Mellén et al., 2012; Song et al., 2013). In human embryonic and mouse stem cells, 5hmC is preferentially enriched within 0.5–2 kb upstream and downstream of the transcription start site (TSS) and at the border of CpG island promoters (Ehrlich & Ehrlich, 2014). The accumulation of 5hmC within gene bodies is generally associated with transcriptional activation (Cui et al., 2020). Dysregulation of DNA methylation machinery, including proteins responsible for establishing, maintaining, and removing methylation, has been implicated in various neurological disorders (Li & Liu, 2011; Tan & Shi, 2012), underscoring the potential importance of DNA methylation in neural development and function. Our previous study identified a cerebellum-specific 5hmC distribution pattern in monkey brains, suggesting its importance in establishing and maintaining cerebellar identity and tissue-specific functions (Xu et al., 2022). Given the essential interplay between 5mC and 5hmC in epigenetic regulation (Kato & Iwamoto, 2014; Li et al., 2019; Ponnaluri et al., 2017), the current study systematically profiled the 5mC and 5hmC landscapes across multiple brain regions in rhesus monkeys to elucidate their functional significance in primate neurobiology. Our findings revealed a distinct cerebellum-specific coordination between 5mC and 5hmC, along with the selective expression of schizophrenia-associated genes in the cerebellum. These results reinforce the notion that DNA methylation modifications are intricately coordinated to govern gene expression and shape cerebellar identity.

MATERIALS AND METHODS

Animal tissues

Brain tissues were obtained from male rhesus monkeys of varying ages during previous studies on genetically modified monkey models (Yang et al., 2019, 2015; Yin et al., 2019). Sample collection and handling followed the same procedures as those used for 5hmC profiling (Xu et al., 2022). All animal procedures were approved by the Institutional Animal Care and Use Committee at Guangdong Landao Biotechnology Co. Ltd., Guangzhou (Approval number LDACU 20220225-02).

Genomic DNA and RNA isolation

Brain tissues were homogenized on ice and digested overnight at 55°C with proteinase K (0.667 μg/μL, Sigma Aldrich, USA) in a lysis buffer containing 100 mmol/L Tris-HCl (pH 8.5), 5 mmol/L EDTA, 0.2% SDS (vol/vol), and 200 mmol/L NaCl). RNAse (40 μg) was then added and incubated at 37°C for 1 h. DNA extraction was performed using phenol:chloroform:isoamyl alcohol (25:24:1) saturated with 10 mmol/L Tris (pH 8.0) and 1 mmol/L EDTA (Sigma Aldrich, USA). Lysates were centrifuged at 12 000 r/min (centrifugal radius of 8 cm) for 10 min at room temperature, and the aqueous phase was precipitated with an equal volume of isopropanol. The resulting pellet was washed with 75% ethanol, air-dried, and resuspended in nuclease-free water. Total RNA was purified from each tissue using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions. DNA and RNA concentrations were measured using a Qubit 3.0 fluorometer (Invitrogen, USA).

Methylated DNA immunoprecipitation (MeDIP)

Genome-wide 5mC profiling was conducted using MeDIP (Szulwach et al., 2010). Briefly, 5 µg of genomic DNA was sonicated to approximately 200 bp using a Bioruptor in 1× TE buffer. Sheared DNA underwent end repair, dA-tailing and adaptor ligation using the NEBNext® Ultra™ II DNA Library Prep Kit (NEB, USA) per the manufacturer’s protocols. The DNA was denatured for 10 min at 100°C in a dry heat block, then immediately placed on ice for 10 min. Immunoprecipitation (IP) was performed in 600 µL of IP buffer (100 mmol/L Tris-HCL, pH 7.4, 150 mmol/L NaCl, and 0.05% Triton X-100) with 1 µg of 5-methycytidine antibody (BI-MECY-0500; Eurogentec, Belgium) or 1 µg of normal mouse IgG (Millipore, USA) at 4°C with rotation overnight. The antibody-DNA complexes were pulled down using 50 µL of Dynabeads Protein G (Invitrogen, USA) incubated at 4°C for 2 h with rotation. The beads were collected with a 1.5 mL microcentrifuge tube holder magnet, washed three times in IP buffer at room temperature for 10 min each with rotation, and resuspended in 100 µL of proteinase K digestion buffer (1× TE and 0.25% SDS and 2.5 mg/mL proteinase K). Digestion was performed at 55°C for 2 h in a thermomixer (Eppendorf, Germany) set at 1 000 r/min (centrifugal radius of 8 cm). Following digestion, beads were magnetically separated, and DNA was extracted from the supernatant using a MinElute PCR Purification Kit (Qiagen, USA). The eluted DNA, used for polymerase chain reaction (PCR) analysis, was further processed using AMPure beads (Beckman Coulter, USA) before quantification with an Agilent 2100 Bioanalyzer (USA). High-quality libraries were sequenced on the Illumina HiSeq X10 platform (150 bp paired-end reads) by Shenzhen Acegen Technology Co. Ltd. (China).

Reverse transcription-quantitative real-time PCR (RT-qPCR)

Total RNA (1 μg) was reverse-transcribed into cDNA using a PrimeScript^TM RT Reagent Kit (TaKaRa, Japan) according to the manufacturer’s instructions. Quantitative PCR (qPCR) was performed in triplicate for each sample using the QuantiNova SYBR Green PCR Kit (Qiagen, USA) in 20 µL reactions. Amplifications were carried out on a Bio-Rad CFX Maestro 1.0 system under the following conditions: initial denaturation at 95°C for 2 min, followed by amplification at 95°C for 10 s and 60°C for 30 s for 40 cycles. Gene expression levels were normalized to RPL32 signals and fold enrichment was calculated using 2^−dCt, where dCt=Ct (5-hmC enriched)−Ct (input). The results were plotted with GraphPad Prism.

Western blot analysis

Monkey brain tissues were homogenized in ice-cold RIPA buffer supplemented with a protease inhibitor cocktail. The lysates were incubated on ice at 4°C for 30 min, sonicated, and centrifuged at 12 000 r/min (centrifugal radius of 8 cm) for 10 min at 4°C. Supernatants containing total protein (30 µg per sample) were separated using sodium dodecyl-sulfate polyacrylamide gel electrophoresis (SDS-PAGE), then transferred to polyvinylidene fluoride (PVDF) membranes (Millipore, USA). The membranes were blocked with 5% milk in phosphate-buffered saline (PBS) for 1 h at room temperature, then incubated at 4°C with primary antibodies diluted in 3% bovine serum albumin (BSA) and PBS overnight. After three washes in PBS, the membranes were incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies in 5% milk and PBS for 1 h at room temperature. Following three washes in PBST (0.1% Tween-20), protein signals were detected using enhanced chemiluminescence (ECL Prime, GE Healthcare, USA).

Immunohistochemistry

Monkey brains were perfused with PBS, followed by 4% paraformaldehyde in PBS at 4°C overnight. Fixed tissues were cryoprotected in 30% sucrose for 48 h, then sectioned into 20 or 40 μm slices using a cryostat (Leica CM1850, Germany) at −20°C. Tissue sections were blocked in 4% donkey serum with 0.2% Triton X-100 and 3% BSA in PBS for 1 h. Sections (20 μm) were incubated with primary antibodies in the same buffer at 4°C overnight. After washing with 1× PBS, the sections were incubated in fluorescent secondary antibodies. Fluorescent images were acquired using a Zeiss microscope (Carl Zeiss Imaging, Axiovert 200 MOT, Germany) equipped with either a 40× or 63× objective lens (LD-Achroplan 40×/0.6 or 63×/0.75). Images were captured using a digital camera (Hamamatsu, Orca-100, Japan) and analyzed using Openlab software v.5.0.2 (Improvision, UK).

Target bisulfite sequencing (TBS)

Genomic DNA was subjected to bisulfite conversion using the MethylCode^TM Bisulfite Conversion Kit (Invitrogen, USA) following the supplier’s instructions. Target-specific primers for TBS amplification were designed using the MethPrimer online tool (https://www.urogene.org/cgi-bin/methprimer/methprimer.cgi). Bisulfite-treated DNA was then used as a template for PCR amplification with TaKaRa EpiTaq™ HS (TaKaRa, Japan). PCR amplicons were purified using Ampure XP beads and eluted in 20 μL of H₂O. DNA concentrations were quantified using a Qubit High Sensitivity Kit (Invitrogen, USA), and equimolar amounts of amplicons from each sample were pooled for library preparation. Pooled libraries underwent MiSeq deep sequencing (100× or above) following standard procedures recommended by Illumina. Image analysis and base calling were performed using standard Illumina pipelines. BSMAP v.2.7.3 (Xi & Li, 2009) was used to map target sequences and detect methylation levels at various sites on the amplicon. Methylated cytosines at CpG sites were identified, and a t-test was applied to detect each amplicon site in each comparison group.

RNA sequencing (RNA-seq)

RNA libraries were constructed using the NEBNext® Ultra™ II Directional RNA Library Prep Kit (NEB, USA) for Illumina following the manufacturer’s instructions. Library quality and concentration were assessed using an Agilent 2100 Bioanalyzer. High-quality libraries were sequenced on the Illumina HiSeq X10 platform (150 bp paired-end reads; USA) by Shenzhen Acegen Technology Co. Ltd. (China).

Sequencing read quality control and alignment

5hmC sequencing data were previously published (Xu et al., 2022). Demultiplexed sequencing reads from MeDIP-seq, BS-seq, and RNA-seq libraries were filtered and trimmed to remove the low-quality bases and adaptor sequences using Trimmomatic v.0.36 (Bolger et al., 2014) with the following parameters: adapter.fa:2:30:7:1:TRUE; LEADING:3; SLIDINGWINDOW: 4:15; TRAILING:3; and MINLEN:36. For the MeDIP libraries, trimmed reads were aligned to the rhesus macaque reference genome (rheMac8) using Burrows-Wheeler Aligner (BWA) software with default parameters. After filtering duplicate reads, mapped read pairs were extracted and converted to BED format for downstream analysis.

Visualization of read densities

To visualize genome-wide read densities, MeDIP-seq read coverage was normalized to one million reads per library and computed at 50 bp intervals across the genome. Normalized read density profiles were stored in Wiggle (WIG) format, and genomic region screenshots were taken using the Integrative Genomics Viewer (IGV) genome browser.

Identification of enriched regions and tissue- and age-specific DMRs

DMRs were identified using MACS (v.1.4.9) for peak calling, comparing 5mC-enriched regions against input controls for each sample. Tissue- and age-specific DMRs were detected using PePr (v.1.1.1.18) with the following parameters: --file-format bed --peaktype sharp --difftest TRUE --keep-max-dup 2-threshold 1e-5-normalization inter-group. PePr uses a sliding window approach and models read count distributions across replicates and between groups using a negative binomial framework (Zhang et al., 2014). To examine gene-level methylation patterns, hypermethylated and hypomethylated DMR-associated genes were visualized using Venn diagrams generated in R. DMR genes related to cerebellar disorders were analyzed using the Cytoscape platform.

RESULTS

Cerebellum-specific 5mC pattern in rhesus monkey brains

Dynamic 5mC states have been reported across different biological systems, including the brain and blood (Hackett & Surani, 2013; Luo et al., 2018; Suelves et al., 2016). To characterize region-specific DNA methylation patterns in the primate brain, genome-wide 5mC distribution was analyzed in four brain regions, including the cortex, cerebellum, striatum and hippocampus, across three age groups: juvenile (2 years), young (8 years), and old (17 years). MeDIP-seq was employed to enrich for genome-wide 5mC, and high-quality sequencing reads were aligned to the rhesus monkey reference genome (rheMac2) using BWA software. Quality assessment confirmed high sequencing accuracy, with Q20>97% and Q30>93% (Supplementary Table S1). Genome-wide 5mC distribution was examined within 2 kb upstream and downstream of the TSS and transcriptional end site (TES). Notably, 5mC was enriched near the TSS region and reduced near the TES region. Furthermore, 5mC accumulation was more pronounced in older monkeys compared to juvenile and young monkeys (Supplementary Figure S1, turquoise line), resembling the age-dependent increase in 5hmC reported previously (Xu et al., 2022). Additionally, 5mC levels were relatively low upstream of TSS, where transcription initiation predominantly occurs, but were elevated upstream of TES, exceeding levels observed at TSS (Supplementary Figure S1). These findings suggest a potential regulatory role for 5mC in transcription termination.

Our previous study revealed a distinct age-dependent 5hmC distribution in the cerebellum of rhesus monkey brains (Xu et al., 2022). Here, principal component analysis (PCA) of genome-wide 5mC mapping from the analyzed brain regions identified a unique clustering pattern in the cerebellum (Figure 1A), closely resembling the previously observed 5hmC distribution. This distinct 5mC pattern was established early in development and remained stable across aging (Figure 1A). Given that 5mC serves as the direct precursor to 5hmC, its enrichment in the cerebellum may contribute to the formation of the cerebellum-specific 5hmC landscape.

Figure  1.  Cerebellum-specific 5mC patterns and DMRs in rhesus monkey brains

A: PCA of 5mC in all samples and in 2-, 8-, and 17-year-old groups. B: Cerebellum-specific DMRs in 2-, 8-, and 17-year-old groups. DMRs for the three groups were 8 497, 26 446, and 40 354, respectively. C: Volcano plots and GO functional analysis of genes with cerebellum-specific DMRs in gene body regions. D: Volcano plots and GO functional analysis of genes with cerebellum-specific DMRs in promoter regions. Cere: Cerebellum; Ctx: Cortex; Hipp: Hippocampus; Stri: Striatum.

DownLoad: Full-Size Img PowerPoint

Age-dependent increase in cerebellum-specific DMRs and their enrichment in metabolic pathways

Genome-wide 5mC profiling revealed both age-associated 5mC accumulation and a distinct cerebellum-specific methylation pattern. To further categorize these epigenetic dynamics, DMRs were identified in the cerebellum and compared with those in other regions (Figure 1B). The number of cerebellum-specific DMRs increased substantially with age, reaching 8 497, 26 446, and 40 354, in juvenile, young, and old monkeys, respectively (Table 1). This trend is consistent with our previous findings, which showed the highest number of DhMRs in the cerebellum of aged monkeys (Xu et al., 2022). Given that 5mC serves as the direct precursor of 5hmC, these results suggest that 5mC and 5hmC are coordinately regulated during aging in the rhesus monkey brain. To further characterize cerebellum-specific DMRs, hypermethylated and hypomethylated genes were quantified, followed by Gene Ontology (GO) analysis. Results showed that most genes were in a hypermethylated state. Notably, within gene bodies, 490 genes were hypermethylated and 252 genes were hypomethylated (Figure 1C). In the promoter regions, 628 genes were hypermethylated and 155 genes were hypomethylated (Figure 1D). Functional enrichment analysis indicated that most DMRs were enriched in metabolic pathways, including primary and macromolecular metabolism. In particular, DMRs were markedly enriched in pathways related to gene expression regulation (Figure 1C, D). Overall, these findings suggest that dynamic cerebellum-specific DMRs may contribute to the aging process by modulating gene expression and metabolic activity.

Table  1.  Cerebellum-specific and brain-specific DMRs in three groups

DMRs	Cerebellum-specific	Brain-specific
Juvenile	8497	12293
Young	26446	16551
Old	40354	15859

 | Show Table

DownLoad: CSV

Distinct gene expression patterns in the cerebellum

DNA methylation plays a crucial role in gene regulation, with extensive evidence supporting the impact of 5mC on transcriptional activity (Gibney & Nolan, 2010; Wu & Morris, 2001). To comprehensively evaluate the relationship between DNA methylation and gene expression, RNA-seq was performed across four brain regions of rhesus monkeys. Surprisingly, PCA revealed a distinct transcriptional profile in the cerebellum, which was completely separated from the other three regions (Figure 2A). Despite this regional specificity, gene expression patterns did not show age-dependent differences, with a high degree of overlap observed among the three age groups (Figure 2B). These results indicate that the cerebellum possesses a unique transcriptional signature, mirroring the previously observed 5mC and 5hmC distribution patterns, strongly implicating DNA methylation as a key regulatory mechanism in cerebellar identity. To further dissect the transcriptional landscape of the cerebellum, differentially expressed genes (DEGs) were quantified across age groups. The number of down-regulated genes exceeded the number of up-regulated genes in all three groups, with 1 791 down-regulated and 1 041 up-regulated in juveniles, 3 018 down-regulated and 1 608 up-regulated in the young group, and 1 490 down-regulated and 822 up-regulated in the old group (Supplementary Figure S2). Furthermore, GO analysis revealed that DEGs were significantly enriched in pathways related to environmental information processing and organismal systems, particularly signaling transduction and neuronal function (Figure 2C). As differential gene expression patterns underpin tissue identity (Fagerberg et al., 2014; Sonawane et al., 2017), these findings align with the Genotype-Tissue Expression (GTEx) project, which established a strong link between gene expression patterns and tissue specification (GTEx Consortium, 2017). The cerebellum-specific transcriptional signature observed in rhesus monkeys suggests that DNA methylation, through its regulatory influence on gene expression, may play a fundamental role in defining cerebellar identity.

Figure  2.  Unique gene expression patterns in the cerebellum

A: PCA of cerebellum-specific gene expression patterns in rhesus monkey brains. B: PCA of three age groups. C: DEGs were enriched in environmental information processing and organismal systems. Most DEGs were enriched in pathways related to signaling transduction and neuronal functions. Y-axis shows GO terms, x-axis shows significance of P-values and groups. Cere: Cerebellum; Ctx: Cortex; Hipp: Hippocampus; Stri: Striatum.

DownLoad: Full-Size Img PowerPoint

Genomic distribution of 5mC and 5hmC in the primate brain

To further investigate the involvement of DNA methylation in gene regulation, the distribution of 5mC across different genomic components was analyzed in four brain regions. Compared to other regions, the cerebellum exhibited significantly lower 5mC levels in promoters, CpG islands, and gene bodies, while no notable differences were observed in repetitive elements (Figure 3A). The most pronounced reduction in 5mC was found at CpG islands, implying a regulatory function in cerebellum-specific gene expression. Beyond its intermediary role in active demethylation, 5hmC also functions as an epigenetic mark associated with transcriptional activation (Guo et al., 2011; Liu et al., 2022; Williams et al., 2011). Previous studies have demonstrated that 5hmC is enriched in active enhancers and gene bodies and is involved in transcriptional regulation (Cui et al., 2020; Wu et al., 2011). Analysis of 5hmC distribution in different genomic components demonstrated higher abundance in CpG islands and gene bodies in the cerebellum than in the other brain regions, exhibiting an inverse relationship to 5mC abundance at the same regions. However, 5mC and 5hmC displayed the same low abundance at promoters in the cerebellum (Figure 3B).

Figure  3.  Comparison of 5mC and 5hmC distribution in genomic components and cerebellum-specific DhMR enrichment

A: 5mC read density in promoter, CpG island, gene body, and repeat regions in four brain regions. B: 5hmC read density in promoter, CpG island, gene body, and repeat regions in four brain regions. C: Hypermethylated-DhMR and hypomethylated-DhMR genes were enriched in schizophrenia-related pathways. Boxes represent genes related to diseases and cerebellum functions. Hyper-DhMR genes (red) and hypo-DhMR genes (green) are shown. Cere: Cerebellum; Ctx: Cortex; Hipp: Hippocampus; Stri: Striatum.

DownLoad: Full-Size Img PowerPoint

Cerebellum-specific DhMRs are enriched in schizophrenia-related pathways

To explore the functional significance of cerebellum-specific 5hmC clustering, enrichment analysis of DhMRs was performed. Results showed that genes harboring cerebellum-specific hyper- or hypomethylated DhMRs were closely associated with pathological conditions, particularly schizophrenia-related pathways, with 38 genes enriched in these networks (Figure 3C). Integration of RNA-seq data further revealed distinct gene expression patterns among these 38 genes. Among the 16 genes exhibiting hypermethylated DhMRs, 12 (75%) were down-regulated, while among the 22 genes with hypomethylated DhMRs, 14 (64%) were up-regulated (Supplementary Table S2). To validate these transcriptomic findings, qPCR was performed on a subset of genes exhibiting the most significant expression changes, including glutamate metabotropic receptor 4 (GRM4), glutamate ionotropic receptor NMDA type subunit 2C (GRIN2C), amyotrophic lateral sclerosis-2 C-terminal like (ALS2CL), paired box 6 (PAX6), guanine nucleotide-binding protein G subunit beta-3 (GNB3), glial fibrillary acidic protein (GFAP), plasmolipin (PLLP), and cholinergic receptor muscarinic 4 (CHRM4). The qPCR results were consistent with the RNA-seq findings (Figure 4A). To further corroborate these observations at the protein level, western blotting (Figure 4B; Supplementary Figure S3) and immunostaining (Figure 4C) were performed for representative genes, confirming expression trends observed at the transcriptomic level. These findings indicate that cerebellum-specific DhMRs influence the regulation of schizophrenia-related genes, with most showing negative correlations. Collectively, these results support a role for 5hmC in modulating cerebellar functions relevant to neuropsychiatric disorders.

Figure  4.  Schizophrenia-related gene expression in the cerebellum

A: Transcriptional levels of GRM4, GRIN2C, ALS2CL, PAX6, GNB3, GFAP, PLLP, and CHRM4. qPCR was applied to quantify the relative expression of cerebellum-specific genes (n=3). Cere: Cerebellum, Ctx: Cortex, Stri: Striatum, Hipp: Hippocampus. B: Protein levels of GRIN2C, PAX6, ALS2CL, GFAP, and CHRM4 in different brain regions (n=4) based on western blot analysis. C: Immunofluorescence staining of cortex and cerebellum of a 9-year-old rhesus monkey with antibodies against GRIN2C and PAX6. Purkinje cells were immunopositive for CD-28K, granule cells were positive for NeuN. D: 5mC and 5hmC levels in promoter regions of PAX6 and GRIN2C in the cerebellum and other brain regions (referred to here as the brain). Cere: Cerebellum. Data were analyzed by one-way ANOVA, with results presented as mean±SEM.***: P<0.001.

DownLoad: Full-Size Img PowerPoint

Cerebellum-specific DNA methylation correlates with gene expression

Accumulating evidence suggests that 5mC and 5hmC function in a coordinated manner to regulate gene expression (Ponnaluri et al., 2017; Singh & Edwards, 2021). To evaluate the contribution of 5mC and 5hmC to cerebellum-specific gene expression, the correlation between 5mC and 5hmC levels in promoter and gene body regions and gene expression was examined. Venn diagram analysis showed that both 5mC and 5hmC could exhibit either negative or positive correlations with gene expression, depending on their genomic localization (Figure 5A, B). These findings underscore the complexity of 5mC and 5hmC interactions in regulating gene expression, highlighting their context-dependent roles. Further quantitative analysis showed that 5mC was predominantly associated with negative correlations in both promoter and gene body regions, whereas 5hmC showed more positive correlations, particularly within gene bodies (Table 2). Venn diagram analysis identified 8 719 genes where 5mC at both the promoter and gene body was negatively correlated with gene expression (Figure 5A). In contrast, 6 458 genes with 5hmC in both regions were positively correlated with gene expression (Figure 5B).

Table  2.  Correlation between DNA methylation and gene expression

DNA methylation	Positive	Negative
5mC at promoter	9247	11662
5mC at genebody	9589	13013
5hmC at promoter	9363	7305
5hmC at genebody	13717	7600

 | Show Table

DownLoad: CSV

Figure  5.  Coordinated regulation of 5mC and 5hmC on gene expression

A, B: Cerebellum-specific 5mC (A) and 5hmC (B) in promoter or gene body regions that overlapped with gene expression changes. Pos: Positive correlation between DNA methylation level and gene expression level. Neg: Negative correlation between DNA methylation and gene expression level. C: GO functional analysis of 5mC-negatively correlated genes. D: GO functional analysis of 5hmC-positively correlated genes. Y-axis shows GO terms, x-axis shows significance test.

DownLoad: Full-Size Img PowerPoint

To explore the functional significance of these epigenetic modifications, GO analysis was performed on genes negatively correlated with 5mC and positively correlated with 5hmC. Results indicated that a substantial number of genes were enriched in processes related to binding and channel activity (Figure 5C, D). Notably, both 5mC-negative and 5hmC-positive regulatory pathways were associated with fundamental biological processes, including DNA biosynthesis, metabolic regulation, and signal transduction. These findings indicate that 5mC and 5hmC play pivotal regulatory roles in modulating core cellular functions. Furthermore, the results imply that cerebellum-specific differential DNA methylation contributes to the establishment of tissue-specific gene expression patterns, reinforcing its role in defining cerebellar identity.

DNA methylation in the promoter region dominates cerebellum-specific gene expression differences

Transcription factors recognize and bind to promoter regions to initiate transcription, a process that is epigenetically regulated by DNA methylation. DNA methylation affects gene expression by blocking the binding of transcription factors or other proteins involved in the transcription process (Buck-Koehntop & Defossez, 2013; Zhu et al., 2016), leading to the repression of gene expression. Conversely, DNA methylation can also facilitate gene activation by serving as a binding site for transcriptional activators (Dhar et al., 2021; Moore et al., 2013). To investigate the contribution of promoter DNA methylation to cerebellum-specific transcriptional regulation, BS-seq was conducted on 29 genes from the most significantly altered pathways, including neuroactive ligand-receptor interactions, axon guidance, and the MAPK signaling pathway (Table 3). Analysis revealed that up to 70% of genes exhibited a negative correlation between DNA methylation and gene expression. Representative genes, such as fibroblast growth factor 5 (FGF5), FMS-like tyrosine kinase 3 (FLT3), cell division cycle 25 homolog B (CDC25B), ephrin type-A receptor 1 (EPHA1), actin-binding LIM protein 1 (ABLIM1), and mitogen-activated protein kinase kinase kinase 6 (MAP3K6), showed significantly reduced DNA methylation in their promoters, especially within 1 000 bp upstream of the TSS (Figure 6A). In contrast, genes such as dopamine D1 receptor (DRD1), cholinergic receptor muscarinic 5 (CHRM5), and somatostatin receptor 5 (SSTR5), showed the opposite patterns, with elevated methylation levels exceeding 50% within 500 bp upstream of the TSS (Figure 6B). The qPCR results confirmed these expression patterns, further supporting the cerebellum-specific relationship between promoter methylation and gene regulation in a non-human primate (Figure 6A, B). These results highlight a negative correlation between DNA methylation in promoter regions and gene expression in the cerebellum of non-human primates.

Table  3.  Cerebellum-specific expressed genes

Pathways	DGEs
Upregulation
MAPK signaling pathways	FGF5, NTF3, DUSP5, FLT3, CDC25B, MAP3K6, CACANA1A
Axon guidance	EPHA8, PIK3R3, PAK7, EPHA1, ABLIM1, SRGAP2, PARD6A
Downregulation
Neuroactive ligand-receptor interaction	DRDRR1, CHRM1, CHRM3, CHRM4, CHRM5, SSTR5, APLNR, HTR6
Axon guidance	NGEF, EPHB6, CAMK2A, MET, SLIT1, SEMA5B, TRPC4

 | Show Table

DownLoad: CSV

Figure  6.  DNA methylation in promoter regions of cerebellum-specific genes

A: Reduced DNA methylation 1 000 bp upstream of promoter was associated with higher expression of cerebellum-specific genes. B: Increased DNA methylation 1 000 bp upstream of promoter was associated with lower expression of cerebellum-specific genes. X-axis represents promoter position. TSS site was 0. Y-axis indicates DNA methylation level. qPCR results confirmed cerebellum-specific gene expression. Data were analyzed by one-way ANOVA, with results presented as mean±SEM. ^*: P<0.05; ^**: P<0.01; ^***: P<0.001; ^****: P<0.0001. Cere: Cerebellum, Brain: other brain regions, Ctx: Cortex, Stri: Striatum, Hipp: Hippocampus.

DownLoad: Full-Size Img PowerPoint

To further validate the negative correlation between DNA methylation and gene expression, the promoter methylation status of PAX6 and GRIN2C was examined. These genes were highly expressed in the cerebellum, and analysis revealed lower 5mC levels in the cerebellum compared to other brain regions, whereas 5hmC levels were elevated (Figure 4D). These findings highlight the negative correlation between DNA methylation and gene expression and strongly support the dominant role of 5mC-mediated transcriptional repression in defining cerebellar gene expression patterns.

Coordinated regulation of 5hmC in promoter regions

Since traditional BS-seq cannot distinguish between 5mC and 5hmC, oxidative BS-seq (oxBS-seq) was employed to determine the specific contribution of 5hmC to gene expression. OxBS-seq introduces an oxidation step that converts 5hmC to 5-formylcytosine (5fC), which is subsequently deaminated by BS treatment and ultimately converted to uracil. This approach enables precise discrimination between 5hmC and 5mC. Analysis of promoter regions revealed that 5hmC levels were consistently lower than 5mC levels, regardless of cerebellum-specific gene expression patterns (Figure 7). Interestingly, two distinct modes of coordinated regulation between 5hmC and 5mC were identified. In the first mode, 5hmC and 5mC exhibited inverse dynamics. As shown in Figure 7A, 5mC levels increased continuously from 1 000 bp upstream to the TSS, with significantly higher levels in the cerebellum compared to other brain regions. Conversely, 5hmC levels decreased in the same promoter regions, reaching the lowest levels in the cerebellum. This pattern was associated with reduced gene expression in the cerebellum, affecting genes such as CHRM1, CHRM3, CHRM4, and semaphorin 5B (SEMA5B). In the second mode, 5hmC and 5mC showed the same dynamic trend, both decreasing from 1 000 bp upstream to the TSS (Figure 7B). This pattern correlated with increased gene expression in the cerebellum, as observed in dual specificity phosphatase 5 (DUSP5) and neurotrophin 3 (NTF3) (Figure 7B). These findings indicate that the regulation of gene expression by 5mC and 5hmC is complex and occurs in the same or opposite directions, depending on the genomic context and promoter architecture.

Figure  7.  Two modes of coordinated regulation of gene expression by 5mC and 5hmC

A: 5mC and 5hmC regulated genes in opposite directions, accompanied by lower gene expression. B: 5mC and 5hmC regulated genes in the same direction, accompanied by higher gene expression. X-axis represents promoter position. TSS site was 0. Y-axis indicates 5mC or 5hmC levels. qPCR results confirmed cerebellum-specific gene expression. Data were analyzed by one-way ANOVA, with results presented as mean±SEM. ^*: P<0.05; ^**: P<0.01; ^***: P<0.001; ^****: P<0.0001. Cere: Cerebellum, Brain: other brain regions, Ctx: Cortex, Stri: Striatum, Hipp: Hippocampus.

DownLoad: Full-Size Img PowerPoint

DISCUSSION

DNA methylation is a fundamental epigenetic modification that regulates gene expression and cellular function, shaping tissue-specific transcriptional landscapes (Greenberg & Bourc'his, 2019; Smith & Meissner, 2013). Characterizing DNA methylation patterns across different tissues and cell types is essential for understanding gene regulation in diverse biological contexts. The present study identified a unique 5hmC landscape in the cerebellum of non-human primates, alongside cerebellum-specific 5mC patterns and transcriptional signatures. Notably, cerebellum-specific genes harboring DhMRs were strongly associated with schizophrenia, a disorder in which cerebellar dysfunction has been implicated. Further analysis revealed that 5mC within 1 000 bp upstream of promoter regions played a dominant role in shaping cerebellum-specific gene expression, with two distinct modes of coordinated regulation by 5hmC. Collectively, these findings suggest that DNA methylation, through both 5mC and 5hmC, contributes to cerebellar identity by regulating genes essential for cerebellar function and disease susceptibility.

Epigenetic modifications are highly dynamic and responsive to environmental influences, with multiple studies on human brain DNA methylation—both pre- and postmortem—revealing substantial variability. However, acquiring fresh human brain samples that meet stringent experimental criteria for regional and age-specific comparisons remains a significant challenge. Non-human primates serve as invaluable models for studying aging and disease mechanisms due to their high genetic, anatomical, and physiological similarity to humans. Specifically, their brain architecture closely resembles that of humans, making them well-suited for neuroscience research. In our previous study, we demonstrated significant epigenetic conservation between monkeys and humans (Xu et al., 2022). The identification of cerebellum-specific DNA methylation patterns in rhesus monkeys and the coordinated regulation of gene expression by 5mC and 5hmC provide new insights into the epigenetic regulation of human cerebellar function and disease pathology. However, challenges remain in directly modeling human neurological disorders using non-human primates, limiting the ability to fully elucidate the relationship between DNA methylation and specific disease phenotypes.

Different tissues in the body acquire specialized identities through unique cellular compositions and gene expression patterns (Breschi et al., 2020; Germanguz et al., 2016). These identities are established during development through a complex interplay of genetic and epigenetic mechanisms (Barrero et al., 2010; Basu & Tiwari, 2021; Berdasco & Esteller, 2010). Epigenetic mechanisms, which play a significant role in determining tissue identity by controlling gene expression, are typically stable and heritable, enabling cells to preserve their tissue-specific characteristics throughout development into adulthood. DNA methylation, a dynamic epigenetic process, plays a crucial role in the establishment and maintenance of tissue identity via gene expression regulation (Bird, 2002). Whole-genome 5hmC profiling across 19 human tissues revealed that tissue-specific DhMRs are enriched in tissue-specific transcription factors, suggesting 5hmC is a fundamental regulator affecting tissue-specific gene expression and function (He et al., 2021). However, despite these insights, the contribution of 5hmC to brain tissue identity has been largely unexplored. By analyzing DNA methylation patterns across multiple brain regions in 12 rhesus monkeys, our study demonstrated the requirement of 5hmC and 5mC for the establishment and preservation of tissue identity. These findings highlight the importance of considering both forms of DNA methylation in understanding the epigenetic mechanisms governing tissue-specific identity in the brain.

The observed cerebellum-specific 5mC and 5hmC patterns were closely linked to unique gene expression profiles. Genes harboring cerebellum-specific DMRs were predominantly associated with metabolic pathways, consistent with previous findings showing enrichment in signaling transduction and metabolic pathways in cerebellum-specific DhMRs (Xu et al., 2022). The convergence of these 5mC and 5hmC enrichment analyses highlights the complementary roles of both forms of DNA methylation in orchestrating cerebellar function and development, reinforcing the notion that 5mC and 5hmC collectively regulate key physiological activities essential for cerebellar identity and specialization.

The cerebellum has long been recognized for its role in coordinating motor activity, balance, and learning (D'Angelo et al., 2011). However, increasing evidence suggests it also plays a crucial role in cognitive function (Jacobi et al., 2021; Schmahmann, 2019). Impaired cerebellar function has been observed in patients with schizophrenia, manifesting as cognitive and attentional deficits, abnormalities in cerebellar morphology, and altered blood flow regulation (Andreasen & Pierson, 2008; Picard et al., 2008). These observations highlight the growing recognition of the cerebellum as a key cognitive organ with potential involvement in neurological and psychiatric disorders such as schizophrenia.

While schizophrenia is known to involve multiple brain areas, including the lateral ventricles, prefrontal lobes, and occipital lobes, as well as gray and white matter abnormalities (Cole & Cooper, 1975), the role of epigenetic regulation within the cerebellum remains largely unexplored. Notably, 5hmC has been shown to be more abundant in neurons than in other cell types (Kriaucionis & Heintz, 2009; Szulwach et al., 2011). Our findings revealed that 5hmC levels were particularly enriched in the cerebellum, suggesting that granule neurons—the predominant excitatory neuronal population in this region—may harbor higher levels of 5hmC compared to other neuronal subtypes in different brain regions. Granule cells serve as the primary conduit of information throughout the cerebellar cortex, integrating inputs from the central nervous system (Aldinger et al., 2021; Luo et al., 2022). This raises the possibility that elevated 5hmC levels in granule neurons may play a role in regulating the cerebellum-cerebral network, with potential implications for schizophrenia-related cognitive dysfunction.

This study revealed for the first time that cerebellum-specific 5hmC was highly correlated with the expression of schizophrenia-related genes. Among these, GRIN2C and PAX6 were validated as hypomethylated DhMRs, exhibiting higher transcript and protein levels in the cerebellum of rhesus monkeys. Previous preclinical studies have established that neurodevelopmental, genetic, and environmental factors contribute to glutamate dysfunction, a key pathological feature of schizophrenia (Egerton et al., 2020; Egerton & Stone, 2012). GRIN2C encodes a subtype of ionotropic glutamate receptor, and frameshift or nonsense mutations in this gene have been linked to increased risk of schizophrenia (Tarabeux et al., 2011). Similarly, PAX6 functions as a transcription factor critical for early neural differentiation and adult neuroplasticity by modulating promoter activity (Marsich et al., 2003; Tutukova et al., 2021). Mutations in PAX6 cause abnormalities in the regulation of early developmental transcription by reducing promoter activity, a process implicated in the pathogenesis of schizophrenia (Stöber et al., 1999).

Collectively, these findings suggest that cerebellum-specific 5hmC patterns may be linked to gene expression changes relevant to schizophrenia. The observed relationship between 5hmC enrichment and increased expression of schizophrenia-related genes, such as GRIN2C and PAX6, provide a potential mechanistic link between cerebellar function and schizophrenia pathophysiology. This underscores the importance of tissue-specific epigenetic modifications, such as 5hmC, in both physiological processes and pathological conditions, such as those seen in schizophrenia, within specific tissues like the cerebellum.

DNA methylation is a crucial mechanism in gene expression regulation (Semick et al., 2019), with 5mC in promoter regions typically associated with transcriptional repression, while 5hmC at enhancer regions is linked to target gene activation (Breiling & Lyko, 2015; Guibert & Weber, 2013; Johnson et al., 2016; Stroud et al., 2011). Although these regulatory functions are well established, the coordinated interplay between 5mC and 5hmC in transcriptional control remains incompletely understood. Our findings showed that 5mC in promoter regions was strongly correlated with cerebellum-specific gene expression patterns, reinforcing its crucial role in transcriptional regulation. However, the interaction between 5mC and 5hmC across both promoter and gene body regions indicated a more complex regulatory network than previously appreciated. Using the ME-class computational tool, Schlosberg et al. (2019) reported that the cumulative levels of 5mC+5hmC within 2 kb of the TSS are linked to distinct gene expression patterns. Our current findings align with and support the conclusions of Schlosberg et al. (2019), identifying two distinct modes of coordinated regulation. In one mode, 5mC increased while 5hmC decreased, correlating with lower expression of cerebellum-specific genes, while in the other mode, both 5mC and 5hmC decreased, which was associated with higher gene expression. These results suggest that 5mC and 5hmC act in a coordinated manner rather than independently, at least in certain genomic contexts. These findings also suggest that this epigenetic interplay may contribute to the regulation of disease-related genes, offering new insights into the mechanisms that shape transcriptional dynamics in the cerebellum.