
==== Front
Epigenetics
Epigenetics
Epigenetics
1559-2294
1559-2308
Taylor & Francis

39306700
10.1080/15592294.2024.2393945
2393945
Version of Record
Research Article
Research Article
Implementation of the Methyl-Seq platform to identify tissue- and sex-specific DNA methylation differences in the rat epigenome
O. H. COX ET AL.
EPIGENETICS
Cox Olivia H. a
Seifuddin Fayaz a
Guo Jeffrey a
Pirooznia Mehdi a
Boersma Gretha J. b
Wang Josh c
Tamashiro Kellie L.K. a
https://orcid.org/0000-0002-3454-8774
Lee Richard S. a
a Mood Disorders Center, Department of Psychiatry and Behavioral Sciences, The Johns Hopkins University School of Medicine , Baltimore, USA
b GGZ Drenthe Mental Health Institute, Department of Forensic Psychiatry , Assen, The Netherlands
c Agilent Technologies, Inc ., Santa Clara, USA
CONTACT Richard S. Lee rlee8@jhmi.edu Department of Psychiatry and Behavioral Sciences, Johns Hopkins University School of Medicine, 720 Rutland Ave, Ross 1068, Baltimore, MD 21205, USA
22 9 2024
2024
22 9 2024
19 1 2393945Integra21 9 2024
Integra21 9 2024
26 2 2024
23 7 2024
13 8 2024
© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.
2024
The Author(s)
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

ABSTRACT

Epigenomic annotations for the rat lag far behind those of human and mouse, despite the rat’s immense utility in pharmacological and behavioral studies and the need to understand their epigenetic mechanisms. We have designed a targeted-enrichment method followed by next-generation sequencing (Methyl-Seq) to identify DNA methylation (DNAm) signatures across the rat genome. The design reflected an attempt to create a more comprehensive investigation of the rat epigenome, as it included promoters, CpG islands, and island shores of all RefSeq genes. In this study, we implemented the rat Methyl-Seq platform and tested its ability to distinguish differentially methylated regions (DMRs) among three different tissue types, three distinct brain regions, and, in the hippocampus, between males and females. These comparisons yielded DNAm differences of differing magnitudes, many of which were independently validated by bisulfite pyrosequencing, including autosomal regions that were predicted to show the least degree of difference in DNAm between males and females. Quantitative reverse transcription PCR revealed that most genes associated with the DMRs showed tissue-, brain region-, and sex-specific differences in expression. In particular, we found evidence for sex-specific DNAm and expression differences at Tubb6, Lrrn2, Tex26, and Sox5l1, all of which play important roles in neurodevelopment and have been implicated in studies examining sex differences. Our results demonstrate the utility of the rat Methyl-Seq platform and suggest the presence of DNAm differences between the male and female hippocampus. The rat Methyl-Seq has the potential to provide epigenomic insights into pharmacological and behavioral studies performed in the rat.

KEYWORDS

Rat
Epigenetics
DNA methylation
sex differences
Methyl-seq
Charles T. Bauer Foundation (RSL) This study was funded by the James Wah Mood Disorders Scholar Fund via the Charles T. Bauer Foundation (RSL).
==== Body
pmcIntroduction

Of the three major species of mammals used for biomedical research, namely human, mouse, and rat, the genome for the rat is the most poorly annotated. The human genome has had the benefit of being the most clinically relevant, and hence was sequenced by the Human Genome Project in 2003, followed shortly by the ENCODE (Encyclopedia of DNA Elements) project that identified functional elements in the human genome [1–3]. The mouse genome was sequenced one year prior to the human genome by the International Mouse Genome Sequencing Consortium, which serves as a testament to the usefulness of Mus musculus as a mammalian model of human diseases [4]. Functional elements within the mouse genome were also published soon after its genome [5], with both human and murine ENCODE annotations continuing to grow with the incorporation of additional sets of functional elements over the years. As the third mammal to be sequenced, Rattus norvegicus or the Norway rat was sequenced in 2004 by the Rat Genome Sequencing Consortium [6]. However, unlike the first two genomes, an ENCODE counterpart for the rat does not exist. In fact, all iterations of its genome assembly available on the UCSC genome browser, from rn3 assembly in 2003 to rn7 in 2020, do not have incorporated functional or experimental data for gene regulation and expression.

One of the main reasons for the paucity of annotations of functional or regulatory elements for the rat is that the rat has been studied primarily for pharmacology, metabolism, neuroendocrinology, and behavior, and these endeavors did not necessarily require investigations into rat genetics. That distinction belonged to mice since their genome was much more amenable to manipulation by introducing knockouts, gene replacement, and transgenes to help elucidate causal roles of candidate genes in human diseases. However, recent advances in rat genetics have introduced knockout and transgenic models, and in-depth rat pharmacological studies, for instance, have frequently required investigations into the rat genome [7–9]. With increases in the capacity and reductions in the cost of next generation sequencing and the increasingly multidisciplinary nature of animal model research, it has become imperative to couple rat studies with those investigating functional elements within its genome.

To our knowledge, there have been only a few top-down designs that have leveraged the rat genome despite the existence of many iterations of its genome assembly. One of the first to do so was the rat version of the powerful CHARM array that showcased 2.1 M features across ~ 44K regions with genome-weighted smoothing of consecutive CpGs to identify differentially methylated regions (DMRs) between groups of samples [10]. We used this platform to identify hundreds of tissue- and brain region specific- DMRs across the rat genome [11]. Unfortunately, the rat CHARM is no longer commercially available. A more recent design involves an array platform that interrogates 36K CpGs that are well conserved across human, mouse, and rat [12], as well as the BeadChip mouse array that has been tested to determine methylation levels in the rat for about 5% of the probes (14K out of 300K) [13].

Unlike these platforms whose designs have incorporated the rat genome, most studies have utilized enrichment methods that are agnostic about the location of DNA methylation. These include but are not limited to MeDIP (methylated DNA Immunoprecipitation), hMeDIP (hydroxymethylated DNA Immunoprecipitation), and RRBS (reduced representation bisulfite sequencing). These techniques provide the means to enrich or reduce the amount of sequenced DNA by using antibodies against methylated or hydroxymethylated DNA (MeDIP and hMeDIP) or restriction enzymes that digest DNA at or near CpG dinucleotides, e.g., MspI. Sequencing libraries are constructed, and reads are aligned against the rat genome to provide the genomic context. These methods have been used to elucidate various conditions and genomic organizations in the rat genome, including a high resolution methylomic map of coding regions [14], hydroxymethylated genes that are associated with drug abuse [8], and phenotypic aging [15]. Despite the utility of these methods, they have certain limitations, which include a bias toward GC rich sequences and the inability to examine specific genes of interest.

Our current study reflects an attempt to implement a sequencing-based capture platform that uses the rat genome to obtain methylation information from all RefSeq genes without having to perform whole genome bisulfite sequencing (WGBS). We compare several peripheral and brain tissues, different regions of the brain, and brains between males and females to validate this platform.

Materials and methods

Animals and tissue collection

Male and female Sprague Dawley rats (10 weeks old, N = 8 per sex, Charles River Laboratories) were housed in polycarbonate rat cages in a temperature-controlled suite on a 12-h/12-h light-dark cycle with light onset at 0600 h. All animals were provided with ad libitum access to water and standard rodent chow (Envigo Harlan Teklad 2018) for 1 week after arrival. Rats were euthanized by decapitation, and two peripheral tissues (liver and blood) and three brain regions (cortex, hippocampus, and hypothalamus) were dissected. We selected these tissues because most of them except for spleen were investigated on the previously published rat CHARM platform, and we anticipated similar results using rat Methyl-Seq [11]. Importantly, we used blood and liver for peripheral tissues and cortex, hippocampus, and hypothalamus for brain regions since the rat has been used for many pharmacological and behavioral studies. Except for the hippocampus, all tissues used for Methyl-Seq were from male rats to minimize confounding factors that may arise from X chromosome methylation differences between sexes. For the brain, an adult rat brain matrix was used to obtain coronal slices that contained the brain regions of interest. The cortex was excised from a 2.5 mm thick coronal section corresponding to Bregma 5.2 to 2.7 mm. A 2.7 mm thick slice corresponding to Bregma −1.8 to −4.5 mm contained the hippocampus and the hypothalamus. For blood, 3 mL of the ACK Lysing Buffer (Quality Biological, Gaithersburg, MD) was added to 1 mL of trunk blood and incubated for 5 min to lyse the red blood cells prior to centrifugation. Dissected tissues and pelleted blood cells were stored at −80°C until processing for Methyl-Seq analysis. A second cohort of 10-week old rats (N = 8) was similarly euthanized, and tissues were collected for validation by bisulfite pyrosequencing and realtime quantitative reverse transcription PCR (RT-qPCR) of genes associated with DMRs. A third cohort of male and female rats (N = 4 per sex) were used for validation of additional DMRs and determination of levels of DNA methyltransferases 1, 3A, and 3B (Dnmt1, Dnmt3a, and Dnmt3b). All procedures were approved by the Institutional Animal Care and Use Committee at the Johns Hopkins University School of Medicine and were performed in accordance with guidelines established in the National Research Council’s Guide for the Care and Use of Laboratory Animals.

Nucleic acid extraction

All frozen tissues were stored at − 80°C. Genomic DNA (gDNA) and total RNA were extracted from tissues using the DNeasy Blood and Tissue Kit and RNeasy Mini Kit, respectively, according to the manufacturer’s instructions (Qiagen, Germantown, MD). Genomic DNA was quantified using a Qubit Fluorometer (Thermo Fisher Scientific, Waltham, MA), and RNA was quantified by 2200 TapeStation (Agilent Technologies Inc., Santa Clara, CA) according to the manufacturers’ instructions. RNA integrity number (RIN) was greater than 8.5 for all brain samples and greater than 7 for blood and liver tissues.

SureSelect Rat Methyl-Seq platform design

Several rat genome assemblies available on the UCSC Genome Browser do not contain annotations of experimental data of regulatory regions. Therefore, we focused on including promoters (±1 Kb of transcription start site) of RefSeq genes, CpG islands, island shores (±1 Kb flanking CpG islands), and > 30,000 GC-rich regions from a previously designed rat CHARM platform for DNA methylation [10,11]. Even after expunging redundant sequences, overall targeted regions exceeded 300 Mb in size, which was deemed too large compared to the mouse (109 Mb) and human (84 Mb) platforms (Agilent Technologies, Inc., Santa Clara, CA) [16]. To further reduce the target size, we sampled alternating 500 bp followed by 1 Kb of omitted sequence across targeted regions that were greater than 5 Kb. The final target size consisted of the following: 111.0 Mb size, 228.8 K unique loci, ~2.3 million on-target CpGs, and an average region size of 594 bp (Table 1). The rat Methyl-Seq kit consists of consecutive, overlapping RNA ‘baits’ or probes that target complementary DNA sequences, and it is commercially available as part of the SureSelect Target Enrichment System (Agilent Technologies Inc.).Table 1. Genomic regions targeted by the rat methyl-seq platform.

Source	Site Classification	Number of Targets	Total Bases Covered (bp)	
rat.Nov.2004.rn4.refGene	Promoters (-/+1kb)	17,620	17,620,000	
rat.Nov.2004.m4.cpgIsland
GgfAndyMasked	CpG Islands	88,554	42,699,212	
rat.Nov.2004.m4.cpgIslandGgf
AndyMasked	Shores
(CpG Islands -/+1kb)	88,554	177,108,000	
Rat CHARM	GC-rich sequences	32,889	79,820,886	
PubMed	GR-Targets	1,183	274,007	
Total Target Area (redundant)	 	228,800	317,522,105	
Final Unique Target Area	 	228,800	111,038,241	

Construction of the Rat Methyl-Seq sequencing library

A detailed protocol and demonstration of library preparation for the Methyl-Seq platform have been published [17]. Briefly, 2 μgs of gDNA from each tissue were sheared using a Covaris sonicator (Covaris, Woburn, MA) to yield 170 ~ 230bp fragments. These sheared fragments were end-repaired, 3’-adenylated, and further ligated with methylated primers. Primer-ligated DNA fragments were then denatured and hybridized to biotinylated, plus-strand DNA-complementary RNA library baits and precipitated from the solution using streptavidin-coated magnetic beads. Following RNase-digestion of the ‘baits,’ captured DNA was bisulfite-converted using the EZ DNA Methylation Gold Kit (Zymo Research, Irvine, CA). Bisulfite-converted DNA samples were PCR-amplified using sample-specific indexed (‘barcoding’) primers to allow for multiplexing and sequenced on Illumina Hi-Seq3000. Four samples were loaded on each lane.

Analysis of Sequencing Data

FASTQC version 0.11.3 was used for quality control of all the paired-end reads to assess per sequence base quality, per tile sequence quality, per sequence quality scores, per base sequence content, per sequence GC content, per base N content, sequence length distribution, sequence duplication levels, overrepresented sequences, and adapter and kmer content. Reads were trimmed using Trim Galore v0.3.7. Default parameters were used, and one or more base pairs was trimmed off at the end of all paired-end reads to improve paired-end mapping. Reads were mapped to the 2004 Rattus norvegicus assembly (Baylor 3.4/rn4), which was produced by the Rat Genome Reference Consortium, using Bismark version 0.13.0 with Bowtie 2 version 2.1.0 [18,19]. Briefly, for alignment purposes, Bismark converts all C’s to T’s (in forward reads) and all G’s to A’s (in reverse reads) prior to mapping and maps these in silico converted reads to both a C-to-T and G-to-A in silico-converted genome. After successful alignment it replaces the T’s and A’s back to their original bases in all converted reads and compares it to the original reference genome to deduce methylated cytosines. Default parameters were used with the exception that ‘bowtie2 and 1 mismatch’ was allowed during the alignment. After running Bismark, PCR duplicates were removed from the mapped reads using the ‘deduplicate bismark’ routine. Post-alignment quality control was performed using SAMtools version 0.1.19 and BamUtil version 1.0.12 [20]. The Default Bismark methylation extractor routine was used with the exception of –paired-end, –no-overlap, and minimum coverage of at least 1 read to extract all CpGs in individual samples. BSseq was used to analyze the CpG level data across each sample [21]. The package inputs the following data from alignment: genomic coordinates, M matrix or % methylation of each CpG dinucleotide, and Cov matrix or the total number of reads covering each CpG. The package then performs the following analyses: use BSmooth function to smooth raw methylation estimates (M-values) across CpGs, compute t-statistics between groups of samples using the BSmooth.tstat function, and establish threshold levels of the t-statistics to identify Differentially Methylated Regions (DMRs) using the function dmrFinder. For smoothing, default parameters were used with the exception of the smoothing window size set to 500, and the minimum number of CpGs within the smoothing window set to 20. For computing the t-statistics between two groups of tissues, CpGs with at least 10× coverage across all samples were used. The t-statistics was not smoothed or corrected. To determine the threshold of the t-statistics to use in the dmrFinder, quantiles for the t-statistics for the entire genome was used. For annotating DMRs, regions with at least three consecutive CpGs showing at least 10% mean methylation difference between two tissue groups were included. Post alignment QC statistics are shown in Supplementary Table S1. M-values for each sample were imported into the dist function, followed by hclust in R using default parameters to compute hierarchical clustering.

Bisulfite PCR and pyrosequencing

Select DMRs identified from bioinformatic analysis of different tissue, brain region, and male vs. female comparisons were further replicated and validated by bisulfite pyrosequencing, which measures methylation variation at > 90% precision [22]. The ability of the pyrosequencer to discriminate among DNA with different methylation levels were tested by treating genomic DNA from different tissues (blood, cortex, liver, hippocampus, and hypothalamus) twice with Phi29 polymerase or bacterial SssI methylase (New England Biolabs, Ipswich, MA) to generate hypo- or hypermethylated template DNA, respectively. DNA concentrations of hypo- and hypermethylated DNA were measured and appropriately mixed with each other to produce 0%, 20%, 40%, 60%, 80%, and 100% methylated template DNA. Two pyrosequencing primers each from two randomly selected DMRs (E2f7 and Zic1) and targeting 2 CpGs were used to test a total of 4 CpGs. Results demonstrate the ability of the pyrosequencer to distinguish differentially methylated regions with good accuracy (R2 >0.99, Supplementary Figure S1a-d). Design and implementation of pyrosequencing assays for rat have been described previously [11], and primers used for each genomic region are included in Supplementary Table S2.

Gene ontology and KEGG pathway analysis

Genes associated with DMRs between the male and female hippocampus were analyzed by DAVID to identify processes and pathways that may differ between the sexes [23]. To calculate enrichment, we included only genes included in the rat Methyl-Seq library.

Quantitative reverse transcription PCR (RT-qPCR)

To determine whether any of the DMRs are functionally significant, we tested transcript levels of genes nearest to each candidate DMR. For each RNA sample, the QuantiTect Reverse Transcription Kit (Qiagen) was used to generate cDNAs. Negative reverse transcriptase samples were used to ensure the absence of contaminating genomic DNA. All reactions were carried out in triplicate using 1× SYBR Green Master Mix (Thermo Fisher Scientific), forward and reverse primers for each transcript, and 30 ng of cDNA template in a total volume of 20 μl. In addition to DMR-associated transcripts, expression levels of the housekeeping gene Actb (β-actin) were used for normalization. RT-qPCR was performed on an Applied Biosystems QuantStudio 5 Real-time PCR System (Thermo Fisher Scientific) under standard PCR conditions (50°C for 2 min; 95°C for 10 min; and 60°C for 1 min for 40 cycles). Each sample was tested in triplicate, and each replicate was checked to ensure that the threshold cycle (Ct) values were all within 0.25 Ct of the other two replicates. For the determination of relative expression values, the −ΔΔCt method was used, where triplicate Ct values for each rat sample were averaged and subtracted from those derived from Actb [24]. The Ct difference for a calibrator sample was subtracted from those of the test samples, and the resulting −ΔΔCt values were raised to a power of 2 to determine normalized relative expression. Most of the primers were designed using NCBI Primer-BLAST [25], and primer pairs that produced more than one amplicon peak during melting curve analysis were redesigned. Primers already optimized for RT-qPCR for the rat DNA methyltransferases Dnmt1, Dnmt3a, and Dnmt3b were purchased from IDTDNA (Coralville, IA). A list of primers is included in Supplementary Table S2.

Results

Design of and sequencing on the rat methyl-seq platform

We optimized our rat Methyl-Seq design to include as much of the relatively unannotated rat genome as possible, while at the same time limiting the size to be similar to that of the mouse Methyl-Seq (110 Mb). In the end, the design included promoters, CpG islands, island shores, and GC-rich sequences to comprise 2.3 M on-target CpGs, 229K unique targets, and 111 Mb capture size. The breakdown of target numbers and bases are shown in Table 1.

To test the ability of the platform to detect differential DNAm levels among peripheral and brain tissues, we constructed sequencing libraries from the following16 samples: 2 liver, 2 blood, 2 cortex, 2 hypothalamus, and 8 hippocampus (4 male and 4 female). The number of samples sequenced were informed by our previous work on a rat array platform where we reported large DNAm differences among peripheral tissues, more subtle DNAm differences among brain regions, and unverifiable brain DNAm differences between the sexes [11]. To increase the power to detect smaller effect sizes, we examined 4 male and 4 female hippocampal samples. Hierarchical clustering based on raw methylation values of each sample shows similarity among the brain tissues, especially between cortex and hypothalamus, as well as inaccurate discrimination between male and female hippocampal samples (Figure 1). Figure 1. Hierarchical clustering was performed on raw DNA methylation values obtained from sequencing 16 tissues: 2 blood, 2 liver, 2 hypothalamus, 2 cortex, and 8 hippocampus samples. The hippocampus samples are further divided into 4 females and 4 males.

DNA methylation differences among peripheral and brain tissues

To test whether the rat Methyl-Seq can accurately identify differentially methylated regions (DMRs) among different types of tissues, we compared the blood, liver, and hippocampus (Supplementary Table S3). First, a comparison of DNAm patterns between the blood and hippocampus yielded 2,628 DMRs that collectively across all DMRs showed significantly higher average DNAm percentage in blood compared to hippocampus (71.7% vs. 30.7%, p = 1.7×10−455). Individually, 90.0% (2,364) of the DMRs showed higher DNA methylation in blood than hippocampus. On average, the DNAm percentage difference at each DMR was 53.9%. Similarly, a comparison of blood and liver yielded 1,722 DMRs, with the mean DNAm percentage at 76.8% for blood and 39.4% for liver (p = 1.5×10−414), with average DNAm percentage differences at each DMR of 47.3%. Also, 91.1% of DMRs (1,569) showed higher DNAm in blood compared to liver. In contrast, there were similar overall methylation levels between hippocampus and liver (55.5% vs. 55.0%, p = 0.43), with 2,463 (51.7%) of the 4,763 DMRs showing higher methylation levels in hippocampus compared to liver. Despite the similar overall DNAm levels, the average DNAm percentage difference at each DMR was 44.3%. The breakdown of the DNAm percentages is shown in Table 2.Table 2. Tissue comparisons and DMR characteristics.

Tissue Comparisons	Avg. % DNA Methylation	P-value	# of DMRs	DNAm/DMR Comparison*	
Blood vs. Hippocampus	Blood: 71.7%	Hippo: 30.7%	1.7 × 10−455	2,628	2,364 DMRs Blood > Hippo	
Blood vs. Liver	Blood: 76.8%	Liver : 39.4%	1.5 × 10−414	1,722	1,569 DMRs Blood > Liver	
Hippocampus vs. Liver	Hippo: 55.5%	Liver: 55.0%	0.43	4,763	2,463 DMRs Hippo > Liver	
Cortex vs. Hippocampus	Cortex: 65.3%	Hippo: 48.3%	1.9 × 10−283	3,619	3,337 DMRs Cortex > Hippo	
Cortex vs. Hypothalamus	Cortex: 54.2%	Hypo: 78.5%	9.8 × 10−184	1,128	29 DMRs Cortex > Hypo	
Hippocampus vs. Hypothalamus	Hippo: 43.2%	Hypo: 68.1%	2.8 × 10−276	2,098	147 DMRs Hippo > Hypo	
Hippocampus
X chromosome
Female vs. Male	Female: 38.6%	Male: 13.8%	7.8 × 10−115	672	642 DMRs Female > Male	
Hippocampus Autosomes
Female vs. Male	Female: 68.0%	Male: 57.3%	2.8×10−25	673	599 DMRs Female > Male	
*DNAm/DMR column shows the number of DMRs whose DNAm levels across each DMR is higher in one tissue than those of the other tissue.

For each tissue comparison, two DMRs with lower DNAm in one tissue and vice versa were selected for validation by bisulfite pyrosequencing. For the comparison between blood and hippocampus, we selected Ephrin Type-B Receptor 1 (Ephb1) and GRB2 Related Adaptor Protein 2 (Grap2). EPHB1 regulates cell proliferation and migration in the hippocampus [26], and its associated DMR showed lower average methylation across 10 CpGs in the hippocampus compared to blood (42.9% vs. 95.4, Figure 2a). In contrast, Grap2, which is involved in leukocyte-specific protein tyrosine kinase signaling [27], showed lower methylation levels in blood compared to hippocampus (3.3% vs. 91.6%, Figure 2d). Other comparisons showed expected results, with CAMPATH-1 Antigen (Cd52), a gene expressed in PBMCs [28], showing lower mean DNAm percentage in blood compared to liver (8.9% vs. 56.5%, Figure 2b), and Cyp17a1, which encodes the cytochrome P450 17A1 enzyme expressed in the liver [29], showing lower mean DNAm percentage in liver compared to blood (32.8% vs. 87.2%, Figure 2e). A gene that encodes another cytochrome P450 enzyme, Cyp8b1, that is involved in bile synthesis in the liver [30], showed lower methylation in liver compared to hippocampus (21.6% vs. 65.4%, Figure 2c). In the hippocampus, SH3 And Multiple Ankyrin Repeat domains 2 (Shank2), which exhibits brain-specific expression and is associated with autism [31], showed lower mean DNAm than in liver (5.0% vs. 64.7%, Figure 2f). Another cohort of animals was used to test DMRs showing the highest and lowest DNAm levels for each tissue comparison (Supplementary Figure S2a-f). These validation results demonstrate the ability of the Methyl-Seq platform to distinguish between tissues where differences are greater than 40%. Figure 2. Candidate DMRs among blood (N = 4), liver (N = 4), and hippocampus (N = 4, males only) were independently validated by bisulfite pyrosequencing. DMR validation results are displayed by the two tissues being compared followed by the name of the gene closest to each DMR. Results for (a) Ephb1, (b) Cd52, (c) Cyp8b1, (d) Grap2, (e) Cyp17a1, and (f) Shank2 are shown. Bar graphs are represented as mean ± SEM. ***p<.001.

DNA methylation differences among three brain regions

Once we confirmed the platform’s ability to distinguish among different types of tissues, we asked whether it could also distinguish DNAm patterns in the same tissue type, i.e., the brain, but in different regions (cortex, hippocampus, and hypothalamus) that have distinct functions (Supplementary Table S4). For DMRs between cortex and hippocampus, average DNAm levels for both tissues were 65.3% and 48.3%, respectively, with 3,337 out of 3,619 DMRs (92.2%) showing higher methylation in cortex (p = 1.9×10−283, Table 2). The average difference in DNAm percentage at each DMR between the two tissues was 20.3%. A comparison of DNAm between cortex and hypothalamus showed lower % methylation in cortex compared to hypothalamus (54.2% vs. 78.5%, p = 9.8×10−184) with only 29 out of 1,128 DMRs (2.6%) showing higher DNAm levels in cortex. The average difference in DNAm percentage at each DMR between the cortex and hypothalamus was 25.4%. We then characterized DMRs between hippocampus and hypothalamus. Hippocampal tissues showed an average of 43.2% DNAm compared to 68.1% in hypothalamus. Similar to the cortex vs. hypothalamus comparison, there were only 147 DMRs out of 2,098 (7.0%) where hippocampus showed DMRs with higher methylation than hypothalamus. On average, difference in DNAm at each DMR between the two brain regions were 27.8%.

We then selected two candidate DMRs from each comparison for further validation by pyrosequencing. For the cortex vs. hippocampus comparison, we tested Dlgap2, a cortex-expressed gene that encodes the SAP90/PSD-95-associated protein 2 [32], which showed lower methylation in cortex compared to hippocampus (60.4% vs. 89.6%, Figure 3a). In contrast, Smad5 (Mothers Against Decapentaplegic Homolog 5), a member of receptor proteins that binds bone morphogenetic proteins in the brain [33], showed relatively higher DNAm in the cortex than in the hippocampus (60.1% vs. 36.9%, Figure 3d). For the cortex vs. hypothalamus comparison, we tested Tcf4 (Transcription Factor 4), which has been implicated in schizophrenia in the context of cortical development [34], and Cadps (Calcium Dependent Secretion Activator), which showed differential expression in the cortex vs. hypothalamus in the canine brain [35]. The Tcf4-associated DMR in the cortex had methylation levels that were approximately half that in the hypothalamus (45.4% vs. 94.3%, Figure 3b), whereas the Cadps-associated DMR showed similar methylation levels between both brain regions (86.4% for cortex and 80.2% for hypothalamus, Figure 3e). We next examined DNAm differences for two candidate DMRs associated with Cit and Mycn. The Cit gene encodes the Citron Kinase whose activity is required for postnatal neurogenesis in the hippocampus and cortex [36,37], and DNAm levels were higher in the hippocampus compared to the hypothalamus (77.9% vs. 61.1%, Figure 3c). The Mycn gene encodes a member of the Myc family of oncogenes that is overexpressed in neuroblastomas but is also crucial for neurogenesis and oligodendrogenesis [38]. Lower DNA methylation levels were observed in the hippocampus compared to the hypothalamus (51.0% vs. 77.1%, Figure 3f). The rat Methyl-Seq predicted DMRs among the three brain region comparisons whose %DNAm differences were on average smaller in magnitude (less than 30%) compared to those among different tissues. However, there were many DMRs whose DNAm differences exceeded 50%. DMRs showing the greatest DNAm difference in either direction, e.g., the lowest and highest DNAm levels in the hippocampus when compared to cortex, are shown in Supplementary Figure S3a-f. Figure 3. Candidate DMRs among cortex (N = 4), hippocampus (N = 4, males only), and hypothalamus (N = 4) were independently validated by bisulfite pyrosequencing. DMR validation results are displayed by the two brain regions being compared followed by the name of the gene closest to each DMR. Results for (a) Dlgap2, (b) Tcf4, (c) Cit, (d) Smad5, (e) Cadps, and (f) Mycn are shown. Bar graphs are represented as mean ± SEM. **p < .01 and ***p < .001.

DNA methylation differences between males and females

The dendrogram in Figure 1 was not successful in fully distinguishing the DNAm levels between the male and female hippocampus and suggested that the overall DNAm levels across the genome may be similar. Despite this possibility, we examined DNAm differences between male and female rats with the DMR list segregated by autosomes and X chromosomes (Supplementary Table S5). For X chromosome DMRs, female hippocampus samples showed 38.6% DNAm on average compared to 13.8% in males (p = 7.8 × 10−115). The average difference in DNAm percentage at each DMR between the sexes was 26.5%, with 642 of the 672 DMRs (95.5%) having higher DNAm levels in females than males as expected, due to the methylated, inactive second X chromosome in females (Table 2). We then explored DNAm differences at autosomal regions between the sexes. DNAm levels of DMRs at autosomes were 68.0% for females and 57.3% in males, with an overall difference of only 10.7% (p = 2.8 × 10−25). The average difference in DNAm percentage at each autosomal DMR between the sexes was 13.7%. Interestingly, 599 of the 673 autosomal DMRs (89.0%) showed higher DNAm in females.

Sex-specific DMRs were also subjected to validation by pyrosequencing. For X chromosomal DMRs, we tested two out of five DMRs associated with the Xist/Tsix locus (LOC680227 in Supplementary Table S5). The locus consisted of the Xist (X inactive specific transcript), a 17-kb long non-coding RNA that plays a pivotal role in the silencing of the inactive X chromosome in females [39], and its anti-sense transcript Tsix responsible for the prevention of Xist accumulation on the active X chromosome in both males and females [40]. Twelve CpGs in the promoter of Xist showed mean DNAm of 49.5% in females compared to 90.1% in males (Figure 4a). In contrast, six CpGs in the promoter CpG island of Tsix showed mean DNAm of 29.1% in females compared to 3.8% in males (Figure 4b). Other DMRs showed expected DNAm patterns consistent with the presence of substantial methylation only in females, presumably on the inactive X chromosome. For instance, at Pgk1, which encodes the glycolytic enzyme Phosphoglycerate Kinase 1, female hippocampal samples showed a mean DNAm level of 34.6%, whereas male samples showed a much lower mean DNAm level of 2.3% (Figure 4c). Other DMRs showed methylation patterns that deviated from those expected between the sexes. The Bcor gene encodes a BCL6 corepressor protein involved in embryonic development [41]. One of several DMRs associated with Bcor showed lower mean DNAm levels in females compared to males (34.0% vs. 51.3%, Figure 4d). Similarly, a DMR associated with Dmd, the gene that encodes dystrophin whose loss of function is associated with Duchenne muscular dystrophy [42], showed only subtle mean DNAm differences between females vs. males (76.4% vs. 86.5%, Figure 4e). Figure 4. Candidate X-linked DMRs in the hippocampus between females (N = 4) and males (N = 4) were independently validated by bisulfite pyrosequencing. Results for (a) Xist, (b) Tsix, (c) Pgk1, (d) Bcor, and (e) Dmd are shown. Bar graphs are represented as mean ± SEM. *p < .05 and ***p < .001.

We then validated autosomal DMRs predicted by comparing DNAm patterns between male and female hippocampal samples. We first examined the gene that encodes Tubulin Beta 6 Class V (Tubb6), based on the ability of this gene to respond to androgen signaling [43]. Pyrosequencing analysis of its intronic DMR on chromosome 18 showed increased DNAm at seven out of nine CpGs in females (Figure 5a). Similarly, another DMR that overlaps two exon/intron boundaries in the Plekhg5 locus (Pleckstrin Homology and RhoGEF Domain Containing G5) showed increased DNA methylation across five out of ten CpGs (Figure 5b). Interestingly, the human PLEKHG5 was implicated in a population study of ALS, a disease with sex bias against males [44,45]. In contrast, two DMRs associated with Lrrn2 (Leucine Rich Repeat Neuronal 2, Figure 5c) and Nrap (Nebulin Related Anchoring Protein, Figure 5d) showed higher DNAm in males. Lrrn2 was implicated in a large transcriptome-wide association study of prostate cancer of over 140,000 individuals, whereas Nrap was identified in an Alzheimer’s Disease study as a transcript with sex-specific expression in the parahippocampal gyrus [46]. Additional DMRs include regions associated with Pou1f1 (POU Class 1 Homeobox 1), Tex26 (Testis Expressed 26), Magoh (Mago Homolog), Sox5l1 (also known as SRY-box Transcription Factor 5, Sox5), Thrb (Thyroid Hormone Receptor Beta), and Ccdc41 (Coiled-coil domain-containing protein 41, Figures 5e-j). All four CpGs in the DMR associated with Pou1f1 showed increased DNAm in males (Figure 5e). In the somatotrophs of the anterior pituitary, loss of function of the leptin receptor has been shown to reduce levels of POU1F1 proteins in females but not males [47]. Tex26 and Sox5l1 were chosen for their membership in families of genes involved in sex determination and neurodevelopment [48–50]. Tex26 was recently identified as a pituitary expressed gene associated with neuroticism in males [51]. Two of the four CpGs in the DMR associated with Tex26 showed increased DNA methylation in females, consistent with its male-specific expression (Figure 5f). In contrast, three of the five CpGs in the DMR associated with Sox5l1 showed increased DNA in males (Figure 5g). Only one of the four CpGs in the DMR associated with Thrb showed higher DNAm in females (Figure 5h). Interestingly, protein levels of THRB1/2 were higher in male compared to female placentas of mothers with gestational diabetes mellitus [52]. Other DMRs include those associated with Magoh, whose CpGs showed higher DNAm in males (Figure 5i) and those of Ccdc41 that showed higher DNAm levels in females (Figure 5j). Magoh was one of many genes that exhibit increased expression and alternative splicing in the male gonads [53], whereas Ccdc41 was highly expressed in fetal ovaries [54].

Figure 5. Candidate autosomal DMRs in the hippocampus between females (N = 4) and males (N = 4) were independently validated by bisulfite pyrosequencing. Results are displayed by the chromosomal location of each DMR along with the name of its nearby associated gene: (a) Tubb6, (b) Plekhg5, (c) Lrnn2, (d) Nrap, (e) Pou1f1, (f) Tex26, (g) Sox5l1, (h) Thrb, (i) Magoh, and (j) Ccdc41. Bar graphs are represented as mean ± SEM. *p<.05, **p<.01, and ***p<.001.

Gene ontology (GO) and pathway analysis of sex-specific DMRs

We asked whether DMRs between the male and female hippocampus could be categorized into specific cellular processes and pathways. GO analysis identified relevant processes under the Cellular Component (CC) subcategory, including pre- and postsynaptic membranes as well as GABAergic and glutamatergic synapses (Table 3). Under the Molecular Function subcategory, there were multiple processes involving DNA binding, transcription, and chromatin binding. In contrast, KEGG analysis identified calcium signaling as the only FDR-significant pathway (q = .04).Table 3. Gene Ontology (GO) analysis of sex-specific DMRs.

Term in Cellular Component Subcategory	Count	%	P Value	Fold
Enrichment	FDR	
GO:0005634~nucleus	300	34.2	9.4E-09	1.3	5.9E-06	
GO:0045121~membrane raft	35	4.0	4.8E-07	2.6	1.5E-04	
GO:0005737~cytoplasm	301	34.4	3.7E-06	1.2	7.8E-04	
GO:0045202~synapse	60	6.8	1.3E-05	1.8	0.002	
GO:0098982~GABA-ergic synapse	18	2.1	4.8E-05	3.2	0.006	
GO:0043025~neuronal cell body	50	5.7	6.7E-05	1.8	0.007	
GO:0014069~postsynaptic density	34	3.9	7.4E-05	2.1	0.007	
GO:0005654~nucleoplasm	170	19.4	1.6E-04	1.3	0.012	
GO:0045211~postsynaptic membrane	24	2.7	3.3E-04	2.3	0.023	
GO:0042734~presynaptic membrane	19	2.2	6.1E-04	2.5	0.038	
GO:0030425~dendrite	44	5.0	6.6E-04	1.7	0.038	
GO:0000785~chromatin	33	3.8	7.6E-04	1.9	0.040	
GO:0098978~glutamatergic synapse	46	5.3	0.001	1.7	0.050	
Term in Molecular Function Subcategory	Count	%	P Value	Fold
Enrichment	FDR	
GO:0003682~chromatin binding	56	6.4	3.0E-09	2.4	2.9E-06	
GO:0043565~sequence-specific DNA binding	45	5.1	1.3E-07	2.4	3.6E-05	
GO:0003677~DNA binding	82	9.4	1.4E-07	1.8	3.6E-05	
GO:0005515~protein binding	117	13.4	1.5E-07	1.6	3.6E-05	
GO:1990837~ss ds DNA binding	47	5.4	3.7E-06	2.1	7.4E-04	
GO:0003700~tf activity, ss DNA binding	42	4.8	1.4E-05	2.1	0.002	
GO:0042802~identical protein binding	112	12.8	1.5E-04	1.4	0.021	
GO:0000978~RNA polymerase II core promoter proximal region ss DNA binding	69	7.9	1.7E-04	1.6	0.021	
GO:0003714~transcription corepressor activity	20	2.3	2.0E-04	2.6	0.022	
GO:0003713~transcription coactivator activity	23	2.6	2.8E-04	2.4	0.025	
GO:0000977~RNA polymerase II regulatory region ss DNA binding	32	3.7	2.8E-04	2.0	0.025	
GO:0004879~RNA polymerase II transcription factor activity, ligand-activated ss DNA binding	10	1.1	3.2E-04	4.5	0.026	
GO:0001228~transcriptional activator activity, RNA polymerase II transcription regulatory region ss binding	39	4.5	4.0E-04	1.8	0.030	
GO:0061629~RNA polymerase II ss DNA binding tf binding	21	2.4	5.8E-04	2.4	0.041	
GO:0046972~histone acetyltransferase activity (H4-K16 specific)	4	0.5	6.9E-04	19.1	0.045	
GO:0046332~SMAD binding	10	1.1	7.3E-04	4.0	0.045	
ss – sequence specific, ds – double stranded, tf – transcription factor

Functional assessment of dmr-associated genes by RT-qPCR

Finally, we measured gene expression of peripheral and brain tissues to determine whether any of the validated DMRs were associated with gene function. RT-qPCR of EphB1, Cd52, Cyp17a1, Cyp8b1, and Shank2 all showed relatively low levels in tissues in which their DNAm levels were significantly higher than those compared. For instance, the DMR associated with Cyp8b1 had shown much higher DNAm levels in the hippocampus than liver (65.4% vs. 21.6%, Figure 2c). Consistently, the hippocampus showed significantly lower expression levels when compared to liver. In contrast, the DMR associated with Shank2 had shown higher DNAm levels in the liver compared to hippocampus (64.7% vs. 5.0%, Figure 2f). This pattern was associated with relatively low levels in the liver compared to hippocampus (Figure 6a). Testing DMRs among brain tissues showed smaller expression differences between two brain regions that were compared. For example, the DMR associated with Dlgap2 had shown lower methylation in the cortex compared to the hippocampus (60.4% vs. 89.6%, Figure 3a), whereas the DMR associated with Smad5 had shown higher methylation in the cortex compared to the hippocampus (60.1% vs. 36.9%, Figure 3d). Consistently, expression levels were higher in the cortex for Dlgap2 and lower for Smad5 compared to levels of same genes in the hippocampus (Figure 6b). There was no expression difference in Cadps between the cortex and hypothalamus, despite its association with a DMR that showed significantly lower DNAm in the hypothalamus (Figure 3e). Mycn showed a significantly higher expression in the hypothalamus compared to hippocampus, despite showing higher DNAm in the hypothalamus (Figure 3f). Comparison of expression levels of DMR-associated genes in three brain regions yielded three differentially expressed genes (Dlgap2, Smad5, and Tcf4) that were consistent with their DNAm differences against other brain regions. Figure 6. RT-qPCR was performed on genes associated with tissue- and brain region-specific DMRs. RT-qPCR results are displayed by the two tissues (a) or brain regions (b) being compared along with the name of the gene closest to each DMR. Bar graphs are represented as mean ± SEM. *p<.05 and ***p<.001. N = 4 samples per group. It should be noted that expression levels close to 0 does not necessarily translate to zero expression levels but rather significantly lower expression levels in one tissue compared to another.

X-linked genes were also examined. As expected, Xist and Tsix, which are involved in silencing and activation of X chromosomes, respectively, showed completely reciprocal pattern of expression between male and female hippocampal samples (Figure 7a). However, expression levels were similar for not only Pgk1 that showed the expected higher DNAm in females, but also for Bcor and Dmd that unexpectedly showed higher DNAm patterns in males. We also investigated autosomal regions that showed differential methylation and found that several showed differential expression between males and females. Significant differences in expression were observed for Tubb6 and Lrrn2, while nonsignificant differences, albeit consistent with DNAm patterns for the two respective genes, were observed for Plekhg5 and Nrap (Figure 7b). Additional genes that showed sex-specific differences in expression include Tex26, Sox5l1, and Thrb. Pou1f1, Magoh, and Ccdc41 showed no significant differences between the sexes (Figure 7c). Figure 7. RT-qPCR was performed on genes associated with sex-specific DMRs. Results are displayed by chromosomal locations of X-linked (a) or autosomal (b and c) gene expression levels being compared. Bar graphs are represented as mean ± SEM. *p <. 05, **p < .01, and ***p < .001. N = 4 animals per group. It should be noted that expression levels close to 0 does not necessarily translate to zero expression levels but rather significantly lower expression levels in one tissue compared to another.

Figure 7. (Continued).

Finally, we measured the expression of Dnmt1, Dnmt3a, and Dnmt3b to determine whether unequal number of DMRs between different tissue, brain region, and sex comparisons can be explained by differential expression of DNA methyltransferases. We observed higher levels of Dnmt1 and Dnmt3a in the hippocampus and liver compared to blood (Supplementary Figure S4a and S4b). Similar levels were observed for Dnmt1 and Dnmt3a for hippocampus and liver, while more than 4.1-fold increase in the expression level of Dnmt3b was observed in the liver compared to hippocampus (p = 0.008, Supplementary Figure S4c). The cortex showed significantly higher levels of Dnmt1 (2.7-fold, p = 0.005) and Dnmt3a (3.1-fold, p = 0.009) compared to hippocampus (Supplementary Figure S4d) and showed higher levels of Dnmt3a (1.9 fold, p = 0.04) compared to hypothalamus (Supplementary Figure S4e). Between the hypothalamus and hippocampus, we observed significantly higher levels of Dnmt1 (2.3-fold, p = 2.7×10−5) and Dnmt3a (1.7-fold, p = 0.001, Supplementary Figure S4f) in the hypothalamus. No significant differences in the expression levels of the three methyltransferases were observed between the male and female hippocampus (Supplementary Figure S4g).

Discussion

In this study, we sought to methodically test whether an enrichment-based sequencing platform against a poorly annotated genome can reliably detect methylation differences across tissue types. This is important because the rat genome lacks annotations of regulatory elements such as those from GENCODE, Ensembl, and ORegAnno as well as additional experimental annotations derived from DNAse I hypersensitive sites and cancer- and tissue-specific DMRs that are available for the mouse and human genomes. Given these significant limitations, we focused on promoters, CpG islands, and island shores. The island shores that flank CpG islands and contain less GC density than CpG islands were chosen based on findings from the human CHARM platform that implicated these regions as ones that preferentially undergo cancer-specific DNAm changes and whose associated genes are differentially expressed between tissues [10]. In addition, we added more than 30,000 GC-rich regions from the previous array-based rat CHARM platform, the only other platform designed for the rat but is no longer commercially available, that identified hundreds of tissue and brain region-specific DMRs [11].

The rat CHARM platform design was derived from its human version that featured several advantages over other methods commonly used at that time, including MeDIP and HELP (HpaII tiny fragment enrichment by ligation-mediated PCR) assays [55]. The authors reported that the sensitivities of MeDIP and HELP were heavily dependent on CpG density, which affected the ability to identify DMRs between HCT116 vs. HCT116 DNMT1/3 KO cell lines across a common testing platform. Since both methods rely on the enrichment of CpGs across the genome, they did not allow genomic target selection. Once next-generation sequencers became more widely used, many of the methods were transitioned from array platforms to sequencing. Techniques such as RRBS (reduced representation bisulfite sequencing) enabled genome-wide investigations based on the presence of MspI recognition sites [56]. However, RRBS is limited by the size selection process as well as the inability to target specific genes. At the same time, whole genome bisulfite sequencing (WGBS), while considered the gold standard, requires substantial read depth to determine accurate DNA methylation percentage, which can be costly at the whole genome level. Further, a noteworthy WGBS study reported that 70–80% of the sequencing reads across WGBS data provided little relevant information about CpG methylation [57]. The most widely used methylation platform to date is the Illumina BeadChip array. However, the recent mouse BeadChip interrogates 285K CpGs compared to the 2.3 million CpGs covered by the rat Methyl-Seq. Further, there is no BeadChip platform for the rat epigenome. Therefore, our goal of developing a gene-centric methylation platform for the rat that strikes a balance between comprehensive coverage and sensitivity is well warranted.

Using the rat Methyl-Seq platform, we first identified DNAm differences among blood, hippocampus, and liver, followed by differences among cortex, hippocampus, and hypothalamus. The magnitude of DNAm differences among different types of tissues, ranging from 44 to 54% among blood, hippocampus, and liver, was similar to that predicted by the previous rat CHARM experiment. These relatively large DNAm differences likely reflect different sets of genes that must be expressed to perform tissue-specific tasks such as mounting the immune response, regulating mood, or detoxification. At the same time, there must be a set of genes that are not specific to any tissue and are involved in housekeeping processes. The magnitude of DNAm differences in the same brain tissue but among different regions, ranging from 20 to 28%, was also similar to that predicted by the rat CHARM experiment and about half of the range of DNAm differences observed among different tissues [11]. The smaller DNAm differences among the brain regions likely stem from their common function in neurotransmission, i.e., by utilizing a similar set of genes such as ion channels and neurotransmitter transporters and receptors, but also due to the presence of glial cell populations that likely perform very similar functions in these different regions. The difference in DNAm among the cortex, hippocampus, and hypothalamus have likely arisen from their more nuanced, differentiated functions such as sensory integration, learning and memory, or homeostasis maintenance, respectively. For sex-specific DMRs, the rat CHARM array provided a list of only autosomal DMRs between males and females, as X-linked loci were precluded from the array design. However, we could not independently validate any of the autosomal DMRs, as pyrosequencing showed no differences in DNAm (data not shown). On the other hand, the rat Methyl-Seq approach was able to identify regions that could be validated by pyrosequencing.

The magnitude of X-linked DNAm differences between the sexes were on average 27%, whereas autosomal differences were 11%. The X-linked differences were substantially lower than the 46% difference at gene promoters observed between human male and female neutrophil DNA in one study [58], which may be more consistent with the average DNAm between a hypomethylated, active X chromosome and a hypermethylated, inactive X chromosome in females. However, we note that the mean DNAm levels for males and females are 39% and 14%, respectively (Table 2), which reflects CpGs with substantial DNAm in males. In fact, although the median DNAm percentage for the male X chromosome was 4%, more than 50 out of 672 DMRs have mean DNAm levels greater than 50% (Supplementary Table S5), including those at Bcor and Dmd (Figure 4d,e) that contributed to the elevated DNAm levels. It is also possible that X-linked methylation levels may deviate from the theoretical baseline of 0% DNAm for males and 50% for females across different tissues. Such deviation may be present in sexually dimorphic organs such as the brain where sex-determining and neurodevelopmental genes on the X chromosomes may play a greater role in sex bias observed in some psychiatric and neurological disorders. However, an in-depth survey of X chromosomal methylation across multiple tissues and organs are needed to answer this question.

One of the most interesting findings from this platform validation study is the ability of the rat Methyl-Seq to detect autosomal DMRs between sexes. We identified just as many DMRs on autosomes as on the X chromosome. Also, most of the autosomal DMRs (89%) showed higher DNAm in females, only slightly less than on the X chromosome (95.5%), and it is currently not clear why. It is possible that some of the epigenetic machinery involved in X inactivation in females, e.g., Xist expression and H3K27 methylation, may influence epigenetic patterns at autosomal loci that are near the inactive X chromosome. However, it is widely established that Xist RNA accumulates in cis, with some exceptions [59,60]. There was one notable difference between X-linked and autosomal DMRs. While X-linked DMRs were hypomethylated in males (14%), autosomal DMRs showed a more robust magnitude of DNAm in males (57%). The higher mean DNAm values of autosomal DMRs reflect the absence of chromosome-wide activation or silencing mechanisms as found on the X chromosome. As we have sampled only a fraction of each gene, it is highly likely that there are many more sex-specific DMRs throughout the genome that when identified can further elucidate an underlying mechanism of these sex-specific DMRs.

Another notable finding from our study is the association between DNA methylation and gene expression. DMRs observed among blood, hippocampus, and liver and further validated by pyrosequencing were associated with relatively higher gene expression from the tissue with lower DNAm. Similarly, DMRs among the cortex, hippocampus, and hypothalamus mostly showed a negative relationship between DNAm and gene expression, i.e., higher DNAm was generally associated lower gene expression. The only exception occurred at the DMR associated with Cadps, where gene expression levels were similar between the cortex and hypothalamus despite the methylation differences between these brain regions. We also note that gene expression differences were not as drastic as those among the blood, hippocampus, and liver. This observation is not surprising given the smaller magnitude of difference in DNAm among the three brain regions (20 to 28%) compared to different tissues (44 to 54%). Similar to DNA methylation differences, expression differences in the brain regions may be smaller due to the presence of glial cells in all brain regions.

On the X chromosome, we observed drastic sex-specific differences in gene expression at Xist and Tsix. These patterns are consistent with similarly drastic differences in DNAm at these loci. However, expression levels of three X-linked genes, one with expected DNAm patterns (Pgk1) and two with higher DNAm in males (Bcor and Dmd), showed no sex preference. Expression patterns for Bcor and Dmd are unexpected given the significantly higher DNAm levels in males and potential gene silencing. It is likely that these DMRs do not play direct roles in gene expression or may be involved in regulating alternative transcripts. Regardless, all three genes underwent proper dosage compensation and resulted in similar expression patterns between sexes. At autosomal DMRs, only a subset of genes, namely Tubb6, Lrrn2, Tex26, and Sox5l1, was associated with differences in expression. For Magoh, Thrb, and Ccdc41, we note that the DMRs were located more than 35,000 bps away from the gene, raising the possibility that their ability to regulate nearby genes may not be strong (Supplementary Table S5). Also, given that these DNAm differences were less than 11% on average, it is no surprise that expression differences were also subtle. It is also possible that many of these DMRs do not participate in establishing sex-specific expression of their associated genes. We note that not all DMRs are always associated with regulation of gene expression. DMRs can represent remnants from a previous developmental period when sex hormone signaling may have contributed to their formation. These may influence sex-specific expression of their associated genes following a surge in sex hormones. It is reminiscent of loss of DNAm due to glucocorticoid exposure, its persistence in the absence of glucocorticoids, and its influence on gene function only after another bout of glucocorticoid signaling [61]. Regardless, we identified autosomal regions that harbor sex-specific differences in DNA methylation and are associated with genes that show sex preference in expression. A more comprehensive transcriptomic method especially following exposure to sex hormones might better delineate the DNAm-expression relationship between males and females.

We also asked whether expression levels of DNA methyltransferases can explain the skewed number of DMRs in one tissue compared to another. Although we observed significant differences in the levels of these genes in various tissue and brain region comparisons, we did not observe a consistent expression pattern for any of the three genes that can explain the skewed number of DMRs in these comparisons. For instance, Dnmt1 and Dnmt3a showed higher levels of expression in both hippocampus and liver when compared to blood. However, results in Table 2 show higher average DNAm levels in blood as well as a significant number of blood DMRs showing higher DNAm levels compared to both hippocampus and liver. In the brain, the cortex shows higher Dnmt1 and Dnmt3a levels, higher average DNAm levels in across the DMRs, and most of the DMRs showing higher methylation compared to the hippocampus. This observation also applies to the hippocampus-hypothalamus comparison where the hypothalamus shows higher expression levels of both Dnmt1 and Dnmt3a. However, the potential role of these Dnmts in establishing brain region-specific DMR patterns is contradicted in the cortex-hypothalamus comparison where Dnmt3a is expressed at higher levels in the cortex that shows lower average DNAm levels across the DMRs and a much smaller number of DMRs that have higher DNAm than the hypothalamus. While interesting, expression levels of Dnmts in these tissues do not provide compelling evidence for their exclusive role in establishing these DMRs.

The current study was conducted to determine whether data generated using this platform can be validated by a more sensitive approach, and we did this by comparing sets of samples that we knew would show differences in methylation. Results from this study support the findings from several studies that have already explored different aspects of the rat epigenome using the current platform. An early pilot study by our group involved a protocol video that employed a rat model of stress to identify epigenetic changes associated with stress hormone exposure [17]. Others have since then used the platform to investigate the epigenetic effects of carcinogens [62], teratogens and neurotoxins [63], and X-irradiation [64]. Together, we have established the rat Methyl-Seq platform as a useful tool for studying the molecular and epigenetic basis of physiology and behavior in rats.

Despite some of the strengths highlighted above, our study has some limitations. First, we used only 2 samples each from blood, liver, cortex, and hypothalamus. The small sample size was due to our a posteriori evidence from the rat CHARM platform where two samples provided sufficient power to identify hundreds of tissue-specific and brain region-specific DMRs. Using additional samples will likely identify more DMRs. Also, the Methyl-Seq platform does not distinguish between methylcytosine (mC) and hydroxymethylcytosine (hmC), as the bisulfite conversion process does not affect either residue, and both are read as cytosines. Looking at hmC modifications would be important for comparing brain tissues since they are most abundant in the brain [65,66]. However, most of the DMRs, including several of the sex-specific DMRs, showed a consistent inverse relationship between DNA methylation and gene expression. Given the positive association of hmCs with gene expression, we would suspect and test for the presence of hmCs if any DMRs showed such a relationship [67]. Other limitations include limited read depth in the current study (mean read depth of 31.5X for all CpGs with at least 10 reads). We might have identified additional sex-specific DMRs with greater read depth. We also acknowledge the fact that the methylome of the outbred Sprague Dawley rat used in the current study differs from those of other commonly used strains of rats in terms of strain-specific DNAm and underlying genetic variations, both of which can impact tissue-, brain region-, and sex-specific methylation. Further, we have yet to explore other epigenetic mechanisms, namely histone modifications and non-coding RNAs, that can potentially confer tissue-, brain region-, and sex-specific gene expression. Despite these limitations, we believe that the current platform is a promising start to begin to unravel the poorly annotated rat genome.

In summary, adaptation of the Methyl-Seq enrichment kit for the rat genome provides a robust platform to identify DNA methylation differences that can be validated by independent, more sensitive methods, i.e., bisulfite pyrosequencing. We successfully used the Methyl-Seq platform to validate methylation differences among tissues and distinct regions of the brain. We additionally demonstrate that Methyl-Seq can detect more subtle sex-specific differences in DNAm that arise in the same tissue from male versus female animals. Deposition of additional regulatory regions, by using other approaches as mentioned above, can make the platform a more powerful tool to investigate the rat epigenome.

Supplementary Material

Supplemental Material

Supp Fig 2.jpeg

Supplementary Table 3 Tissue DMRs.xlsx

Supp Fig 3.jpeg

Supp Fig 1.jpeg

Supplementary Table 1 rat MethylSeq post QC.xlsx

Supplementary Table 4 Brain DMRs.xlsx

Supplementary Table 2 List of Primers.xlsx

Supplementary Table 5 Male_Female DMRs.xlsx

Supp Fig 4.jpeg

Acknowledgments

The authors thank Lindsey K. Macias and Jasmine Shakir for their technical assistance.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The FASTQ data that support the findings of this study are openly available at the NIH Sequence Read Archive (SRA) after 16 August 2024, at https://www.ncbi.nlm.nih.gov/sra, reference number PRJNA997727.

Abbreviations and acronyms

DNAm (DNA methylation); DMR (differentially methylated region)

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2024.2393945
==== Refs
References

[1] Green ED, Watson JD, Collins FS. Human genome project: twenty-five years of big biology. Nature. 2015;526 (7571 ):29–22. PMC5101944 doi:10.1038/526029a 26432225
[2] Collins FS, Morgan M, Patrinos A. The human genome project: lessons from large-scale biology. Science. 2003;300 (5617 ):286–290. doi: 10.1126/science.1084564 12690187
[3] Consortium EP. The ENCODE (ENCyclopedia of DNA elements) project. Science. 2004;306 :636–640.15499007
[4] Genome Sequencing CM, Waterston RH, Lindblad-Toh K, et al. Initial sequencing and comparative analysis of the mouse genome. Nature. 2002;420 :520–562.12466850
[5] Mouse EC, Stamatoyannopoulos JA, Snyder M, et al. An encyclopedia of mouse DNA elements (mouse ENCODE). Genome Biol. 2012;13 (8 ):418. PMC3491367. doi:10.1186/gb-2012-13-8-418 22889292
[6] Gibbs RA, Weinstock GM, Metzker ML, et al. Genome sequence of the Brown Norway rat yields insights into mammalian evolution. Nature. 2004;428 (6982 ):493–521. doi: 10.1038/nature02426 15057822
[7] Huang G, Ashton C, Kumbhani DS, et al. Genetic manipulations in the rat: progress and prospects. Current Opinion In Nephrology And Hypertension. 2011;20 (4 ):391–399. PMC3857098 doi: 10.1097/MNH.0b013e328347768a 21546835
[8] Cadet JL, Brannock C, Krasnova IN, et al. Genome-wide DNA hydroxymethylation identifies potassium channels in the nucleus accumbens as discriminators of methamphetamine addiction and abstinence. Mol Psychiatry. 2017;22 :1196–1204. PMC7405865.27046646
[9] Lee RS, Pirooznia M, Guintivano J, et al. Search for common targets of lithium and valproic acid identifies novel epigenetic effects of lithium on the rat leptin receptor gene. Transl Psychiatry. 2015;5 (7 ):e600. PMC5068731. doi:10.1038/tp.2015.90 26171981
[10] Irizarry RA, Ladd-Acosta C, Wen B, et al. The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific CpG island shores. Nat Genet. 2009;41 (2 ):178–186. PMC2729128. doi:10.1038/ng.298 19151715
[11] Lee RS, Tamashiro KL, Aryee MJ, et al. Adaptation of the CHARM DNA methylation platform for the rat genome reveals novel brain region-specific differences. Epigenetics. 2011;6 (11 ):1378–1390. PMC3242812. doi:10.4161/epi.6.11.18072 22048247
[12] Arneson A, Haghani A, Thompson MJ, et al. A mammalian methylation array for profiling methylation levels at conserved sequences. Nat Commun. 2022;13 (1 ):783. PMC8831611 (publication number WO2020150705) related to this work for which A.A. B.B. J.E. and S.H. are named inventors. S.H. is a founder of the non-profit Epigenetic Clock Development Foundation, which has licensed several patents from his employer UC Regents, and distributes the mammalian methylation array. Bret Barnes is an employee for Illumina Inc which manufactures the mammalian methylation array. The remaining authors declare no competing interests.doi:10.1038/s41467-022-28355-z 35145108
[13] Zhou W, Hinoue T, Barnes B, et al. DNA methylation dynamics and dysregulation delineated by high-throughput profiling in the mouse. Cell Genom. 2022;2 1 100144. doi:10.1016/j.xgen.2022.100144. PMC9306256.35873672
[14] Sati S, Tanwar VS, Kumar KA, et al. High resolution methylome map of rat indicates role of intragenic DNA methylation in identification of coding region. PLOS ONE. 2012;7 (2 ):e31621. PMC3280313. doi:10.1371/journal.pone.0031621 22355382
[15] Levine M, McDevitt RA, Meer M, et al. A rat epigenetic clock recapitulates phenotypic aging and co-localizes with heterochromatin. Elife. 2020;9 . PMC7661040. doi:10.7554/eLife.59201
[16] Hing B, Ramos E, Braun P, et al. Adaptation of the targeted capture methyl-seq platform for the mouse genome identifies novel tissue-specific DNA methylation patterns of genes involved in neurodevelopment. Epigenetics. 2015;10 (7 ):581–596. PMC4622595. doi:10.1080/15592294.2015.1045179 25985232
[17] Carey JL, Cox OH, Seifuddin F, et al. A rat methyl-seq platform to identify epigenetic changes associated with stress exposure. J Vis Exp. 2018. PMC6235597. 140 ). doi: 10.3791/58617-v
[18] Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-seq applications. Bioinformatics. 2011;27 (11 ):1571–1572. PMC3102221. doi:10.1093/bioinformatics/btr167 21493656
[19] Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9 :357–359. PMC3322381 22388286
[20] Li H, Handsaker B, Wysoker A, et al. Genome project data processing S. The sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25 (16 ):2078–2079. PMC2723002. doi:10.1093/bioinformatics/btp352 19505943
[21] Hansen KD, Langmead B, Irizarry RA. Bsmooth: from whole genome bisulfite sequencing reads to differentially methylated regions. Genome Biol. 2012;13 :R83. PMC3491411 23034175
[22] Colella S, Shen L, Baggerly KA, et al. Sensitive and quantitative universal Pyrosequencing™ methylation analysis of CpG sites. Biotechniques. 2003;35 (1 ):146–150. doi: 10.2144/03351md01 12866414
[23] Huang DW, Sherman BT, Tan Q, et al. DAVID bioinformatics resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 2007;35 (suppl_2 ):W169–75. PMC1933169. doi:10.1093/nar/gkm415 17576678
[24] Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25 (4 ):402–408. doi: 10.1006/meth.2001.1262 11846609
[25] Ye J, Coulouris G, Zaretskaya I, et al. Primer-blast: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics. 2012;13 (1 ):134. PMC3412702. doi:10.1186/1471-2105-13-134 22708584
[26] Chumley MJ, Catchpole T, Silvany RE, et al. EphB receptors regulate stem/progenitor cell proliferation, migration, and polarity during hippocampal neurogenesis. J Neurosci. 2007;27 (49 ):13481–13490. PMC6673089 doi:10.1523/JNEUROSCI.4158-07.2007 18057206
[27] Ma W, Xia C, Ling P, et al. Leukocyte-specific adaptor protein Grap2 interacts with hematopoietic progenitor kinase 1 (HPK1) to activate JNK signaling pathway in T lymphocytes. Oncogene. 2001;20 (14 ):1703–1714. doi: 10.1038/sj.onc.1204224 11313918
[28] Rao SP, Sancho J, Campos-Rivera J, et al. Human peripheral blood mononuclear cells exhibit heterogeneous CD52 expression levels and show differential sensitivity to alemtuzumab mediated cytolysis. PLOS ONE. 2012;7 (6 ):e39416. PMC3382607 receive compensation. Alemtuzumab is in clinical development by Genzyme. This does not alter the authors’ adherence to all the PLoS ONE policies on sharing the data and materials. doi:10.1371/journal.pone.0039416 22761788
[29] Milona A, Massafra V, Vos H, et al. Steroidogenic control of liver metabolism through a nuclear receptor-network. Mol Metab. 2019;30 :221–229. doi: 10.1016/j.molmet.2019.09.007 31767173
[30] Zhong S, Chevre R, Castano Mayan D, et al. Haploinsufficiency of CYP8B1 associates with increased insulin sensitivity in humans. J Clin Invest. 2022;132 (21 ):132. PMC9621133. doi:10.1172/JCI152961
[31] Berkel S, Tang W, Trevino M, et al. Inherited and de novo SHANK2 variants associated with autism spectrum disorder impair neuronal morphogenesis and physiology. Hum Mol Genet. 2012;21 :344–357. PMC3276277 21994763
[32] Jiang-Xie LF, Liao HM, Chen CH, et al. Autism-associated gene Dlgap2 mutant mice demonstrate exacerbated aggressive behaviors and orbitofrontal cortex deficits. Molecular Autism. 2014;5 (1 ):32. PMC4113140. doi:10.1186/2040-2392-5-32 25071926
[33] Rios I, Alvarez-Rodriguez R, Marti E, et al. Bmp2 antagonizes sonic hedgehog-mediated proliferation of cerebellar granule neurones through Smad5 signalling. Development. 2004;131 (13 ):3159–3168. doi: 10.1242/dev.01188 15197161
[34] Li H, Zhu Y, Morozov YM, et al. Disruption of TCF4 regulatory networks leads to abnormal cortical development and mental disabilities. Mol Psychiatry. 2019;24 (8 ):1235–1246. doi: 10.1038/s41380-019-0353-0 30705426
[35] Roy M, Kim N, Kim K, et al. Analysis of the canine brain transcriptome with an emphasis on the hypothalamus and cerebral cortex. Mamm Genome. 2013;24 (11–12 ):484–499. doi: 10.1007/s00335-013-9480-0 24202129
[36] Ackman JB, Ramos RL, Sarkisian MR, et al. Citron kinase is required for postnatal neurogenesis in the hippocampus. Dev Neurosci. 2007;29 (1–2 ):113–123. doi: 10.1159/000096216 17148954
[37] Bianchi FT, Gai M, Berto GE, et al. Of rings and spines: the multiple facets of Citron proteins in neural development. Small GTPases. 2020;11 (2 ):122–130. PMC7053930 doi:10.1080/21541248.2017.1374325 29185861
[38] Chen J, Guan Z. Function of oncogene mycn in adult neurogenesis and Oligodendrogenesis. Mol Neurobiol. 2022;59 :77–92. PMC8786763 34625907
[39] Cerase A, Pintacuda G, Tattermusch A, et al. Xist localization and function: new insights from multiple levels. Genome Biol. 2015;16 (1 ):166.PMC4539689 doi:10.1186/s13059-015-0733-y 26282267
[40] Lee JT, Davidow LS, Warshawsky D. Tsix, a gene antisense to xist at the X-inactivation centre. Nat Genet. 1999;21 (4 ):400–404. doi: 10.1038/7734 10192391
[41] Wamstad JA, Corcoran CM, Keating AM, et al. Role of the transcriptional corepressor bcor in embryonic stem cell differentiation and early embryonic development. PLOS ONE. 2008;3 (7 ):e2814. doi:10.1371/journal.pone.0002814 18795143
[42] Duan D, Goemans N, Takeda S, et al. Duchenne muscular dystrophy. Nat Rev Dis Primers. 2021;7 (1 ):13. doi: 10.1038/s41572-021-00248-3 33602943
[43] Mariani M, Zannoni GF, Sioletic S, et al. Gender influences the class III and V β-tubulin ability to predict poor outcome in colorectal cancer. Clinical Cancer Research. 2012;18 (10 ):2964–2975. doi: 10.1158/1078-0432.CCR-11-2318 22438565
[44] Manjaly ZR, Scott KM, Abhinav K, et al. The sex ratio in amyotrophic lateral sclerosis: a population based study. Amyotroph Lateral Scler. 2010;11 :439–442. PMC6485484.20225930
[45] Ozoguz A, Uyan O, Birdal G, et al. The distinct genetic pattern of ALS in Turkey and novel mutations. Neurobiology Aging. 2015;36 (4 ):1764 e9–e18. PMC6591733. doi:10.1016/j.neurobiolaging.2014.12.032
[46] Sun LL, Yang SL, Sun H, et al. Molecular differences in Alzheimer’s disease between male and female patients determined by integrative network analysis. J Cell Mol Med. 2019;23 (1 ):47–58. PMC6307813. doi:10.1111/jcmm.13852 30394676
[47] Odle AK, Allensworth-James ML, Akhter N, et al. Tropic role for leptin in the somatotrope as a regulator of POU1F1 and POU1F1-dependent hormones. Endocrinology. 2016;157 (10 ):3958–3971. doi: 10.1210/en.2016-1472 27571135
[48] de la Rocha Am, Sampron N, Alonso MM, et al. Role of SOX family of transcription factors in central nervous system tumors. Am J Cancer Res. 2014;4 :312–324. PMC4106650 25057435
[49] Matos B, Publicover SJ, Castro LFC, et al. Brain and testis: more alike than previously thought? Open Biol. 2021;11 (6 ):200322. PMC8169208.doi:10.1098/rsob.200322 34062096
[50] Schartl M, Schories S, Wakamatsu Y, et al. Sox5 is involved in germ-cell regulation and sex determination in medaka following co-option of nested transposable elements. BMC Biol. 2018;16 (1 ):16. PMC5789577. doi:10.1186/s12915-018-0485-8 29378592
[51] Wendt FR, Pathak GA, Singh K, et al. Sex-specific genetic and transcriptomic liability to neuroticism. Biol Psychiatry. 2023;93 (3 ):243–252. doi: 10.1016/j.biopsych.2022.07.019 36244801
[52] Knabl J, de Maiziere L, Huttenbrenner R, et al. Cell type- and sex-specific dysregulation of thyroid hormone receptors in placentas in gestational diabetes mellitus. Int J Mol Sci. 2020;21 (11 ):21. PMC7313460. doi:10.3390/ijms21114056
[53] Planells B, Gomez-Redondo I, Pericuesta E, et al. Differential isoform expression and alternative splicing in sex determination in mice. BMC Genomics. 2019;20 (1 ):202. PMC6419433. doi:10.1186/s12864-019-5572-x 30871468
[54] Chen H, Palmer JS, Thiagarajan RD, et al. Identification of novel markers of mouse fetal ovary development. PLOS ONE. 2012;7 (7 ):e41683. PMC3406020 doi:10.1371/journal.pone.0041683 22844512
[55] Irizarry RA, Ladd-Acosta C, Carvalho B, et al. Comprehensive high-throughput arrays for relative methylation (CHARM). Genome Res. 2008;18 :780–790. PMC2336799 18316654
[56] Meissner A, Mikkelsen TS, Gu H, et al. Genome-scale DNA methylation maps of pluripotent and differentiated cells. Nature. 2008;454 (7205 ):766–770. PMC2896277 doi:10.1038/nature07107 18600261
[57] Ziller MJ, Gu H, Muller F, et al. Charting a dynamic DNA methylation landscape of the human genome. Nature. 2013;500 (7463 ):477–481. PMC3821869 doi:10.1038/nature12433 23925113
[58] Yasukochi Y, Maruyama O, Mahajan MC, et al. X chromosome-wide analyses of genomic DNA methylation states and gene expression in male and female neutrophils. Proc Natl Acad Sci USA. 2010;107 (8 ):3704–3709. PMC2840519. doi:10.1073/pnas.0914812107 20133578
[59] Brockdorff N. Localized accumulation of xist RNA in X chromosome inactivation. Open Biol. 2019;9 (12 ):190213. PMC6936258 doi:10.1098/rsob.190213 31795917
[60] Jeon Y, Lee JT. YY1 tethers xist RNA to the inactive X nucleation center. Cell. 2011;146 :119–133. PMC3150513 21729784
[61] Cox OH, Song HY, Garrison-Desany HM, et al. Characterization of glucocorticoid-induced loss of DNA methylation of the stress-response gene Fkbp5 in neuronal cells. Epigenetics. 2021;16 (12 ):1377–1397. PMC8813076 doi:10.1080/15592294.2020.1864169 33319620
[62] Ito Y, Nakajima K, Masubuchi Y, et al. Expression characteristics of genes hypermethylated and downregulated in rat liver specific to Nongenotoxic Hepatocarcinogens. Toxicological Sciences. 2019;169 (1 ):122–136. PMC6484883 doi: 10.1093/toxsci/kfz027 30690589
[63] Kikuchi S, Takahashi Y, Ojiro R, et al. Identification of gene targets of developmental neurotoxicity focusing on DNA hypermethylation involved in irreversible disruption of hippocampal neurogenesis in rats. J Appl Toxicol. 2021;41 :1021–1037. PMC8247304 33150595
[64] Sallam M, Mysara M, Benotmane MA, et al. DNA methylation alterations in fractionally irradiated rats and breast cancer patients receiving radiotherapy. IJMS. 2022;23 . PMC9783664. 24 ):16214. doi: 10.3390/ijms232416214 36555856
[65] Kriaucionis S, Heintz N. The nuclear DNA base 5-hydroxymethylcytosine is present in Purkinje neurons and the brain. Science. 2009;324 (5929 ):929–930. PMC3263819 doi:10.1126/science.1169786 19372393
[66] Globisch D, Munzel M, Muller M, et al. Tissue distribution of 5-hydroxymethylcytosine and search for active demethylation intermediates. PLoS One. 2010;5 (12 ):e15367. PMC3009720 does not alter the authors’ adherence to the policies of PLOS ONE. doi:10.1371/journal.pone.0015367 21203455
[67] Perera A, Eisen D, Wagner M, et al. TET3 is recruited by REST for context-specific hydroxymethylation and induction of gene expression. Cell Rep. 2015;11 (2 ):283–294. doi: 10.1016/j.celrep.2015.03.020 25843715
