Demethylation and derepression of genomic retroelements in the skeletal muscles of aged mice

Abstract Changes in DNA methylation influence the aging process and contribute to aging phenotypes, but few studies have been conducted on DNA methylation changes in conjunction with skeletal muscle aging. We explored the DNA methylation changes in a variety of retroelement families throughout aging (at 2, 20, and 28 months of age) in murine skeletal muscles by methyl‐binding domain sequencing (MBD‐seq). The two following contrasting patterns were observed among the members of each repeat family in superaged mice: (a) hypermethylation in weakly methylated retroelement copies and (b) hypomethylation in copies with relatively stronger methylation levels, representing a pattern of “regression toward the mean” within a single retroelement family. Interestingly, these patterns depended on the sizes of the copies. While the majority of the elements showed a slight increase in methylation, the larger copies (>5 kb) displayed evident demethylation. All these changes were not observed in T cells. RNA sequencing revealed a global derepression of retroelements during the late phase of aging (between 20 and 28 months of age), which temporally coincided with retroelement demethylation. Following this methylation drift trend of “regression toward the mean,” aging tended to progressively lose the preexisting methylation differences and local patterns in the genomic regions that had been elaborately established during the early period of development.

murine skeletal muscles by methyl-binding domain sequencing (MBD-seq). The two following contrasting patterns were observed among the members of each repeat family in superaged mice: (a) hypermethylation in weakly methylated retroelement copies and (b) hypomethylation in copies with relatively stronger methylation levels, representing a pattern of "regression toward the mean" within a single retroelement family. Interestingly, these patterns depended on the sizes of the copies. While the majority of the elements showed a slight increase in methylation, the larger copies (>5 kb) displayed evident demethylation. All these changes were not observed in T cells. RNA sequencing revealed a global derepression of retroelements during the late phase of aging (between 20 and 28 months of age), which temporally coincided with retroelement demethylation. Following this methylation drift trend of "regression toward the mean," aging tended to progressively lose the preexisting methylation differences and local patterns in the genomic regions that had been elaborately established during the early period of development.
Comprehensive DNA methylation analyses have been facilitated by genome-wide sequencing technology that enables the detection of numerous methylation alterations that are present in aged tissues.
However, methylome research on skeletal muscle aging in mammals has been scarce, and no consensus on the relationship between them has been reached. Zykovich et al. (2014) compared the DNA methylation patterns of human skeletal muscles from healthy young (~20 years of age) and aged (~80 years of age) males and observed predominant hypermethylation patterns in the aged skeletal muscles. Day et al. (2013) explored DNA methylation changes with age in various tissues and found that skeletal muscle tissues exhibited the strongest association with tissue-specific genes and the least overlap in the list of age-associated CpG dinucleotides (CpGs) with other tissues. Jin et al. (2014) compared the skeletal muscle methylation of young and middle-aged pigs and found that there was a global hypomethylation trend in middle-aged pigs. In contrast to the former two human studies that utilized bead chip arrays, which are used to perform limited examinations of the selected CpGs around the genic regions, the latter pig study utilized a genome-wide analysis technique called methylated DNA immunoprecipitation (MeDIP), but their analyses were also restricted only to the genic regions.
Consequently, all of these skeletal muscle methylome studies on age-associated methylation drift lacked information about the nongenic regions, including retroelements.
Aging is associated with the loss of repressive heterochromatin integrity, leading to the abnormal activation of gene expression in those regions Villeponteau, 1997). This heterochromatin loss with age, a hallmark of aging occurring in species ranging from yeast to humans (Haithcock et al., 2005;Larson et al., 2012;Villeponteau, 1997), makes the genome unstable, and this genome instability is furthered by the increased expression and transposition of a variety of retroelements (Pal & Tyler, 2016). Sedivy et al. proposed an "aging-by-transposition" model, in which the retroelements and their transposases are the primary causes of structural dysregulation in the genome to manifest the aging phenotypes . It is notable that retroelement expression and transposition are also associated with cancer development (Kang, 2018;Lee et al., 2012) and neurodegenerative diseases (Reilly, Faulkner, Dubnau, Ponomarev, & Gage, 2013;. With some exceptions, aging in mammals is more commonly associated with CpG hypomethylation, particularly at repeats (Bjornsson et al., 2008;Bollati et al., 2009;Bormann et al., 2016;Horvath, 2013;Jintaridth & Mutirangura, 2010). Since the loss of DNA methylation at genomic repeats will intensify the risk of transposition, it may be, at least partly, causally related to the loss of heterochromatin during aging. However, the age-associated methylomic changes at genomic repeats have not been extensively studied. In the present study, we analyzed the methylomes of aged skeletal muscle tissues with a keen interest in retroelement sequences; our aim was to evaluate the epigenetic stability during skeletal muscle aging. We used methyl-CpG-binding domain (MBD)-based capture sequencing, which is relevant to look into whole genomic repeats and other nongenic chromosomal regions that have never been explored before in conjunction with skeletal muscle aging. We observed that most retroelement copies, particularly those large-sized copies (>5 kb), were evidently demethylated in the skeletal muscles of superaged mice, and the clustering analysis uncovered two contrasting patterns of DNA methylation shifts within the related genomic repeats. We consider these global, age-related methylation changes to be one of the epigenomic characteristics of skeletal muscle aging.

| Statistics of MBD capture sequencing
We examined the epigenetic stability of skeletal muscle DNA during aging. The genomic DNA was extracted from the skeletal muscles of 2-month-old (young, 2 m), 20-month-old (aged, 20 m), and 28-month-old (superaged, 28 m; n = 3 each group) mice. With the genomic DNAs pooled by age, the MBD-based capture was performed to enrich methylated DNA fragments. Using the resulting DNA fragments, we performed high-throughput Illumina sequencing and obtained 62 million reads per age group on average, and the alignment efficiency rates ranged from 95.6% to 97.4% (Table S1).
The average GC content of the reads was ~65% in all groups (data not shown). This high GC content indicates the successful enrichment of methylated DNAs during MBD capture. The T-cell MBD-seq data, which were included for comparison, showed similar features.
T cells were chosen because they are implicated in chronic inflammation, which is known as the most important contributor to sarcopenia (Coelho Junior et al., 2016;Wilson, Jackson, Sapey, & Lord, 2017), a muscle disease (muscle failure) rooted in adverse muscle changes that accrue across a lifetime (Cruz-Jentoft et al., 2019).

| A perturbation in the methylomic structures of superaged skeletal muscle DNA
Box plots show the distributions of the read counts in 500-bp bins and reveal that there are increased mean methylation levels and decreased variance among the genomic regions in 28-m muscle ( Figure 1a). Histograms for the distribution of bins by fold-change values between the age groups showed a shift to hypermethylation in the 28-m muscle ( Figure S1). In the scatter plots comparing the methylation levels between the age groups, there was a perturbation in the methylomic structure as the slope of regression line was noticeably reduced in the matches of the 28-m muscle versus the 2-m or 20-m muscle (Figure 1b). This change was not evident in T cells. Figure 1c shows a representative genome browser view of a randomly chosen genomic region. Basically, both skeletal muscles and T cells possessed similar methylation profiles in all age groups, but notwithstanding the resemblance, the local events of the loss of methylation (LOM) and the gain of methylation (GOM) appeared only in the skeletal muscles of superaged mice.
To measure the extent of the methylation shift during aging, we examined the genomic regions with methylation levels at all values: top-ranked (high-level methylation; ~1 × 10 5 bins), bottom-ranked (low-level methylation; ~1 × 10 6 bins), or intermediate (~10 read counts; ~1 × 10 5 bins) after sorting the bins in 2-m muscles by their methylation levels. The high-methylation regions in the skeletal muscle showed an obvious LOM that was intensified with increasing age (p ~ 0, two-sample t test; Figure 1d). Conversely, both the low and moderate methylation regions showed GOM. These alterations accelerated with age; in the high-methylation regions, the rate

| Two contrasting patterns of methylation drift among the members of each repeat family
Genomic repeats are known to commonly undergo demethylation during aging (Pal & Tyler, 2016), which is in contrast to our result of the overall increase in DNA methylation in the 28-m skeletal muscle.
We examined the methylation states of various genomic repeats, including endogenous retroviruses (ERVs; ERV1, ERVK, ERVL, and MaLR) and long and short interspersed nuclear elements (LINEs and SINEs). We observed an age-dependent increase in the DNA methylation levels in all of these repeat families ( Figure S2). Clustering by the k-means algorithm (k = 2) classified the members in each repeat family into highly (cluster 1) and weakly (cluster 2) methylated copies (Figure 2). We found that the elements in cluster 2 of the 28-m muscle were hypermethylated. More interestingly, in all repeat families in the muscle, the two clusters displayed an opposite pattern of methylation changes: an LOM in cluster 1 and a GOM in cluster 2. This contrasting change was especially prominent in the ERV1, ERVK, LINE1, and SINE B1, B2, and B4 sequences. On the other hand, the ERVL and ERVL-MaLR sequences exhibited a GOM pattern only without a recognizable LOM, the significance of which is unknown at present. In contrast to the skeletal muscle tissues, the F I G U R E 1 A bimodal pattern of change in global DNA methylation states during murine skeletal muscle aging. (a) Box plots showing distributions of read counts (log 2 scale) in 500-bp bins obtained by MBD capture sequencing of skeletal muscle and T cells in 2-, 20-, and 28-month-old mice (2, 20, and 28 m, respectively; n = 3 each). (b) DNA methylome correlations between age groups. Each x-and y-axis indicates the read counts in log 2 scale. Coefficient of determination (R 2 ) and the equation for each linear regression line (blue) are denoted. (c) The browser view of DNA methylation of a randomly chosen genomic region (chr1:74,540,000-74,670,000) in mice of different ages. Arrows and brackets indicate the observed loci of demethylation (loss of methylation sites, LOM) and noisy methylation (gain of methylation sites, GOM), respectively, in the genome browser. Notably, CpG islands (CGIs) are devoid of methylation signals. (d) A bimodal pattern of methylation change during skeletal muscle aging. Fold changes on the log 2 scale of methylation levels of 20-m and 28-m muscles relative to 2-m muscle in top-ranked high-level methylation (n = 100,000), mid-ranked moderate-level methylation (n = 100,000; average read count ≅ 10), and bottom-ranked low-level methylation regions (n = 1,048,575) were calculated after sorting the genomic regions (bins) in 2-m muscle by methylation level. Error bars: standard deviation. nd, not determined; r, Pearson's correlation T cells stably maintained their methylation states in these genomic repeats regardless of age.

| A sharp demethylation of large retroelements in skeletal muscle
Since the retroelement copies in the ERV and LINE1 families are present in diverse sizes (ranging from tens to thousands of base pairs) across the genome, we interrogated the length distribution of the retroelements across the genome. The sizes of the fulllength ERVs (except ERVL-MaLR) and LINE1 vary and are approximately 5-7 kb on average (Crichton, Dunican, Maclennan, Meehan, & Adams, 2014). However, the majority of LINE1 and ERV elements in the genome are shorter than 1 kb (Figure 3a), indicating that the short copies were predominant in the methylation heatmap shown in Figure 2 and Figure S2. This finding prompted us to reanalyze those retroelements by size and to search for a size-dependent manifestation of demethylation in relation to the functionality of large retroelements. We subdivided the genomic copies of each retroelement family into three groups as follows: those >1, >3, and >5 kb. As illustrated in Figure 3b-f, unique methylation patterns were unearthed in the methylation heatmaps of large retroelement fractions, with alternating peaks and valleys of DNA methylation along the retroelement bodies, particularly around the start and end loci. Furthermore, we observed that these sharp methylation patterns became dull and blunt in 28-m skeletal muscle by LOM, unambiguously in cluster 1 of the ERV1, ERVK, and LINE1 families and moderately in ERVLs. In contrast, cluster 2 of ERVL and LINE1s showed GOM. ERVL-MaLR copies displayed GOM events only. Some LINE1 families including L1Md_A, L1Md_Gf, and L1Md_Tf are evolutionarily active and capable of independent mobilization with full-length transcripts (Hardies et al., 2000;Sookdeo, Hepp, McClure, & Boissinot, 2013). When we inspected age-related methylation dynamics in these three LINE1 families, we found that both the initial methylation level and the degree of LOM were far more prominent in these families than the rest LINE1s ( Figure 3g, see also Figure S3). It raises a possibility that the age-associated LINE1 demethylation can practically affect the genome stability in skeletal muscle through expression and transposition. We confirmed the age-related demethylation at these retroelement loci in an independent MBD-seq analysis ( Figure S4). T cells exhibited the same methylation patterns, which suggests the preservation of the methylation patterns of retroelement families across the cell types, although T cells did not show any recognizable changes in their DNA methylation in superaged mice ( Figure S5); these results agree with previous results that LINE1 methylation did not vary with age in human blood samples (Bollati et al., 2009;Talens et al., 2012).
F I G U R E 2 Age-associated changes in global DNA methylation in retroelements and single-copy sequences. Heatmaps of DNA methylation densities are shown within the start (S) and the end (E) sites of various retroelement families. Each group was subdivided by methylation level into two unsupervised clusters-cluster 1 (higher methylation) and cluster 2 (lower methylation)-and the mean methylation levels of the two clusters (blue and red lines, respectively) were examined separately, as shown in the graphs at the bottom. The numbers (n) of the members of each group are indicated. Schematic structures of ERVs and LINE1s are presented above The heatmap of DNA methylation demonstrated that the retroelements had distinct methylation patterns of their own: Methylation peaks were found to be (a) strong in the 5′ long terminal repeat (LTR) and weak at the 3′ LTR in ERV1s; (b) weak at the 5′ LTR and strong at the 3′ LTR in ERVLs; (c) strong at both ends in ERVKs with the 5′ peak usually the stronger; and (d) as single distinguished peaks in the 5′ untranslated region (UTR) in LINE1s. To test whether these unique methylation patterns indeed manifested themselves from individual copies of retroelements, we extracted only the large LINE1 (>5.0 kb) and ERV (>3.0 kb) copies from the re-peatMasker (mm10) track and uploaded them to the UCSC Genome Browser. The Genome Browser view in Figure 4a proved that the respective retroelements faithfully followed the specific patterns of their own families identified in the heatmap profiles with the peak heights lower in the 28-m muscles than the peak heights in the younger muscles. The T cells manifested the same methylation patterns across the retroelements, but an age-associated change in DNA methylation was not observed ( Figure S6). Meanwhile, it was unlikely that the methylation peak simply mirrored the distribution of CpG sites along the retroelements because the CpG dinucleotides were relatively evenly dispersed within each retroelement (ERVK, ERVL, and LINE1; Figure 4b). One exception was the ERV1 family, in which the CpGs were the most abundant and biased to the 5′ half region, and this CpG distribution resembled the peak positions in the methylation profile of ERV1s.

| Two different phases of retroelement expression during skeletal muscle aging
Since DNA demethylation at large retroelement sequences is related to the transcriptional derepression of the retroelements in superaged muscles, we investigated the genome-wide retroelement expression. To estimate the expression levels of the retroelements, RNA-seq was performed with the same batches of muscle tissues that were used for the methylation analysis (Table S1). The count distributions and the total numbers of genes expressed were similar among the samples, and principal component analysis separated the samples into groups of old and young skeletal muscles ( Figure   S7). To test the pertinence of the RNA-seq data, we examined Myh7, the gene of type 1 (slow-twitch) muscle fibers that is well known for its elevated expression in aged muscle (Ohlendieck, 2011); we also examined Tnni1, Myl12a, Myom1, Myh6, and Actn2, which are coexpressed with Myh7 (STRING; https ://string-db.org). Indeed, they were all significantly overrepresented in aged (20 m and 28 m) muscles ( Figure S8a). In contrast, the gene expression levels of Myh4, type 2B (fast-twitch) muscle fiber, and the coexpression genes, such as Tnni2 and Calm2, were not different among the age groups ( Figure   S8b). Myod1, which plays a major role in muscle differentiation (Rudnicki et al., 1993), was significantly underrepresented in old muscles ( Figure S8c), as was demonstrated in a previous study (Tamaki et al., 2000). Therefore, we concluded that our RNA-seq data matched We examined the expression levels of retroelements, particularly those from the ERV1, ERVK, and LINE1 families, which showed obvious demethylation in superaged muscles. The differential expression analysis suggested that the retroelement expression

| Expression patterns of the epi-driver genes that potentially affect the epigenetic states of retroelements
To find the cause of the age-associated fluctuation of retroelement expression, we analyzed the genes that are involved in the regulation of retroelements. Retroelement families are known to be repressed by the histone methyltransferase SETDB1 in cooperation with TRIM28 and HNRNPK (Matsui et al., 2010). As depicted in Figure 5f, F I G U R E 5 Changes in expression levels of retroelements during skeletal muscle aging. (a & b) MA plots for comparing the transcriptomes of whole retroelements between 2-m and 20-m muscles (a) or between 20-m and 28-m muscles (b). Transcriptome data were obtained through RNA-seq with the skeletal muscles of the same mouse groups (2-, 20-, and 28-month-old female mice; n = 3 each group) as used in the methylation analysis, and as a result, ~26 million mapped reads (including multiple mapped reads) on average per sample with 91% mapping efficiency were obtained. The expression levels of retroelements were determined by employing the "repeatMasker" track downloaded from the UCSC Genome Browser. (c & d) Scatter plots for comparing the expression levels of large retroelements between 2-m and 20-m muscles (c) or between 20-m and 28-m muscles (d). In panels (a-d), the numbers of differentially expressed retroelements are indicated by different colors (green for younger samples and red for older ones). (e) Age-associated changes in the transcript levels of various ERV1, ERVK, and LINE1 subfamilies. (f & g) Expression levels of genes involved in the repression of retroelements (f) and in the regulation of DNA methylation (g). Asterisks indicate significant differences (p < .05, Student's t test) from the levels in 2-m muscles. Error bars: standard deviation the mean expression levels of the genes Setdb1, Trim28, and Hnrnpk also fluctuated in a pattern that was well matched to that of retroelement expression because they increased in the 20-m muscles and then diminished in the 28-m muscles. Nevertheless, the extent of the change in the respective genes was not very significant, except for Hnrnpk; this indicates that these genes may not be the cause of the fluctuation of the retroelement expression levels (see Discussion).
Other genes, such as G9a-Glp (Maksakova et al., 2013), Suv39h1 (Bulut-Karslioglu et al., 2014), and Ezh2 (Leeb et al., 2010), whose products are also known to repress retroelements transcriptionally, were not expressed differently between the young and older muscles (data not shown). In addition, we examined the expression levels of genes that function in the establishment, maintenance, and modification/erasure of CpG methylation to test the possibility that the abrupt demethylation in the retroelements was associated with the perturbed expression in this group of genes. The RNA-seq data indicated that DNMT1 and DNMT3A are the DNA methyltransferase enzymes that primarily function in skeletal muscles, while DNMT3B and DNMT3L were barely expressed. Dnmt1 expression did not greatly change in aged muscles, whereas Dnmt3a expression significantly diminished (Figure 5g). The reduced DNMT3A level in the skeletal muscles may be the cause of global demethylation, but it cannot account for the observed methylation noise (i.e. a small and random increase in DNA methylation) across the genome. In addition, in the ten-eleven translocation (TET) protein family, which is the newest group of DNA methylation "editors" and retroelement regulators (de la Rica et al., 2016), Tet1 and Tet3 were significantly downregulated in the 28 m muscles; however, Tet2 was stably expressed. The connection of the reduced expression levels of the Tet genes with aging and DNA methylation changes in skeletal muscles is currently unknown (Tan & Shi, 2012).

| D ISCUSS I ON
We found that both the GOM (hypermethylation) and LOM (hypomethylation) events coexisted among the related members of each repeat family (ERV1, ERVK, LINE1, and SINE B1, B2, and B4) in superaged skeletal muscles (Figures 2 and 3). The discovery of the two contrasting patterns of methylation changes was attributed to the use of analytical tools, such as methylation-based clustering and the size-based fractionation of retroelement copies; this was believed because the analysis results with whole retroelement copies of each family produced just a simple one-way change ( Figure S2). This bidirectional pattern of methylation drift indicates that the age-associated methylomic change may be driven by composite, antithetical forces instead of by a unidirectional flow of event (i.e., hypomethylation) that has been generally thought to sweep over the genomic repeats in aged cells. Following this "regression toward the mean" tendency of methylation drift, aging may tend to progressively reduce the preexisting methylation differences between genomic regions, the pattern of which had been established during muscle development and differentiation.
The gradual expansion of a disorder in the epigenomic structure leads to a perturbation of transcriptomic homeostasis and a steep rise of transcriptional noise; the latter is another landmark of aging and is a primary cause of variability in the gene expression between cells in an isogenic population (Enge et al., 2017). Supposing that DNA methylation is the primary mechanism regulating the transcription of genes in the corresponding regions, the LOM events in regions with heavy methylation may make the promoters more accessible to transcriptional activators, whereas the GOM events in the hypomethylated regions may make the promoters more and more occupied by repressors. This condition most likely results in a bimodal transcriptional drift, where the genes that are weakly expressed in the young animals are increasingly transcribed in older animals; in contrast, the genes that are highly expressed in the young animals are gradually downregulated with age (Min, Park, Jeon, et al., 2018). Our observation of age-associated epigenetic drift therefore lends sound support to the bimodal transcriptional drift theory and is consistent with the general view that aging reflects a stochastic process of increasing disorder. Meanwhile, some studies observed only a weak correlation between DNA methylation (particularly at the promoters) and gene expression changes (Marttila et al., 2015;Slieker, Relton, Gaunt, Slagboom, & Heijmans, 2018;Steegenga et al., 2014). Notably, these studies used the Infinium 450K platform, which enables the identification of DNA methylation changes with a single CpG resolution (Bibikova et al., 2011). The base-level sensitivity notwithstanding, it is unknown how well the selected few CpGs represent the surrounding cluster(s) of CpGs for the overall methylation state and whether the selection is grounded with an experimental basis (i.e., the connection with gene expression). Therefore, there is some difficulty directly linking the methylation of the selected CpGs to the expression of associated genes, and we assume that the difference between their studies and the present study lies in our approach to such a whole-genome methylation analysis as MBDseq. As a drawback of the MBD-seq method, it is semiquantitative in the sense that it does not yield direct estimates of CpG methylation levels. Thus, the previously suggested pipeline of methylome study (Aberg et al., 2012) looks relevant: MBD-seq offers a first pass at global methylation analysis, and then, the identified genomic loci of interest are subject to a base-level validation using technologies, such as pyrosequencing and targeted bisulfite PCR sequencing (TBP-seq; . DNA methylation contributes to the transcriptional repression of retroelements, including the intracisternal A particles (IAPs) of the ERVK family. For example, IAP transcript levels are elevated 50-to 100-fold in murine Dnmt1 knockout embryos, with a parallel demethylation of the IAP promoter (Walsh, Chaillet, & Bestor, 1998). Intracisternal A particles derepression occurs in stem cells of the male germline at undifferentiated stages when genome-wide demethylation occurs. In addition, the induction of IAP expression occurs after treatment with demethylating agents, such as 5-azacytidine (Hojman-Montes de Oca et al., 1984). Retroelements, including Alu and LINE1s, were subject to derepression in cultured cells during senescence (De Cecco et al., 2013). As a possible mechanism for the derepression of retroelements, the authors of that study suggested the acquisition of a permissive chromatin state in senescent cells, which is likely related to a demethylation-induced perturbation of the heterochromatin structure at retroelement loci (Sedivy, Banumathy, & Adams, 2008), as we have demonstrated in this study.
Based on these precedent findings, the derepression of various retroelements could be predictable in superaged skeletal muscle, as we observed demethylation at various retroelements (Figure 3). In reality, the global demethylation event temporally coincides with an increase in the expression of retroelements, including the IAP sequences in skeletal muscles, as the mice advance in age from 20 to 28 months ( Figure 5). Recent papers show that LINE1 expressions increase with age in mouse skeletal muscle (De Cecco et al., 2019;Simon et al., 2019), and our result on the age-associated demethylation of LINE1, especially of those evolutionarily active LINE1 families (Figure 3g), now provides a mechanistic link with these observed increase in LINE1 expression. Meanwhile, the accumulation of retroelement-derived transcripts and their reverse-transcribed products in skeletal muscles may trigger an innate immune response and chronic inflammation through inducing IFN-α and IFN-β gene expressions, as suggested before (Bourque et al., 2018;De Cecco et al., 2019;Kang, 2018;Simon et al., 2019).
Various retroelements were more repressed in the 20-m skeletal muscle than they were in the 2-m muscle; then, these retroelements became derepressed in superaged muscle (Figure 5e). The tighter repression of retroelements in the 20-m skeletal muscle appears to be independent of DNA methylation because there was no apparent methylation change in the retroelements during that time. It would be interesting to determine by what mechanism the retroelements become even further repressed during early aging and what the significance/purpose is in it. In addition to DNA methylation, histone methylation contributes to the repression of genomic retroelements.
Retroelements are controlled by histone H3 lysine 9 methylation (H3K9me) catalyzed by SETDB1 [see (Kang, 2018) and references therein], G9A/GLP (Maksakova et al., 2013), and SUV39H1 and SUV39H2 (Bulut-Karslioglu et al., 2014), and/or by H3 lysine 27 methylation (H3K27me) performed by the PRC2 complex involving EZH2 (Leeb et al., 2010). However, a recent study (Min, Park, Jeon, et al., 2018) has suggested that the expression levels of epi-driver genes encoding these histone-modifying enzymes and other related epigenetic modifiers were relatively stable and had little change during skeletal muscle aging. One seemingly tiny, but nonetheless noticeable, change was that the expression levels of the genes encoding members of the well-known KRAB-ZNF-TRIM28-SETDB1-HNRNPK complex, which represses retroelements, all consistently manifested a negative correlation with the retroelement expression levels (Figure 5f). The same was true for TET2, which was recently reported to silence retroelements by the post-transcriptional destabilization of their transcripts (Guallar et al., 2018). In this regard, we speculate that although the changes in the expression of the respective genes cannot lead to a significant outcome, the sum of their small changes can make a visible difference in the transcript levels of retroelements. Skeletal muscles are made of different types of muscle fibers, and the composition changes with age (Ohlendieck, 2011). Given that one of the major drivers of DNA methylation pattern is cell type, it could be that the altered methylation pattern in the superaged muscle reflects a change in the composition of muscle fiber types with age. We only know that the human type 1 slow-twitch and type 2a fast-twitch muscle fibers had similar levels and distributions of methylated CpGs at promoters and CGIs (Begue, Raue, Jemiolo, & Trappe, 2017). However, there is no information about whether the different muscle fibers undergo unique methylation changes in their own pattern with age. In this study, we noticed that the type 1 slow-twitch muscle fiber genes, such as Myh7, Tnni1,Myl12a,Myom1,Myh6, and Actn2, which are expressed at higher levels in aged muscle (Ohlendieck, 2011), were increasingly expressed in the 20-m and 28-m muscles compared to their expression in the 2-m muscles (seen in Figure S8). However, given no further difference between the 20-m and 28-m muscles, we hypothesize that the muscle fiber composition is not very different between the 20-m and 28-m muscles and, accordingly, that the distinct methylation pattern in superaged muscle is ascribable to the aging process itself, rather than to the age-associated change in muscle composition. Nevertheless, we cannot completely exclude the possibility of age-specific compositional changes in other types of muscle cells that can make differences in the DNA methylation patterns in aged skeletal muscles.
If this were the case, we could tell that the regression to the mean tendency in skeletal muscles is rooted in the presence of the variety of cell types and that the lack of such a tendency in the MAX-purified T cells appears to make sense.
The degree of DNA methylation in retroelements can affect the expression levels of the associated genes. A well-known example is the agouti viable yellow (A vy ) allele. The A vy allele contains an IAP insertion close to the transcription start site of the agouti gene, and its expression level inversely correlates with the CpG methylation state of the IAP sequence (Michaud et al., 1994). In Arabidopsis thaliana, the methylation status of intronic transposable elements was observed to correlate with lower transcription levels of the associated genes with TE insertions (Le, Miyazaki, Takuno, & Saze, 2015). In this study, we observed two categories of large retroelements: the heavily methylated copies being demethylated with aging and the other, poorly methylated copies (cluster 2), which constitute approximately 10%-30% in each retroelement family and naturally do not undergo a methylation change (Figure 3b-3f). Thus, those genes that are intimately associated with the former group of retroelements may be transcriptionally modulated after the demethylation of retroelements during aging. If so, these genes are considered age-sensitive and can serve as marker genes for skeletal muscle aging.

| Isolation of splenic T cells and skeletal muscles
To obtain skeletal muscles and T cells, C57BL/6J female mice at 2 months (n = 3), 20 months (n = 3), and 28 months (n = 3) of ages were sacrificed, and their spleens and skeletal muscles in hindlimbs were surgically removed. Mice about 28 months of age, which correlates with humans ranging from 76 to 78 years of age, can be considered "very old" and generally show a 50% survivorship (The Jackson Laboratory; lifespan as a biomarker). For isolation of T cells, a spleen from each mouse was ground with frosted glass slides followed by a series of syringe homogenization as shrinking the needle sizes (from 18 to 27 gauges), and T cells were isolated by magnetic-activated cell sorting (MACS; Pan T Cell Isolation Kit II, mouse; Miltenyi Biotec) as described elsewhere . For skeletal muscles, the whole muscle tissues were completely frozen in liquid nitrogen, and they were ground to powder using a mortar and a pestle. The powdered tissues were further homogenized using a Biomasher II (DWK Life Sciences) in tissue lysis buffer (ATL buffer; Qiagen).

| Construction of MBD-seq libraries
Genomic DNAs (gDNAs) were extracted from the tissue samples using DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer's instruction. Genomic DNAs obtained from three mice of the same ages were pooled before sonication using Bioruptor Pico (Diagenode) and purified using with AMPure XP beads (Beckman) . For preparation of MBD capture, 7 μl MBD-Biotin proteins per sample from MethylMiner Methylated DNA Enrichment Kit (Thermo) was first incubated with 10 μl of precleared Dynabeads (Thermo) at RT for 1 hr. Methylated DNAs were captured by incubating 1 μg of sonicated gDNA fragments with MBD beads at 4°C for overnight. Beads were washed twice, and captured DNA fragments were eluted using 2 M NaCl and precipitated using ethanol. Enrichment of methylated DNA fragments was checked by PCR for control sets of methylated and nonmethylated DNAs (Thermo) that were premixed as spike-ins. Using the enriched methylated DNA fragments, Illumina sequencing libraries were generated. All the captured DNAs were end-repaired by reacting with T4 DNA polymerase (NEB), Klenow (NEB), and T4 polynucleotide kinase (NEB) in the presence of dNTP at 20°C for 2 hr, and then 3′-end adenylated using Klenow fragment (NEB) at 30°C for 2 hr.
Finally, the end-repaired DNA fragments were ligated with NEBNext adapters (NEB) and were enriched by 18 cycles of PCR. The resulting libraries were subjected to Illumina sequencing with NEXTseq 500. To confirm the MBD-seq result, we repeated the same experiment of MBD capture and sequencing using the same set of pooled gDNAs as in the first-round MBD-seq experiment.
For analysis of RNA-seq results, the trimmed reads were aligned on the reference genome (mm10) using Hisat2 with "--dta" option for splice junction mapping. For estimation of gene/retroelement expression, "refSeq" and "repeatMasker" tracks, respectively, were downloaded from UCSC genome browser in GTF, and based on the genomic coordination in GTF files, reads were counted using "htseq-count" with "intersection-strict" parameter. All count data were merged into one data table using a home-brew bash script, and the raw count data were normalized and used for differential expression analysis with "DESeq2." Genes or retroelements with fold change >2 and p-value <.05 were considered differential expression. Bam files were converted to bigWig files for visualization on UCSC genome browser using "bamCoverage" in "deepTools."

| Validation-retroelement expression
Total RNAs were extracted from the muscle tissue of each age group as described above. Reverse transcription was performed by incubating 1 μg of DNase I-pretreated RNA with Superscript III enzyme (Invitrogen), 20 μM oligo-dT primers (Invitrogen), and 50 ng random hexamers (Invitrogen) at 50°C for 1 hr. Ten ng of the synthesized cDNAs was used for real-time quantitative PCR (QuantStudio3 Real-Time PCR system, ABI) with the specific primers for individual retroelement subfamilies (Table S1). The PCR was performed with a following program, 10 min of predenaturation at 95°C followed by 40 cycles of 95°C for 15 s, 55°C for 15 s, and 72°C for 1 m. Finally, relative expression levels of each retroelement against Gapdh were calculated using QuantStudio Design & Analysis Software (Thermo).

ACK N OWLED G EM ENTS
We thank J.H. Lee and M. Park for technical assistance.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

AUTH O R S ' CO NTR I B UTI O N
BM conducted bioinformatic analyses and constructed MBD-seq libraries. KJ constructed RNA-seq libraries. JSP maintained and provided mice. YKK designed and supervised the study, performed bioinformatic analyses, and wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Gene Expression Omnibus (GEO) with accession number of GSE12 3954.