Genome-wide Methyl-Seq analysis of blood-brain targets of glucocorticoid exposure

ABSTRACT

Chronic exposure to glucocorticoids (GCs) can lead to psychiatric complications through epigenetic mechanisms such as DNA methylation (DNAm). We sought to determine whether epigenetic changes in a peripheral tissue can serve as a surrogate for those in a relatively inaccessible tissue such as the brain. DNA extracted from the hippocampus and blood of mice treated with GCs or vehicle solution was assayed using a genome-wide DNAm platform (Methyl-Seq) to identify differentially methylated regions (DMRs) induced by GC treatment. We observed that ∼70% of the DMRs in both tissues lost methylation following GC treatment. Of the 3,095 DMRs that mapped to the same genes in both tissues, 1,853 DMRs underwent DNAm changes in the same direction. Interestingly, only 209 DMRs (<7%) overlapped in genomic coordinates between the 2 tissues, suggesting tissue-specific differences in GC-targeted loci. Pathway analysis showed that the DMR-associated genes were members of pathways involved in metabolism, immune function, and neurodevelopment. Also, changes in cell type composition of blood and brain were examined by fluorescence-activated cell sorting. Separation of the cortex into neuronal and non-neuronal fractions and the leukocytes into T-cells, B-cells, and neutrophils showed that GC-induced methylation changes primarily occurred in neurons and T-cells, with the blood tissue also undergoing a shift in the proportion of constituent cell types while the proportion of neurons and glia in the brain remained stable. From the current pilot study, we found that despite tissue-specific epigenetic changes and cellular heterogeneity, blood can serve as a surrogate for GC-induced changes in the brain.

KEYWORDS:

• Biomarker
• blood-brain correlation
• Cushing's disease
• DNA methylation
• epigenetics
• glucocorticoids
• tissue heterogeneity

Abbreviations and Acronyms
GC	=	glucocorticoid
DNAm	=	DNA methylation
DMR	=	differentially methylated region
CORT	=	corticosterone

Introduction

Tools that can assess functions of specific genes and processes in the brain through the use of peripheral tissues may be immensely useful for the diagnosis of disease. However, only a few studies have correlated measurements occurring in both tissues.¹ Furthermore, to our knowledge, no study has identified the underlying mechanisms that link the cross-tissue changes caused by disease-causing factors, since most of such factors are either unknown or poorly characterized. Administration of such factors and assessment of their impact on the central nervous system (CNS) and peripheral tissues can provide us with underlying biological principles that can facilitate biomarker discovery.

One such factor is a class of steroid hormones known as glucocorticoids (GCs). Chronic exposure to GCs, either by iatrogenic administration, tumors, or as a part of the neuroendocrine response to stress, leads to numerous health complications, including cardiovascular disease, obesity, and diabetes.^2-5 The CNS is particularly vulnerable to GCs, as prolonged GC exposure is associated with cognitive decline,^6,7 and psychiatric disorders such as anxiety,⁸ depression,⁹ and bipolar disorder.¹⁰ Many studies have examined the physiological, biochemical, and behavioral consequences of exposure to stress and GCs in both animals and humans. In animals, the ability to control the duration and magnitude of GC exposure and the accessibility of relevant tissues make them an ideal model to study the effects of GC exposure. In humans, where studies of stress and GCs and stress-related disorders are the most relevant, studies of the brain have used either postmortem specimens or imaging studies that do not provide (epi)genomic resolution. Clinical studies performed in patients invariably involve the use of peripheral tissues that may or may not provide useful information on GC-induced brain impairments. Efforts to directly examine such impairments are hindered by the highly inaccessible and irreplaceable nature of brain tissues.

Previously, we have demonstrated that the degree of DNA methylation (DNAm) changes in the stress-response gene Fkbp5 in blood correlated strongly with mean 30-day plasma GC levels as well as GC-dependent weight changes in the spleen and visceral fat.¹¹ Interestingly, DNAm changes in the blood can also accurately reflect 30-day GC exposure and dose-dependent alterations in Fkbp5 DNAm and gene expression in the hippocampus,¹² as well as correlating with anxiety-like behavior.¹¹ This demonstrated the feasibility of measuring epigenetic changes in a peripheral tissue to assess GC-dependent epigenetic changes and gene function in the brain, and prompted us to identify additional hippocampal genes whose epigenetic alterations may be reflected in blood. To this end, we used the genome-wide, targeted capture Methyl-Seq platform to identify GC-induced, differentially methylated regions (DMRs) in mouse blood and brain tissues. For the brain, we chose the hippocampus, as exposure to stress and GCs has been shown to hinder hippocampal function,^6,13,14 its volume,^15,16 and its development,^17,18 which in turn can act as risk factors for the development of psychiatric disorders.

Materials and methods

Animals

At 5 weeks of age, male C57BL/6J mice (n = 12 in each group; Jackson Laboratories, Bar Harbor, ME) were given ad libitum access to solutions containing the rodent stress glucocorticoid corticosterone (Sigma-Aldrich, St. Louis, MO; 100 μg/ml with 1% ethanol; “CORT” group) or 1% ethanol (“VEHICLE” group) in place of their normal drinking water for 4 weeks. The method of providing CORT in the drinking water has been established by several studies that have sought to examine the role of corticosterone in mediating the stress response in rodents.^19-22 In the first cohort for Methyl-Seq, 4 pooled samples (n = 2 mice per sample and 2 samples per group) from each group were chosen from animals whose 4-week mean plasma CORT levels were closest to the average of each group. Rest of the individual samples in each group were used for bisulfite pyrosequencing. A second cohort (n = 12 each group) was treated with the same VEHICLE and CORT solutions for 4 weeks to obtain mRNA for expression studies and cortex tissues for fluorescence-activated cell sorting (FACS). Finally, a third cohort (n = 6 in each group) was treated with the same VEHICLE and CORT solutions for one week to obtain white blood cells for FACS. At the end of the treatment period, animals were killed, tissues harvested, and their brains frozen on powdered dry ice and subsequently stored at -80^oC. All procedures were approved by the Institutional Animal Care and Use Committee at 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.

Blood collection

Whole blood (20 μL) was collected weekly (0900 h) from each mouse into heparinized glass capillary tubes through a small nick at the tip of the tail. Blood samples were centrifuged at 4^oC and plasma was collected and frozen at −80^oC for further analysis of total (free and bound) plasma CORT by radioimmunoassay, according to the manufacturer's instructions (MP Biomedicals, Costa Mesa, CA). All samples were run in duplicates with an intra-assay coefficient of variance (CV) of 3.9% and inter-assay CV of 7.5%. After the 4-week treatment period, animals were killed, and trunk blood (∼250 μL) was collected. Trunk blood samples were incubated with 3 volumes of Ack Lysing Buffer (Quality Biological, Gaithersburg, MD) to lyse red blood cells, centrifuged to collect white blood cells, and subsequently processed for genomic DNA and mRNA extraction or FACS. Whole blood cell (WBC) count of lymphocytes and neutrophils were also determined by the Phenotyping Core at the Johns Hopkins Department of Molecular and Comparative Pathobiology. FACS of trunk blood and WBC were performed after one week of CORT treatment.

Dissection of the cortex and dorsal hippocampus

Frozen mouse brains were sectioned using a cryostat, and 400 μm sections were mounted on glass slides. Needles (19 gauge; 0.686 mm inner diameter and 1.086 mm outer diameter) were used to dissect multiple (>10) punches from the cortex (bregma 1.18 through 0.78 mm) and 2–3 punches from the dorsal hippocampus (bregma -0.98 through -2.0 mm). Punched tissues were stored at -80°C in the Tissue & Cell Lysis Buffer (Epicentre Biotechnologies, Madison, WI) or the RLT Buffer (Qiagen, Germantown, MD) until processed for genomic DNA (gDNA) or total RNA, respectively.

Fluorescence-activated cell sorting of the cortex and blood tissues

For the cortex, nuclei were prepared from fresh-frozen cortex tissues by detergent and physical lysing using the Nuclei Extraction Buffer²³ and Dounce homogenization, respectively. Cortex homogenates were layered on top of a 1.8 M sucrose gradient and subjected to ultracentrifugation for 2 h in a Beckman L7–65 Ultracentrifuge with a SW 41 Ti rotor. Pelleted nuclei were washed once with ice-cold PBS and labeled with AlexaFluor-conjugated primary antibodies against the neuron-specific marker NeuN (Millipore, Billerica, MA) in PBS with 0.5% BSA and 10% goat serum. Labeled nuclei were subsequently sorted in the BD FACS Aria at the Johns Hopkins Flow Cytometry Core.²⁴ To assess blood, freshly-collected samples were treated with the Ack Lysing Buffer for red blood cell lysis, washed twice with PBS, and labeled with APC-conjugated CD3e, PerCP-conjugated CD45R, and PE-conjugated Ly-6G/Ly-6C antibodies (BD Biosciences, Franklin Lakes, NJ) to simultaneously label and isolate T-cells, B-cells, and neutrophils, respectively, by FACS. Sorted cells were pelleted and resuspended in the Tissue & Cell Lysis Buffer for subsequent gDNA extraction.

Genomic DNA extraction

Genomic DNA (gDNA) from mouse hippocampus, cortex, and blood was isolated with the Masterpure DNA Purification Kit, according to the manufacturer's instructions (Epicentre Biotechnologies). Concentration of gDNA was determined using the Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA).

mRNA extraction

mRNA from mouse hippocampus, cortex, and blood was isolated with the RNeasy Micro Kit according to the manufacturer's instructions (Qiagen, Germantown, MD). Concentration and quality of the mRNA were determined using the Agilent TapeStation 2200, and all of the RIN (RNA integrity number) was greater 8.5.

SureSelect Mouse Methyl-Seq

SureSelect Target Enrichment System (Mouse Methyl-Seq, Agilent Technologies, Santa Clara, CA) was implemented to identify CORT-induced epigenetic changes in the blood and brain tissues. The Methyl-Seq platform targets regulatory regions across the genome, as described previously.²⁵ Briefly, 1 μg of gDNA from each animal was pooled from 2 animals to form one sample, and 2 samples per tissue were sequenced for each treatment group. Each gDNA sample was sheared using the Covaris sonicator (Covaris, Woburn, MA) to yield 170–230 bp DNA fragments. These fragments were end-repaired, 3’-adenylated, and further ligated with methylated primers. Following hybridization to biotinylated, plus-strand DNA-complementary RNA library “baits," precipitation from the solution using streptavidin-coated magnetic beads, and RNase-digestion of the “baits,” captured DNA was bisulfite-converted using the EZ DNA Methylation Gold Kit (Zymo Research, Irvine, CA). Subsequently, DNA samples were PCR-amplified using sample-specific indexed (“barcoding”) primers to allow for multiplexing and sequenced on the Illumina Hi-Seq2000. 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 guanine-cytosine content, per base N content, sequence length distribution, sequence duplication levels, overrepresented sequences, adaptor content, and kmer content. Reads were trimmed using Trim Galore v0.3.7. Default parameters were used, and one base pair was trimmed off at the end of all paired-end reads to improve paired-end mapping. If adaptor contamination was observed, the standard Illumina adapters were trimmed off at the end of all paired-end reads. FASTQC post-trimming was rerun to perform quality control to ensure that the trimming step did not produce any adverse side effects.

Reads were mapped to the Mus musculus assembly (Genome Reference Consortium Mouse Build 38/ mm10), which was produced by the Mouse Genome Reference Consortium, using Bismark version 0.13.0 with Bowtie 2 version 2.1.0. Briefly, for alignment purposes Bismark converts all Cs to Ts (in forward reads) and all Gs to As (in reverse reads) before 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 Ts and As 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 in the seed” 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. Bismark divides all cytosines into 4 categories: cytosines followed by guanines (CpGs), cytosines followed by non-guanines followed by guanines (CHGs), cytosines followed by at least 2 non-guanines (CHHs), and cytosines followed by Ns (CNs). Analysis for the current study focused on CpGs. 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 the samples. Using the genomic location and matrices consisting of M (methylation) and Cov (read coverage of CpGs) values, BSseq smoothes the M-values across CpGs using the BSmooth function, computes t-statistics between groups of samples using the function BSmooth.tstat, and establishes the threshold levels of the t-statistics to identify 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 the VEHICLE and CORT groups, CpGs with at least 10x coverage across all samples were used. The t-statistics were 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 calculated. The threshold for including consecutive CpGs in a DMR was the 99^th percentile (tail ends) of the t-statistic distribution. Regions with 3 or more CpGs and greater than or equal to 10% mean methylation difference between VEHICLE vs. CORT groups were ranked and displayed by “areaStat” and subjected to further investigation. To assess the statistical uncertainty for each observed DMR, we also used permutation techniques. We permuted the outcome variable 1000 times and re-calculated the t-statistics for each permutation. The dmrFinder was run each time to produce a set of null areas. Empirical P-values were calculated as the fraction of null areas greater than each observed area. For example, an observed area greater than 95% of the areas obtained from the permutation exercise was assigned an empirical P-value of 0.05. To account for the multiplicity problem introduced by the genome-wide screening, false discovery rates (FDRs) were computed based on these P-values.

Pathway analysis

The hypergeometric test for enrichment^26-29 was used to identify KEGG pathways that were enriched in the DMR lists for the hippocampus, blood, and both tissues. Genes associated with all of the DMRs in each list, including genes of DMRs common to both tissues, were compared with genes representing each of the KEGG pathways. The hypergeometric test was used to analyze KEGG genes observed and expected for each DMR list to calculate enrichment and P-value scores. Q-values were determined by the FDR method for multiple-test correction of the enrichment P-values. The same analysis was also performed using a subset of DMRs that were P < 0.05 in the DMR lists for the hippocampus, blood, and both tissues (results shown in the Table S6).

Bisulfite pyrosequencing

DMRs identified from Methyl-Seq experiments were further replicated by bisulfite pyrosequencing of the PCR products, which measures methylation variation at >90% precision.³⁰ Genomic DNA was bisulfite-converted using the EZ DNA Methylation Gold Kit. Design and implementation of pyrosequencing assays have been described previously,^31,32 and the primers used for each assay are included in the Supplementary Tables.

Quantitative real-time PCR

Expression levels of several genes associated with the DMRs identified in hippocampus, cortex, and blood were assessed by quantitative real-time PCR (qPCR). cDNA was generated from the mRNA using the QuantiTect Reverse Transcription Kit (Qiagen). TaqMan probes targeting the transcripts of DMR-associated genes and that of β-actin were used on the 7900HT Fast Real-time PCR System (Applied Biosystems, Foster City, CA). Cycle threshold (Ct) values for each gene, which are indicators of mRNA abundance, were normalized against those of the housekeeping gene β-actin to determine relative expression.

Results

Plasma glucocorticoid levels in mice

To identify common epigenetic targets between the blood and brain tissues, we treated mice for 4 weeks with 100 μg/mL of CORT or VEHICLE solutions (n = 12 each). Mean plasma CORT levels derived from 4 weekly measurements were 538.2 ± 67.5 ng/mL (mean ± SEM) for the CORT-treated group and 50.0 ± 7.3 ng/mL for the VEHICLE-treated group. Weekly levels are shown in Fig. S1. Genomic DNA isolated from the hippocampus and blood of 4 animals from each group was processed for the SureSelect Methyl-Seq Target Enrichment platform followed by high-throughput sequencing on the Illumina HiSeq2000. For each tissue, 2 samples each consisting of pooled gDNA from 2 CORT-treated animals were compared with those of the VEHICLE-treated animals. DNA was pooled since relatively small hole-punched tissues of the hippocampus from one animal were insufficient to provide the recommended amount (∼3.5 μg) for Methyl-Seq. DNA from blood was similarly pooled to provide consistent results using the same animals from which the hippocampal gDNA was pooled. Genomic DNA processed from the remaining 8 animals in each group was used for replication.

Methyl-Seq and Next-generation sequencing

Across 4 blood samples (representing gDNA from 8 animals) that were sequenced, we obtained an average of 40.5 million paired-end reads of which an average of 22.6 million reads were deduplicated and uniquely mapped paired-end reads. This corresponded to 1.7 million individual CpG dinucleotides with an average read-depth coverage of 30X when confining our analysis to CpGs that have at least 10 sequenced reads. Similarly, we obtained an average of 39.9 million paired-end reads from the hippocampal samples of which an average of 24.7 million reads were deduplicated and uniquely mapped. For the hippocampus, we analyzed 1.9 million CpGs with an average read-depth coverage of 36X after at least 10X reads per CpG threshold had been imposed. In addition, the sequenced reads corresponded to an on-target rate of 85% of the targeted genomic regions for blood and 87% for hippocampus, which represented 143 Mb and 149 Mb of the mouse genome, respectively. The metrics for the high-throughput sequencing of the blood and hippocampal samples are shown in Table 1.

Table 1. Read information from Methyl-Seq.

Sequencing Metrics	BLOOD	HIPPOCAMPUS
Paired End Reads (PER)	40,514,128	39,852,125
Uniquely Mapped Paired End Reads (UMPER)	31,349,404	30,021,816
Mapping Efficiency (UMPER/PER)	77.4%	75.4%
Duplicate Reads (% of UMPER)	27.4%	18.2%
Deduplicated UMPER	22,631,126	24,699,490
Average Read Depth Coverage (X) (ARDC)	12.5	15.0
CpGs (N)	5,747,028	6,071,788
ARDC (X) of CpGs	10.5	13.0
CpGs with at least 10 reads (N)	1,689,017	1,943,698
ARDC (X) of CpGs with at least 10 reads	29.5	35.8
ON TARGET CpGs with at least 10 reads (N)	1,473,324	1,647,390
ON TARGET ARDC (X) of CpGs with at least 10 reads	30.3	37.8
ON TARGET (PER with 1 or more Base Pair overlap with Probe Target Regions) (UMPER)	19,293,523	21,559,486
% ON TARGET (of Deduplicated UMPER)	85.2%	87.3%
ON TARGET (Total Bases Mapped) Mb	143.0	149.0
ON TARGET Average Read Depth Coverage (X) (ARDC)	27.0	28.8

Methyl-Seq analysis and DMR validation

Using BSseq and BSmooth, methylation values at each of the CpGs were computed, smoothed estimates of raw methylation values across CpGs for each sample were determined, t-statistics for the smoothed estimates between the treatment groups were calculated, and DMRs were identified. Following BSseq analysis of CpGs with several criteria (at least 10X coverage, 3 or more CpGs, ≥ 10% mean methylation difference between VEHICLE vs. CORT groups and empirical P-value ≤ 0.05), we identified 5,491 and 5,365 DMRs in the blood and hippocampus, respectively. The top 10 DMRs for blood and hippocampus are shown in Table 2. From the extensive list of DMRs, we selected several biologically relevant DMRs to independently validate the Methyl-Seq platform by pyrosequencing. Some of the genes associated with the validated DMRs are involved in tissue-specific development, psychiatric disorders, and/or glucocorticoid signaling, as expected from mice treated with CORT. The goal was to demonstrate the robustness of the Methyl-Seq platform.

Table 2. Top 10 most significant genes regulated by CORT in blood and hippocampus.

Location of Blood DMRs	Gene	Distance from Gene	N of CpGs	DMR Width	areaStat	Vehicle (%)	CORT (%)
Chr17: 28,898,895-28,899,766	4930539E08Rik	in_gene	28	872	205.1	43.3	76.7
Chr10: 17,725,094-17,727,734	Cited2	1,867	32	2,641	182.6	76.1	29.9
Chr8: 121,570,435-121,571,212	Fbxo31	in_gene	34	778	165.6	69	21.6
Chr11: 117,116,679-117,117,586	Sec14l1	in_gene	28	908	162.6	92.3	51.2
Chr7: 29,212,131-29,212,609	Catsperg1	in_gene	23	479	162.6	79.3	36.8
Chr7: 15,963,944-15,964,287	Ehd2	in_gene	23	344	161.5	35.6	71.2
Chr14: 33,189,451-33,190,294	Arhgap22	27,371	29	844	151	31.7	86.3
Chr13: 103,592,675-103,593,009	Srek1	148,939	33	335	150.1	56.5	89.8
Chr15: 82,166,760-82,168,036	Srebf2	in_gene	30	1,277	146.1	45.1	70.6
Chr15: 85,231,797-85,232,719	Fbln1	in_gene	26	923	141.6	84.3	42.9
Location of Hippocampal DMRs	Gene	Distance from Gene	N of CpGs	DMR Width	areaStat	Vehicle (%)	CORT (%)
Chr5: 75,834,926-75,835,856	Kdr	98,343	32	931	111.1	34.2	23.6
Chr12: 12,936,763-12,937,453	Mycn	in_gene	30	691	86.2	47.8	35.9
Chr8: 20,281,184-20,281,505	6820431F20Rik	in_gene	26	322	81.1	59.6	47.4
Chr8: 8,621,164-8,622,036	Efnb2	in_gene	19	873	80.4	56.3	37.9
Chr11: 120,280,582-120,280,938	Bahcc1	in_gene	25	357	77.3	38.1	21.5
Chr13: 74,533,968-74,534,503	Zfp825	53,912	20	536	76	41.8	30.1
Chr4: 141,477,474-141,477,726	Spen	in_gene	16	253	72.3	39.1	19.8
Chr8: 20,596,693-20,597,331	Gm21119	18,233	19	639	64	55.8	42.9
Chr7: 13,006,038-13,006,361	Zbtb45	in_gene	14	324	60.8	61.8	40.3
Chr7: 45,465,579-45,466,110	Bax	in_gene	16	532	59.6	54.9	42.5

For the hippocampus (HPC), we chose 2 DMRs within the EfnB2 and Shank3 loci that had been implicated by Methyl-Seq (Fig. 1A and C). EfnB2 encodes a member of the ephrin family of tyrosine kinase receptors that along with REELIN are crucial for neurodevelopment.³³ EfnB2 is also involved in thymus development and T-cell differentiation.^34,35 Shank3 encodes a postsynaptic scaffold protein that contains multiple ankyrin repeat domains and also plays an important role in neuron function and development.^36,37 Interestingly, numerous genetic and animal model studies have linked this gene to autism.^38-40 Bisulfite pyrosequencing validation of several CpGs in each DMR for VEHICLE vs. CORT animals are shown (n = 8 per group, Fig. 1B and D). For the blood, we chose 2 DMRs within the Nfil3 and Ltbr loci (Fig. 1E and G). Nfil3, or nuclear factor interleukin 3 regulated, encodes a transcription factor that represses Per1/Per2 gene expression to control circadian rhythm,^41,42 and plays important roles in T-cell development and cytokine production in the blood.^43-45 Ltbr encodes the lymphotoxin β receptor that is a member of the tumor necrosis factor superfamily and is involved in lymphoid organ development.^46,47 Interestingly, both Nfil3 and Ltbr are demonstrated targets of GC signaling,^48,49 and validation of their CpGs by bisulfite pyrosequencing is shown in Fig. 1F and H. Exhaustive DMR lists for both blood and hippocampus are shown in Tables S1 and S2, respectively. Also, in light of the subtle epigenetic changes in the hippocampus, we tested a few genomic regions that are not on the DMR list as positive and negative controls for GC exposure. Pyrosequencing results for previously identified DMRs¹² in Fkbp5 intron 5 and intron 1 (positive controls for hippocampus and blood, respectively), as well as promoter regions of Crh, Fkbp4, and Gfap (negative controls) in hippocampus are shown in Figs. S2A-E. Genes associated with the DMRs showed an increase in expression in the CORT-treated animals for Ltbr, Nfil3, and Fkbp5 in the blood (Fig. S3A) and EfnB2, Shank3, and Fkbp5 in the hippocampus (Fig. S3B). On the other hand, we observed no change in the expression of the hippocampal Gfap, consistent with its static DNAm levels between the 2 treatment groups.

Figure 1. Identification and replication of DMRs implicated by Methyl-Seq. (A) Methyl-Seq graphical output of EfnB2 DMR in the hippocampus according to BSmooth. Pink shaded region represents the DMR, gray vertical lines represent CpG dinucleotides, and the dots represent % DNAm estimates for VEHICLE- (blue) or CORT-(red) treated samples. The lines are the smoothed estimates used for DMR calculations. (B) Pyrosequencing replication (n = 8 per group) of sequential CpGs in the DMR of EfnB2. Similar graphical output (C) and pyrosequencing replication (D) for a hippocampal DMR in Shank3. Graphical output (E, G) and pyrosequencing replication results (F, H) are also shown for blood DMRs in Nfil3 and Ltbr, respectively. Bar graphs in the right panels are represented as mean ± SEM. *P<0.05, **P<0.01, and ***P<0.001.

Pathway analysis of epigenetic targets of glucocorticoids in the blood and brain

Given the numerous DMRs that exist between VEHICLE vs. CORT samples in both tissues, we performed KEGG pathway analysis to identify pathways that may be over-represented. Interestingly, metabolic and cancer pathways were the top most significant entries for both tissues (Table 3). Pathways involved in myeloid leukemia and endocytosis were highly significant in the blood DMRs, whereas pathways involved in WNT signaling and axon guidance were highly significant in the hippocampal DMRs. Given the presence of common pathways such as those involved in MAPK, neurotrophin, and insulin signaling, we performed an additional analysis of genes that were commonly represented in both DMR lists. As expected, similar pathways were implicated as in the analysis of individual tissue. The results suggest that a significant subset of genes in each tissue commonly undergo epigenetic changes following GC exposure. Comprehensive tables of all FDR-significant pathways in blood, hippocampus, and both tissues are included in Tables S3–5. We also performed pathway analysis only on DMRs that met nominal P-value significance (P < 0.05), and those results are included in Table S6.

Table 3. Distinct pathways regulated by CORT in blood and hippocampus.

KEGG Pathways in Blood	Genes in Category	Genes Observed	Genes Expected	Enrichment	P-value	Q-value
Metabolic pathways	1,184	152	40.9	3.7	2.7 × 10^−44	2.4 × 10^−40
Pathways in cancer	325	51	11.2	4.5	2.1 × 10^−19	1.9 × 10^−15
Acute myeloid leukemia	57	21	2	10.7	1.2 × 10^−16	1.1 × 10^−12
Neurotrophin signaling	131	27	4.5	6	6.9 × 10^−14	6.4 × 10^−10
Insulin signaling pathway	137	27	4.7	5.7	2.2 × 10^−13	2.0 × 10^−09
Endocytosis	220	33	7.6	4.3	1.6 × 10^−12	1.5 × 10^−08
Adipocytokine signaling	68	18	2.4	7.7	1.1 × 10^−11	1.0 × 10^−07
Phosphatidylinositol signaling	78	19	2.7	7.1	1.5 × 10^−11	1.4 × 10^−07
MAPK signaling pathway	268	35	9.3	3.8	2.0 × 10^−11	1.8 × 10^−07
Purine metabolism	168	27	5.8	4.7	3.3 × 10^−11	3.0 × 10^−07
KEGG Pathways in HPC	Genes in Category	Genes Observed	Genes Expected	Enrichment	P-value	Q-value
Metabolic pathways	1,184	220	78.9	2.8	4.0 × 10^−44	3.7 × 10^−40
Pathways in cancer	325	104	21.7	4.8	7.2 × 10^−43	6.6 × 10^−39
MAPK signaling pathway	268	74	17.9	4.1	2.3 × 10^−26	2.1 × 10^−22
Wnt signaling pathway	154	50	10.3	4.9	1.1 × 10^−21	1.0 × 10^−17
Calcium signaling pathway	178	53	11.9	4.5	6.0 × 10^−21	5.5 × 10^−17
Neurotrophin signaling	131	42	8.7	4.8	3.1 × 10^−18	2.9 × 10^−14
Type II diabetes mellitus	49	23	3.3	7.1	8.8 × 10^−15	8.1 × 10^−11
Axon guidance	131	37	8.7	4.2	2.8 × 10^−14	2.6 × 10^−10
Adipocytokine signaling	68	26	4.5	5.7	6.6 × 10^−14	6.1 × 10^−10
Insulin signaling pathway	137	37	9.1	4.1	1.3 × 10^−13	1.2 × 10^−09
KEGG Pathways in Blood and HPC	Genes in Category	Genes Observed	Genes Expected	Enrichment	P-value	Q-value
Pathways in cancer	325	42	5.81	7.22	2.06E^−23	9.78E^−19
Metabolic pathways	1,184	68	21.18	3.21	6.89E^−17	3.27E^−12
Acute myeloid leukemia	57	17	1.02	16.67	1.23E^−16	5.84E^−12
Pancreatic cancer	71	14	1.27	11.02	2.95E^−11	1.40E^−06
Chronic myeloid leukemia	74	14	1.32	10.58	5.31E^−11	2.52E^−06
Chemokine signaling pathway	185	20	3.31	6.04	1.83E^−10	8.69E^−06
Adipocytokine signaling	68	13	1.22	10.69	2.25E^−10	1.07E^−05
Wnt signaling pathway	154	18	2.75	6.53	4.10E^−10	1.95E^−05
Insulin signaling pathway	137	17	2.45	6.94	4.89E^−10	2.32E^−05
Neurotrophin signaling	131	16	2.34	6.83	1.99E^−09	9.45E^−05

Genomic organization of the DMRs

DMRs identified in the blood and hippocampus were classified based on their locations with respect to genes and CpG islands. More than 50% of the DMRs were located within intronic regions for both tissues, with promoters making up about 10% of the DMRs and intergenic regions (upstream and downstream of genes) making up more than 35% (Fig. 2A). In addition, more than 77% of the blood and hippocampal DMRs were located well beyond the CpG islands, as approximately 20% of the DMRs fell within island shores, commonly defined as regions within 3 kilobases flanking the CpG islands,⁵⁰ and less than 2% of the DMRs were located within the islands themselves (Fig. 2B). DMR locations are consistent with those of previous studies that reported tissue-specific DMRs predominantly in intronic and intergenic regions that are not CpG islands.^31,51 Also, studies that have performed similar analyses with targets of GCs in vitro⁵² and in vivo⁵³ or HPA axis dysregulation in clinical samples⁵⁴ have primarily implicated intronic regions.

Figure 2. Genomic locations of the blood and hippocampal DMRs with respect to genes and CpG islands. (A) Pie charts for both blood and hippocampal DMRs show the percentage of DMRs that fall in promoters, exons, introns, and intergenic regions. (B) Pie charts for both blood and hippocampal DMRs show the percentage of DMRs that fall in CpG islands, shores, and non-CpG-associated regions (Other).

Characteristics of the DMRs

Given the thousands of DMRs identified, we also asked whether they could provide us with insights into GC biology. We first observed that many of the DMRs in both tissues were reduction-of-methylation events in the CORT group, which prompted us to score these events for the entire list. Of the 5,491 DMRs identified in blood, 70.7% were DMRs that lost methylation following CORT treatment. Similarly, 72.0% of the 5,365 DMRs identified in the hippocampus were also loss-of-methylation events. We then retraced our steps in the analytical pipeline, examined percent methylation values determined for individual CpGs between the treatment groups, and constructed a boxplot of all of the CpG methylation values for each chromosome. Loss of DNAm was observed in all of the chromosomes for both blood (Fig. 3A) and hippocampus (Fig. 3B), with the exception of the Y-chromosome in the hippocampus that exhibited overall increase in DNAm. However, we noted that the Y-chromosome in the hippocampus also showed an unusual number of outlier CpGs in the CORT-treated animals, which underwent drastic loss of methylation. We also noted that the relatively small number of mitochondrial CpGs that are represented in Methyl-Seq are largely devoid of DNAm in both tissues. Finally, we asked whether the magnitude of methylation changes between the treatment groups is different between gain vs. loss of methylation events. We constructed a boxplot to visualize the absolute magnitude of the methylation changes segregated by direction. In blood, DMRs that lost methylation with CORT treatment underwent a median loss of 38.3%, whereas DMRs that gained methylation with CORT underwent a median gain of 25.4% (Fig. 3C). The mean differences of the DMRs between the treatment groups were highly significant (P = 2.2 × 10⁻¹⁶). In contrast, we observed a more subtle, but significant difference in the magnitude of change between gain vs. loss of methylation events in the hippocampus. DMRs that lost methylation with CORT treatment underwent a median loss of 14.9%, compared to 14.1% in the DMRs that gained methylation with treatment (Fig. 3D). Despite the similarity in magnitude, the mean differences of the DMRs were still significant (P = 0.001).

Figure 3. Box plots of CORT-induced DNAm changes in blood and hippocampus. (A) Box plot of average methylation values (n = 2 samples or 4 mice per group) for individual CpGs with at least 10 reads. The plots for blood DNA are displayed for each chromosome for VEHICLE-(blue) vs. CORT-(red) treated samples. (B) Similar box plot for the hippocampus is shown. (C) DMRs between VEHICLE- and CORT-treated animals were segregated based on the direction of DNAm change. Box plots of absolute mean % methylation differences for loss-of-methylation and gain-of-methylation DMRs in blood are shown. (D) Similar box plots for the hippocampal DMRs are shown.

Common DMRs between the blood and hippocampus

In light of our previous finding where we observed CORT-induced loss of methylation of Fkbp5 occurring at different loci in blood and brain,^12,32 we asked whether genomic coordinates of across-tissue DMRs that were incident in the same genes overlapped with one another. From the DMR lists, we identified 3,095 DMRs that mapped to the same genes in both blood and hippocampus, with 1,853 DMRs (∼60%) showing DNAm changes in the same direction. We then constructed a histogram of all intragenic blood-hippocampus DMR distances. We found that most of the DMRs did not overlap when the distances were represented as absolute distances in basepairs (Fig. S4A) or as ratios normalized by the length of each RefSeq gene (Fig. S4B). In fact, only 209 DMRs (<7%) overlapped in genomic coordinates between the 2 tissues, and of these, 136 DMRs (65%) changed DNAm in the same direction. We then performed bisulfite pyrosequencing on the tissue-specific DMRs belonging to 3 common genes whose expression was affected by GCs according to literature. Neuron navigator 2 (Nav2, Fig. 4A and B) plays a role in neuronal growth and migration, and its tissue-specific DMRs are 145,528 basepairs apart. Recent work in corneal fibroblasts demonstrates that Nav2 is upregulated by the synthetic GC dexamethasone.⁵⁵ Similar loss-of-methylation patterns were observed in both of the tissue-specific DMRs within Slc38a2 (also known as Sat2 or Snat2, Fig. 4C and D), a widely expressed amino acid transporter that is upregulated by GCs.⁵⁶ Its DMRs are 1,053 bp apart. We also tested the DMRs within Nr1d1 (also known as Rev-erbA-α), a gene involved in the development and circadian regulation of tissues such as the liver, brain, and adipose tissues. Its DMRs are merely 334 bp apart. However, unlike the other 2 genes examined, the pyrosequencing confirmed that the methylation changes in blood and brain occur in opposite directions (Fig. 5E and F). Interestingly, Nr1d1 is differentially regulated by GCs, as studies have reported GC-induced upregulation in the striatum⁵⁷ and downregulation in the liver.⁵⁸ Finally, in the hippocampus, we examined the identical genomic coordinates of the DMRs that were implicated from the blood Methyl-Seq experiment. As peripheral DNA is often used to identify regions that are associated with psychiatric disorders, we sought to test the assumption that the same regions in the brain undergo disease-relevant epigenetic changes. Blood DMRs for Nav2, Slc38a2, and Nr1d1 tested on the hippocampal DNA showed little evidence of DNAm changes as observed in blood DNA (Figs. S5A-C).

Figure 4. Pyrosequencing replication of blood and hippocampal DMRs that fall on identical genes. (A) and (B) Bisulfite pyrosequencing for the DMRs in the Nav2 gene for the hippocampus and blood are shown, respectively. (C) and (D) Bisulfite pyrosequencing for the DMRs in the Slc38a2 gene for the hippocampus and blood are shown, respectively. (E) and (F) Bisulfite pyrosequencing for the DMRs in the Nr1d1 gene for the hippocampus and blood are shown, respectively. Note that DMRs for Nr1d1 in blood (gain of DNAm) and brain (loss of DNAm) are predicted to occur in opposite directions. Bar graphs in all of the panels are represented as mean ± SEM.*P < 0.05, **P < 0.01, and ***P < 0.001.

Figure 5. FACS to delineate the contribution of different cell types to the observed DNAm differences. (A) Cell count by FACS after one-week VEHICLE or CORT treatment. Contribution from each of the cell types are presented as a percentage of the total cell count obtained from T-cells, B-cells, and neutrophils. (B) Intronic CpGs within Fkbp5 previously shown to undergo CORT-induced loss of DNAm in whole blood were tested. Results for one of the CpGs is shown to undergo substantial methylation change in only T-cells. (C) Induction of Fkbp5 in whole blood by CORT treatment. Fkbp5 expression has been normalized by that of Actb (β-actin). (D) Testing of FACS on the mouse cortex. Cortex tissues were separated into neuronal (NeuN+) and non-neuronal (NeuN-) fractions by FACS. Bisulfite pyrosequencing of the conserved Gfap promoter shows substantially lower methylation in the (NeuN-) fraction, consistent with it's expression in the glia. (E) Absence of CORT-induced changes cellular heterogeneity in the cortex. Bisulfite pyrosequencing of the Gfap promoter on unsorted cortex tissues shows no difference in methylation between the VEHICLE- and CORT-treated animals. (F) CORT-induced changes in DNAm in the cortex. CpG-4 within the fifth intron of Fkbp5 undergoes methylation loss in the unsorted cortex tissues. (G) More substantial change in DNAm can be seen in the same intronic region of Fkbp5 in the NeuN+ fraction following FACS. (H) Induction of Fkbp5 in the cortex by CORT treatment. Fkbp5 expression has been normalized by that of Actb (β-actin). Bar graphs in all of the panels are represented as mean ± SEM.*P < 0.05, **P < 0.01, and ***P < 0.001.

Tissue heterogeneity in the brain and blood

Chronic exposure to GCs is associated with significant changes in the populations of lymphocytes, neutrophils, and platelet counts in blood^59-63 as well as astrogliogenesis in the hippocampus,⁶⁴ and mesenchymal-to-epithelial transition in some tissue-derived cell lines.⁶⁵ We asked whether changes in the populations of cell types and different methylation patterns in each cell type played a role in the observed methylation differences between the VEHICLE- and CORT-treated animals.

First, to demonstrate the presence of CORT-induced changes in cell types in the blood, mice were treated for one week, during which the blood underwent substantial changes in DNAm (∼70% of total).¹¹ Trunk blood (n = 8 per group) was then used to perform complete blood count (CBC) to identify CORT-induced changes in the percentage of lymphocytes (%LY), which includes T-cells and B-cells, and neutrophils (%NE) in the blood. We observed a substantial reduction in the %LY from 81.0 ± 2.5% to 41.3 ± 2.7% (mean ± SEM, P = 1.46 × 10⁻⁷) and an increase in the %NE from 12.8 ± 1.4% to 52.0 ± 2.7% (P = 2.23 × 10⁻⁸) in the CORT-treated animals. In a separate experiment, we took the same whole blood from the 2 groups and determined the percent composition of the T-cells, B-cells, and neutrophils by labeling them with cell type specific antibodies and counting the different sorting events during FACS. Among these 3 cell types, we observed reductions in the percentages of T-cells and B-cells and an increase in the percentage of neutrophils in the CORT-treated group (Fig. 5A) similar to changes observed during CBC.

This raised the possibility that some of the observed methylation changes in the blood may be confounded by an influx of neutrophils with lower methylation content.⁶⁶ To measure the effect of heterogeneity on DNAm, DNA was extracted from the sorted cells and subjected to bisulfite pyrosequencing to interrogate the Intron 1 GC response element (GRE) CpGs of Fkbp5 (also tested in Fig. S2B in whole blood). We observed significant loss of methylation in the T-cells, but not in the B-cells or neutrophils (Fig. 5B) and an increase in expression of Fkbp5 in the whole blood (Fig. 5C). Our results demonstrated that while blood underwent significant CORT-induced changes in cell type composition, only T-cells, which are the predominant cell type in the mouse blood,⁶⁷ exhibited GC-induced changes in DNAm.

We then asked whether similar events occur in the brain. Since the hippocampal punches did not yield sufficient amount of tissue for FACS, we instead used punches from the cortex, where larger tissue material can be obtained and where we have also demonstrated CORT-induced transcriptional and epigenetic changes in the Fkbp5 gene (unpublished). Sucrose gradient-treated nuclei were labeled with the NeuN/Fox-3 antibody, FAC-sorted, and processed for bisulfite pyrosequencing. First, there were no significant differences in the ratio of NeuN+ nuclei to the total nuclei sorted between the VEHICLE- and CORT-treated animals (n = 4, 78.2 ± 2.3% vs. 80.5 ± 1.2%, P = 0.43), suggesting that there were only trivial CORT-induced changes in cell type composition in the cortex. To further demonstrate the absence of a change in the cell type proportion in the brain, we also implemented an epigenetic approach. Examination of a highly-conserved promoter region of the gene encoding the glial fibrillary acidic protein (GFAP) in the VEHICLE-treated cortex showed substantial methylation differences across 6 CpGs between NeuN- and NeuN+ fractions (15.4–54.9% DNAm increase in NeuN+, P = 3.5 × 10⁻⁴–1.2 × 10⁻⁶), consistent with expression of Gfap in the NeuN-negative glia (Fig. 5D). No differences in DNAm were observed in Gfap between the VEHICLE- and CORT-treated groups in unsorted cortex tissues (Fig. 5E) or in both NeuN- and NeuN+ fractions (data not shown). If CORT treatment led to changes in the populations of neuronal and non-neuronal fractions, then we would be able to observe substantial changes in the Gfap promoter in the unsorted tissue. On the other hand, a subtle difference in DNAm was observed in an intronic CpG of Fkbp5 before FACS in cortex (Fig. 5F), and the differences became more substantial in the CORT-treated, NeuN+ fractions after FACS (Fig. 5G). Specifically, we observed 7.9% (P = 0.02) DNAm difference at CpG-4 between the VEHICLE- and CORT-treated groups in the unsorted cortex, and FACS of the cortex tissue unmasked greater differences in DNAm at multiple CpGs: CpG-1 (15.6%, P = 0.002); CpG-3 (16.2%, P = 0.002), CpG-4 (15.6%, P = 0.009). The epigenetic changes in the cortex were also associated with an expression change in Fkbp5 (Fig. 5H), similar to that observed in blood. There was an absence of epigenetic and expression changes in Gfap following CORT treatment in the hippocampus (Figs. S2E and S3B, respectively). Taken together, these observations suggest that neuronal fractions undergo CORT-induced DNAm changes without substantial changes in neuronal populations in the cortex.

Discussion

Although there are many studies that have identified and characterized peripheral “markers” of disease-relevant processes, basic underlying principles that govern across-tissue analysis are largely uncharacterized. In this study, we capitalized on our previous finding that GCs directly influence DNAm to identify common epigenetic targets of GCs in the mouse blood and hippocampus. Our goal was to determine whether the same genes underwent GC-induced epigenetic changes in both tissues and in doing so, to assess the feasibility of using blood as a potential surrogate for brain tissue.

To this end, we used a genome-wide methylation platform called Methyl-Seq to identify thousands of GC-induced DMRs across the genome in mice that were treated with CORT. To date, several notable studies have compared tissue-specific epigenetic signatures between tissue types, in the context of aging^68,69 or neurodevelopment.⁵¹ To our knowledge, the current study is one of the first that documents treatment-induced, genome-wide epigenetic changes across tissue types.

We made several observations regarding epigenetics and GC biology. First, based on the direction of methylation change of the DMRs, GCs promote loss-of-methylation events. This is surprising given that GC-bound GR can act as both activator and repressor of transcription in a locus-specific way.^70,71 Also, CORT-induced loss-of-methylation events undergo larger magnitude change in methylation than gain-of-methylation events. Such loss-of-methylation events may be the primary way that GCs facilitate transcriptional activation. For gene repression, especially at pro-inflammatory genes, alternative mechanisms such as modifications of histones⁷² or GRE-independent mechanisms such as co-repression with NF-κB⁷³ may play a bigger role.

Second, pathway analysis of the blood and hippocampal DMRs are consistent with the role of GCs in metabolism, obesity, type-2 diabetes, and the immune system, as metabolic, adipocytokine, insulin, and chemokine signaling pathways have been implicated. Despite the significant contribution of cellular heterogeneity to the observed methylation differences between the VEHICLE- and CORT-treated animals in the blood, implication of similar pathways as those in the brain supports the use of whole blood to identify CORT-induced epigenetic events. Of particular interest are pathways such as WNT, insulin, neurotrophin, and JAK-STAT signaling pathways that play crucial roles in brain function. As such, pathways and DMRs that are common to these tissues first serve to affirm the biological role of GCs at the epigenomic level by implicating anticipated pathways, but more importantly make it feasible to test activity of brain-specific genes that may be similarly regulated in the blood. For instance, the WNT signaling regulates neurogenesis in the hippocampus⁷⁴ and is negatively regulated by GCs.⁷⁵ It may be possible to assess the degree of GC-induced, DMR-mediated inhibition of WNT genes in the hippocampus by assessing the methylation levels of their counterparts in the blood.

Third, another observation made from the DMRs that map to common genes is that only a small percentage of the DMRs overlap in genomic coordinates. This observation has important implications. While these overlapping DMRs have the potential to serve as peripheral markers of GC-exposure in the brain, their relatively low occurrence, at least in our study, may impede their practical use. Also, different tissues undergo GC-induced changes at different loci. Previous studies on identification of tissue-specific DMRs³¹ suggest the presence of tissue-specific methylation differences scattered throughout any given gene, and it is thought that these regions undergo epigenetic changes to influence gene expression, as observed in our study. A more important implication is that across-tissue epigenetic markers may exist for environmental factors, as in our case with GC exposure, but the epigenetic marks may be located at different regions within the gene. In other words, assessing the same genomic region for GC-induced differences in different tissue types may reduce the probability of finding a peripheral marker. Nevertheless, we state that these tissue-specific DMRs can still serve as useful biomarkers, so long as the DNAm changes in one tissue are associated with DNAm or expression changes in the brain.¹²

Lastly, we found that the magnitude of DNAm change was far greater in blood than in the brain. Since reference methylomic information is not yet available for the different leukocytes in the mouse blood as it is for humans, we used FACS to examine the degree to which cellular heterogeneity may play a role in the observed DMRs. In blood, we observed the expected influx of neutrophils and efflux or elimination of the T- and B-cells from circulation. Results obtained from both the whole blood cell count and FACS support previously published studies that have documented demargination of neutrophils⁷⁶ and either efflux or apoptotic destruction of T- and B-cells following GC administration.⁷⁷ Methylation analysis of an intronic region of the Fkbp5 gene previously shown¹¹ to undergo CORT-induced change revealed that only the T-cells were affected. Taken together, the large methylation changes observed at the Fkbp5 locus and elsewhere are likely due to a combination of demargination of neutrophils with low CpG methylation, elimination of T- and B-cells from the circulation, and loss of methylation in specific cell types such as the T-cells. Additional studies are warranted to further examine the interaction of these factors and characterize their roles in the CORT-treated blood.

On the other hand, strong support for GC-induced changes in the proportion of cell types was not observed in the cortex or the hippocampus. Using the more abundant cortex tissue, we first demonstrated that there are substantial methylation differences at the Gfap promoter between neuronal and non-neuronal fractions and that methylation levels of the Gfap promoter in the unsorted cortex tissues did not change in the CORT-treated animals. A significant reduction of DNAm in the CORT-treated samples would have suggested a reduction in the neuronal population that contributed the higher methylation levels in the VEHICLE-treated animals. More importantly, similar observations were made in the unsorted hippocampal tissues. This is noteworthy, since the dentate region of the hippocampus is one of the few regions of the developed brain thought to undergo active neurogenesis.⁷⁴ Our findings justify either physically separating or statistically adjusting for the proportion of the cell types of interest to unmask potential GC-induced epigenetic changes confounded by cellular heterogeneity.

As the current study is exploratory in nature, it has several limitations. First and foremost, we did not measure expression of genes for all of the DMRs that we implicated by Methyl-Seq. Although many of the DMRs in both blood and hippocampus were associated with GC target genes, comprehensive expression studies such as RNA-Seq are necessary to derive the functional significance of each DMR. Unfortunately, our approach of performing hole-punch dissection of the dorsal hippocampus, a subregion highly susceptible to stress and GCs,^78,79 limited the amount of tissue available for RNA-Seq. However, using qPCR, we demonstrated CORT-induced change in the expression levels of several genes that were associated with candidate blood and hippocampal DMRs. Further, in light of the ∼70% of the DMRs showing loss-of-methylation events in both tissues, it is imperative to reconcile this with the number of genes that are up- or downregulated in the same tissues.

Second, DNA from 2 animals was pooled for each sample analyzed by Methyl-Seq. Consequently, this approach affected our ability to perform useful FDR calculations, as 2 pooled DNA samples (from 4 animals) per treatment group per tissue did not provide adequate sample size. Nonetheless, DMR rankings based on the areaStat parameter provided an adequate measure of confidence that allowed us to replicate the results using additional samples. There are advantages to pooling samples, as several studies have used this strategy to improve statistical power and reduce variations that may arise from individual samples.^80-82 In the current study, our primary rationale for pooling samples was to obtain sufficient amount of starting material for Methyl-Seq.

Third, it is not clear whether any of these DMRs become occupied by the GC receptor complex following chronic GC treatment. It is possible that many of these DMRs may arise indirectly by proteins encoded by genes that are the primary targets of GC signaling. However, there are no comprehensive genome-wide data sets for targets of GC signaling for the mouse genome, and many of the GC receptor binding sites identified in the human⁸³ or rat^52,53 genomes do not show strong cross-species conservation. Chromatin immunoprecipitation (ChIP) assays performed on the mouse T-cells and hippocampus are needed to verify the DMRs as such.

Fourth, the animals were exposed to high levels of GCs. The level of GC exposure is similar to that of Cushing's disease patients or patients with autoimmune diseases such as rheumatoid arthritis or lupus who receive GCs. As the GC levels exceed those experienced during stress exposure, it is unlikely that identical magnitude of methylation changes occurs in stressed animals, although it is plausible that the epigenetic signatures (loci) could be similar. Despite the high dose of CORT, the current study serves as a proof-of-concept demonstration that prolonged exposure to GCs, the primary agent in the stress response, can cause genome-wide epigenetic changes. Use of stress paradigms, such as psychosocial stress⁸⁴ or chronic variable stress,⁸⁵ are needed to characterize the epigenetic changes that occur due to stress exposure.

Fifth, treatment with CORT led to a significant shift in the populations of specific cell types and served as a confounding factor in our measurement of DNAm levels in the blood. As demonstrated with methylation changes at the Fkbp5 locus, a significant shift in the cell type composition of blood likely made a significant contribution in the larger magnitude change observed in DMRs that underwent loss-of-methylation vs. gain-of-methylation. Currently, we are unable to determine to what extent the shift in cell type composition may play a role on other DMRs beyond Fkbp5. FACS followed by Methyl-Seq of the major cell types in blood needs to be performed in mice to properly control for the contribution of each cell type in the observed changes in methylation. In addition, our blood Fkbp5 data showed that GC-induced changes occurred only in T-cells. However, this observation may only apply to Fkbp5. It is plausible that at other loci, only B-cells or neutrophils undergo GC-induced changes in methylation. Additional DMRs in different cell types need to be tested to ascertain whether our observation at Fkbp5 can be generalized.

Lastly, we have yet to determine the portion of DMRs that are comprised of hydroxymethylation events, as bisulfite treatment of DNA cannot distinguish methylcytosine (mC) from hydroxymethylcytosine (hmC) without additional modifications. This issue is especially relevant for the hippocampal DMRs, as hmC modifications are highly abundant in the brain.⁸⁶ Additional chemical modifications using potassium perruthenate (KRuO₄, oxBS-Seq) or enzymatic modifications of DNA using β-glucosyltransferases and TET (TAB-Seq) are needed to resolve mC from hmC.⁸⁷

Despite the limitations, we identified several genes where blood could serve as a surrogate for GC-induced epigenetic changes in the brain. This finding warrants additional studies to examine and establish across-tissue epigenomic and transcriptomic correlations in greater detail.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Additional information

Funding

This study was funded by NIH grants MH101392 (RSL), AA020890 (GSW) and MH97150 (KLT), and support from the following awards and foundations: a NARSAD Young Investigator Award, Margaret Ann Price Investigator Fund, the James Wah Mood Disorders Scholar Fund via the Charles T. Bauer Foundation, and the Project Match Foundation (RSL).