Parallel shift of DNA methylation and gene expression toward the mean in mouse spleen with aging

Age-related decrease in global DNA methylation level in mouse spleen

To investigate age-associated methylation drift (AMD), we performed RRBS (Supplementary Figure 1) on genomic DNA extracted from the spleens of 2-, 6-, 12-, and 23-month-old mice (2m, 6m, 12m, and 23m, respectively; n = 3 per group). We expected a genome-wide methylation change with age in the spleen because it is known that the function and microarchitecture of the spleen gradually deteriorate in older mice, leading to age-dependent impairments in functional immunity [30]. RRBS generated approximately thirty million reads per sample (Supplementary Figure 2A), of which 59% and 23% were uniquely and multiply mapped, respectively, after Bismark alignment (Supplementary Figure 2B). Supplementary Figure 2C depicts the quantitative profile of RRBS-captured MspI fragments to the number of inherent CpG sites in those fragments, highlighting the prevalence of MspI fragments nesting 2-7 CpGs as well as a proportional quantitative representation of the reads aligned to them. A significant proportion (>80%) of individual CpG sites had 10x coverage, which we set to the minimum acceptable depth for methylation ratios (Supplementary Figure 2D).

The mean methylation levels of whole CpG sites in samples decreased in older samples (Figure 1A), indicating a decline in global methylation level with age. The intra-group variation of mean methylation level was found to be at its minimum in the 6-month group and increased with age. The decrease was prominent in the 23m age group, which was in agreement with the principal component analysis (PCA) result using M-values (the log2 ratio of the frequencies of methylated CpG versus unmethylated CpG [31]) that revealed the 23m samples far removed from the other younger samples (Figure 1B). Therefore, the results indicated a global, age-associated methylation changes in mice of 23 months of age.

Region-based methylation analysis proves bidirectional AMD pattern

Violin plot using the M-values revealed that while the distribution of methylation ratios and their mean value were largely even and similar among the samples, the interquartile range (IQR (%) = Q3 – Q1), which is a measure of statistical dispersion of data, was shrunk in the 23m samples (Figure 1C and see also Supplementary Figure 3A for their normal distribution), indicating that the spread of the M-values is reduced in older samples. The result suggested the presence of AMDs with competing modes of action operating on different CpG clusters, which prompted us to conduct a “region”-based study for AMD as opposed to a singleton CpG-based approach. The term ‘region’ here encompasses spatially grouped CpG sites within a defined area that share a structural chromatin setting; henceforth, when we mention a region’s methylation level, we refer to the mean methylation of all the CpGs present in that region. We did this by assigning RRBS CpGs to one of three regional categories with respect to genome structure/function: gene promoters, gene bodies, and intergenic regions. Every CpG that was associated with a gene promoter (transcription start site (TSS) ± 1 kb in length) or gene body (encompassing TSS to transcription end site (TES) minus promoter region) was assigned to the gene, and the mean methylation level of the gene-associated CpGs was obtained collectively. The “intergenic block/region” is defined in the Method section. The frequency of regions with varying CpG counts in the three genome categories is depicted in Supplementary Figure 3B. As shown in Figure 1D, the majority of promoters (74% of 15,292 regions with 5 CpGs or more) were weakly methylated (β-value, the proportion of methylated CpG at a given locus, <5%) in 2m samples, whereas intergenic regions (82.4% of 57,679) were highly methylated (>75%) and gene-body regions (n = 16,554 genes with CpG >0.5 ea per kb) showed a mixed, bimodal methylation state. When the mean methylation levels were compared across ages, the promoters showed hypermethylation as age advanced, whereas the gene bodies and intergenic regions showed hypomethylation (Figure 1E). The extent of AMD was noticeably high from 12 months to 23 months. So, our region-based analysis proved a bidirectional AMD pattern in the mouse spleen. The result suggests that the AMD may cluster spatially as seen in cancer [32], refuting the claim that AMD implicates a larger proportion of singleton CpGs that exhibit solitary DNA methylation changes and are hence less likely to have functional significance.

Inter-individual heterogeneities in methylation level increase with age

The mean levels of DNA methylation across different genomic categories, including promoters, gene bodies, and intergenic regions, were found to be variable across chromosomes in this study (Figure 2A). Importantly, the standard deviations of methylation levels between individuals appeared to be relatively similar across chromosomes for each genomic category. We calculated the coefficient of variation (CV, the ratio of the standard deviation to the mean) of methylation levels across each chromosome among samples from the same age group to assess the change of inter-individual variability upon aging. The chromosome-wise CVs were persistently low on average in younger samples but significantly higher in older samples, suggesting a rise in heterogeneity with age (Figure 2B). It has been proposed that epigenetic drift, which leads to inter-individual epigenetic divergence, is determined by a lifelong accumulation of epimutations generated by environmental effects and/or intrinsic stochastic variables [33]. Our results with laboratory mice favor stochastic variables for AMD rather than environmental factors; clustered CpGs displaying epigenetic drift may undergo stochastic processes with aging, most likely due to the vulnerability of the region’s chromatin architecture to change.

AMD pattern at a specific CpG site is related to the density of peripheral CpGs

When promoters were classified by CpG density, those with thin CpGs (5-10 CpGs ea per 2 kb, n = 2,017) were relatively highly methylated (48.0%) at two months and became hypomethylated with age, whereas those with thick CpGs (>100 CpGs, n = 1,311) were lowly methylated (1.9%) and then hypermethylated with age (Figure 3A). It suggests that the AMD pattern is determined by the CpG density of the region. We sought to see if there were different types of AMDs in different CpG-density settings, regardless of genomic structure/function. To determine the CpG densities of adjacent sequences, we counted the number of CpGs present in a 200-bp region on either side of a single CpG (Figure 3B, top). The DNA methylation levels at the same density CpG sites were box-plotted, which revealed three qualitatively distinct groups, one with Q1 quartile values larger than 50%, another with Q3 values less than 50%, and the others otherwise (Figure 3B, bottom). Based on this methylation level-based criterion, the density of each CpG site was classified as low (L, 0–5 CpGs), high (H, >12), or moderate (M, 6-12), and each CpG could be assigned to one of nine classifications with two capitals: HH, HM, HL, MH, MM, ML, LH, LN, and LL. The H-labeled CpGs–HH, HM/MH, and HL/LH–were generally undermethylated whereas the L-labeled CpGs–LL and LM/ML–were heavily methylated (Figure 3C). It is commonly considered that HH-CpGs are primarily found in CpG islands (CGIs) and are typically weakly methylated, whereas LL-CpGs are located outside of CGIs (or in ‘open-ocean,’ as the CGI-related description suggests) and are comparatively highly methylated. The HH and HL groups had the most and the fewest CpGs, respectively, of the 1.4 x 10⁶ CpGs available in our RRBS data (n = 487,612 and 8,822; Figure 3D).

We noticed that CpG groups with different density labels exhibit different AMD patterns (Figure 3E). DNA methylation increased with age for the HH and HM/MH CpGs, but declined for the L-/M-labeled CpGs. In the classifications exhibiting a hypomethylation-type AMD, 2m samples were slightly but consistently less methylated than 6m samples, which may indicate that global DNA methylation patterns in mice are still being established at two months of age [34]. Interestingly, the HL/LH group differed from the other H-labeled groups such as MH/HM group and resembled the non-H-labeled (L/M) group. We guess that the methylation level difference at two months between the HM/MH and HL/LH groups (13.7% vs. 27.7%, respectively) could explain the different AMD patterns. Scatter plots of 6m, 12m, and 23m versus 2m samples demonstrated that while the HH and LL classifications were equivalent until 12 months of age, they split into two distinct subsets at 23 months of age, hypermethylation (HH-CpGs) and hypomethylation groups (LL-CpGs; Supplementary Figure 4).

AMD pattern is determined by the region’s methylation level

While the bulk of promoters were undermethylated, a minor percentage was considerably methylated. We wondered if these highly methylated promoters experienced the same hypermethylation with age as the lowly methylated promoters. We classified the promoters based on their methylation level (0.1 β-value intervals). Ten % (n = 1,508) of promoters were quite substantially methylated (>50%) and showed a hypomethylation rather than hypermethylation with age (Figure 4A). A similar result was observed from the analysis of gene-body regions (Supplementary Figure 5). The result indicates that when certain regions are initially substantially methylated, they tend to be hypomethylated regardless of the genome structure/function (Figure 4B), indicating that the AMD pattern is methylation level dependent.

The regions in each genomic structure/function were grouped based on the correlation of their methylation levels with sample ages (R; 0.1 intervals each, a total of 20 intervals ranging from -1.0 to +1.0). As shown in Figure 4C, there was a negative relationship between the groups’ R-values and their mean methylation levels at 2 months; as the groups had higher methylation levels, their R-values were more negative (or hypomethylated), whereas as the groups had lower methylation levels, they were more positive (or hypermethylated). It indicates that the pattern of AMD is more closely related to the region’s methylation level. This was most noticeable in the intergenic group (coefficient of determination, R², = 0.982). The promoter category groups with positive R-values showed an exceptionally poor correlation (R² = 0.087) with the methylation levels, probably due to their overall low-lying methylation levels with a marginal difference (~5% difference in the intervals of 0.3 to 1.0) between the intervals.

We then looked into promoters that had an unusual pairing of CpG density and methylation level, such as the one with low CpG density and low methylation level and the one with high CpG density and high methylation level. The incidence of hypermethylation was three times that of hypomethylation (454 vs. 141 promoters; |R-value| >0.8; Figure 4D) for promoters with a low CpG density (15>CpG>4) and a low methylation level (<10%). There were only a few promoters (7 of 1258 promoters) with a high CpG density (>100 CpGs) and high methylation level (>20%), insufficient for a reciprocal analysis. On the other hand, for promoters with a low methylation level (<10%, n = 13,076), hypermethylation prevailed over hypomethylation in all CpG-density intervals (Figure 4E). Together, our results indicate that the methylation level of the area is more important than CpG number in determining AMD.

Age-associated transcriptional drift reduces the gene-to-gene disparity of transcriptional activity

The steady expansion of an epigenomic disorder impairs transcriptome homeostasis and increases transcriptional noise. The SiNG-PCRseq approach, which reliably measures transcript levels, particularly of low-set expression level genes [29], was used to investigate the transcriptional activity of genes in the aged spleen. In addition to the spleen samples used for methylation analysis, we included additional spleen samples of ages 2m (n=2), 6m (n=3), and 12m (n=2) for gene expression analysis. We also used a reference sample, consisting of a mixture of mouse and rat genomic DNAs, to correct amplification biases. We used a set of primer pairs that can amplify 424 target exonic sequences from 257 genes of both mouse and rat, including 175 genes identified as epi-drivers [35]. Rat genomic DNA was utilized as a spike-in in multiplexed SiNG-PCRseq to accurately quantify mouse cDNA sequences compared to matched rat gDNA sequences. To ensure equal amplification and to establish detectable flags between the cDNA-derived reads and the rat genomic DNA reads after sequencing, each primer pair was designed to contain one or more nucleotide differences between the amplicons of different species. After SiNG-PCRseq, we acquired 2.2×10⁸ reads from triplicate libraries of 20 samples (3.7×10⁶ ±0.6×10⁶ reads/library), with around 76% of them mapping perfectly to our amplicon reference. The M/R ratio of each exonic target, comparing mouse read counts to rat counts, underwent quality filtering for low read numbers and extensive amplification biases, resulting in a set of 288 exonic targets from 206 genes. After filtering, the data were processed for amplification bias correction and normalization, allowing us to evaluate the expression level of each exonic sequence across all samples (Supplementary File 1). For instance, among DNA methyltransferases (Dnmt) genes, Dnmt1 was less expressed with age whereas Dnmt3b and Dnmt3l were more expressed at 23 months (Supplementary Figure 6).

The violin plots revealed that the distribution of expression levels among the genes was the broadest in the 6m group, then shrunk significantly in the 12m group, and finally narrowed in the 23m group (Figure 5A). As a result, the 6m group had the largest variance in the mean gene expression levels, but it was decreased to less than one-fifth in the 23m group (Figure 5B). Furthermore, as shown by the lowered slope of scatter plot in Figure 5C, those genes that were weakly expressed in young spleen became more transcribed in old, whereas those genes which were relatively well expressed in young tended to diminished in old. We also observed a strong correlation between the third quartile values of CpG methylation levels and the first quartile values of gene expression levels across all samples analyzed (R² = 0.62; Figure 5D, left). This suggests a potential link between the global loss of methylation in highly methylated regions and the increase in transcription of lowly expressed genes. We also found that the number of CpG sites with no methylation decreased with age and was moderately correlated with the first quartile values of gene expression (R² = 0.57; Figure 5E). This implies another link between methylation gain in lowly methylated regions and the decrease of expression levels. We previously identified a similar age-associated bidirectional change in transcriptional activity in mouse CD4-positive T cells [36]. Together, our results suggest that age-associated bidirectional transcriptional drift occurs, resulting in a reduction in disparities in expression levels between strongly and weakly expressed genes. This age-associated transcriptional drift pattern mirrors the AMD pattern, in that genes with a lower-level transcript experience ‘hyper-transcription’ whereas those with a high-level transcript show ‘hypo-transcription’. Meanwhile, there was no obvious association (R² = 0.04) between methylation alteration and expression change (23m/2m) in the epi-driver genes (Figure 5F, see Discussion).

Mouse and nucleic acid preparation

Genomic DNAs and RNAs were simultaneously extracted from the spleen tissues using All Prep® DNA/RNA Mini Kit (Qiagen, Germany) following the manufacturer’s instruction.

Reduced representation bisulfite sequencing (RRBS)

Bisulfite-converted MspI-restricted genomic fragments were molecularly indexed and sequenced using an in-house developed RRBS method (Supplementary Figure 1). In this method, we devised two different adaptors with a CpG overhang in one end. Adaptor-A carries a molecular barcode and the binding site for Illumina R1 sequencing primer. Adaptor-B has the binding site for Illumina R2 sequencing primer. All the cytosines on the Illumina primer binding sites of two adaptors were protected from bisulfite conversion as they were priorly methylated during synthesis. Fragmented genomic DNAs (100ng) were purified using Zymo DNA clean-up Kit (Zymo Research, USA) and ligated with the adaptors (2 pmol/each) at 25° C for 30 min using Blunt/TA Ligase Master Mix (NEB, USA). The recessive part of the ligates was then filled in using Tag DNA polymerase (Solgent, Korea) at 72° C for 30 min, which concomitantly removes the unligated adaptor strand. In this Fill-in step, dCTP was replaced with methyl-dCTP. Bisulfite conversion was performed for the adaptor-attached MspI fragments using the EZ DNA Methylation-Gold™ Kit (Zymo Research, USA) following the manufacturer’s instructions. After bisulfite conversion, MspI fragments were subsequently amplified by two successive rounds of PCR to generate sequencing-ready libraries. First PCR was performed using primers, P1 and P2, for 10 cycles of 95° C for 30sec, 60° C for 30sec, and 72° C for 30 sec. Second PCR was performed using index primers, idxP1, and idxP2, for 10 cycles of 95° C for 30 sec, 60° C for 30 sec, and 72° C for 30 sec. Massive parallel sequencing was performed on Illumina Hiseq platform with 151 bp paired-end reading.

SiNG-PCRseq

Total RNA (1ug) was converted to cDNA using iScript cDNA Synthesis Kits (Bio-Rad, USA). SiNG-PCRseq was performed using the cDNA and the rat genomic DNA as a competitor array for the target genes as described previously. Briefly, we composed 10 primer pools using 424 primer pairs and used each in multiplexed PCR with 15 ng of mouse cDNA and 2 ng of rat gDNA as templates in the following conditions: 95° C for 15 min, 45 cycles of 95° C/20 sec, 57° C /40 sec, 72° C/1 min, followed by a final extension at 72° C for 5 min. The Amplicons were pooled together for each sample and purified using Expin purification kit (GeneAll, Korea), 5’-end phosphorylated using T4 polynucleotide kinase (NEB) by incubating at 37 for 30 min, ligated with NGS adapters using T4 ligase (NEB) at 25 for 2 hours and amplified using idxP1 and idxP2 to generate sequencing libraries with following conditions: after 98° C for 15 min, 20 cycles of 98° C/20 sec, 68° C/30 sec, and 72° C/1 min, followed by a final extension at 72° C for 5 min. Three independent libraries were prepared for each sample to ensure reproducibility. Additionally, a library of genome mixture of mouse and rat was generated in triplicate and used to correct for any amplification biases in the target DNA. Massive parallel sequencing was performed on Illumina Nextseq platform with 151 bp single-end reading.

Sequence analysis

The read data from the RRBS sequencing were processed using trim_galore (version 0.6.6) to remove the adaptors. The cleaned reads were then aligned to the mouse reference genome (mm10) using the Bismark package (version 0.22.3). The methylation status of each CpG site was determined using ad hoc Perl scripts. Clonal copies were identified by the presence of the same molecular barcode and were removed by deduplication. The methylation levels of CpG sites were collected, and only those with at least 10 deduplicated reads in each sample were selected for further analysis. To compare the global and regional methylation patterns between samples, only CpG sites with qualified methylation levels from at least 13 samples were used in the analysis. The promoter region was defined as the region of 2 kilobases surrounding the transcription start site (TSS ± 1 kilobase), the gene body region was defined as the region from the start to end described in the annotation files, and the intergenic region was defined as any other region. To capture the collective influence of neighboring CpG sites, the intergenic CpG block was operationally defined as the region where any two neighboring CpG sites are less than 1 kilobase apart from each other in the intergenic region.

The read data from SiNG-PCRseq were analyzed essentially as described. Briefly, adaptor-trimmed sequences were identified for their target and species origin by aligning them to the target DNA sequences using the BLAST program. The count data for each target sequence from each sample were obtained and subjected to a quality filter. Targets were excluded if any reference library had fewer than 100 read counts or if any sample library had read 300 counts. The count data for each target amplicon from each sample were utilized to determine the mouse fraction. Targets with extensive amplification bias in the reference sample (<0.3 or >0.7) were excluded from the analysis. The amplification biases present in the mouse fractions of the remaining targets were corrected by employing the mouse fractions of the reference sample as previously described. The bias-corrected fractions were then converted to relative quantities of the rat amplicon and normalized by using the average relative quantities of the library. The final expression values of a sample were determined by averaging the relative quantities obtained from the triplicate libraries of the sample.

Data availability

Sequencing data were deposited under GSE203289.