Locus-specific H3K9me3 gain in aged somatic tissues in Caenorhabditis elegans

Epigenetic alterations occur as organisms age, and lead to chromatin deterioration, loss of transcriptional silencing and genomic instability. Dysregulated epigenome has been linked to increased susceptibility to age-related disorders. We aim to characterize the age-dependent changes of the epigenome and, in turn, to understand epigenetic processes that drive aging phenotypes. In this study, we focused on the aging-associated changes in the repressive histone marks H3K9me3 and H3K27me3 in C. elegans. We observed redistribution of of both histone marks, but the changes are more significant for H3K9me3. We further found alteration of heterochromatic boundaries in aged somatic tissues. Interestingly, we discovered that the most significant changes reflected H3K9me3-marked regions that are formed during aging, and are absent in developing worms, which we termed “aging-associated repressive domains” (AARDs). These AARDs preferentially occur in genic regions that are marked by high levels of H3K9me2 and H3K36me2 in larval stages. Interestingly, maintenance of high H3K9me2 levels in these regions have been shown to correlate with longer lifespan. Next, we examined whether the changes in repressive histone marks lead to de-silencing of repetitive DNA elements, as reported for several other organisms. We observed increased expression of active repetitive DNA elements but not global re-activation of silent repeats in old worms, likely due to the distributed nature of repetitive elements in the C. elegans genome. Intriguingly, CELE45, a putative short interspersed nuclear elements (SINE), was greatly overexpressed at old age and upon heat stress. SINEs have been suggested to regulate transcription in response to various cellular stresses in mammals, it is likely that CELE45 RNAs also play roles in stress response and aging in C. elegans. Taken together, our study revealed significant and specific age-dependent changes in repressive histone modifications and repetitive elements, providing important insights into aging biology. Author summary Heterochromatin refers to the portion of the genome that is tightly packed where genes stay silent. Heterochromatin is typically decorated by particular chemical groups called histone modifications, such as trimethylation of lysine 9 or lysine 27 on histone 3 (H3K9me3 or H3K27me3). To understand how the heterochromatin landscape may change from a “youthful” to an “aged” state, we monitored the genome-wide patterns of H3K9me3 and H3K27me3 during aging using the genetic model soil worm C. elegans. We found that while H3K27me3 remained relatively stable with age, H3K9me3 showed profound genome-wide redistribution in aged worms. We observed that new H3K9me3-marked heterochromatin preferentially formed in specific gene-rich regions in aged worms. Interestingly, these particular regions were marked by high levels of three other histone modifications when worms were young. This result suggested that H3K9me3 gain during aging is influenced by the gene-specific landscape of histone modifications established at young age rather than occurs in a stochastic manner. In summary, our study discovered reproducible and gene-specific changes in histone modifications that likely contribute to the aging phenotypes.

5 during aging in whole worms and human cells [14,15,18,23,24]. In fly, the levels of H3K9me3 93 were found to decrease in the aging intestine [25] but to increase in the aging brain [26]. In 94 mouse, the levels of H3K27me3 reduced in the senescent fibroblast cells [23] but elevated in 95 the brain of a mouse model of accelerated aging [27]. Therefore, it remains unclear how age-96 associated changes in H3K9me3 and H3K27me3 levels contribute to heterochromatin 97 deterioration and aging phenotypes. It has been proposed that it is not the net abundance, but 98 rather the genomic distribution and utilization of the histone modifications that regulate the aging 99 processes [28].

100
Caenorhabditis elegans lacks cytologically dense-staining heterochromatin due to its 101 holocentric chromosomes with dispersed centromeres. Repressive H3K9me3 and H3K27me3 in 102 worms are predominantly localized on repeat-enriched distal arms of all five autosomes, where 103 their distributions are largely overlapped [29][30][31]. In contrast, heterochromatin is mostly marked 104 by either H3K9me3 or H3K27me3 in monocentric species, such as fly and human [31].

105
Moreover, H3K27me3-marked regions that lack H3K9me3 can be found in the central regions of 106 autosomes and on the X chromosome in worms. In this study, we investigated how H3K9me3 107 and H3K27me3 change during aging in the somatic tissues of germlineless glp-1 mutants by 108 chromatin immunoprecipitation (ChIP) followed by deep sequencing. We found site-specific loss 109 and gain of H3K9me3 and H3K27me3 in the aged somatic cells. Age-dependent loss of 110 repressive marks was predominantly found in the central regions of autosomes while gain was 111 localized on the distal autosomal arms. Furthermore, the redistribution of H3K9me3 at old age 112 resulted in the formation of new H3K9me3-marked heterochromatic regions. These H3K9me3 113 peak regions that newly formed at old age were preferentially localized at genic regions marked 114 by H3K9me2 and H3K36me1/2 at the juvenile stages. We hypothesized that the H3K9me3 gain 115 at these regions was a result of H3K9me2-to-H3K9me3 conversion and that the trimethylation of 116 H3K9me2 could contribute to structural and functional changes of chromatin at old age. We 117 6 next examined whether somatic aging in the glp-1 mutants is associated with loss of 118 transcriptional suppression of repetitive DNA elements. Although we did find an increase in the 119 levels of RNAs derived from repetitive DNA elements at old age, that nevertheless represented 120 a small fraction (<1%) of the overall transcriptome. Interestingly, we found CELE45, a short 121 interspersed nuclear element (SINE), to be significantly over-represented among the 122 upregulated repeats in aged worms. In murine and human cells, SINE RNAs were induced and 123 implicated in cellular stress response [32][33][34]. To this end, we re-analyzed the publicly available 124 RNA-seq data and found evidence supported that CELE45 was induced upon heat stress in 125 worms. Moreover, overexpression of CELE45 in the aged germlineless glp-1 mutants was not 126 due to the culture temperature at 25°C and/or lack of germline. These results indicated that 127 CELE45 was specifically induced in response to aging.

130
Redistribution of H3K9me3 and H3K27me3 in aged somatic tissues in C. elegans

131
To examine the genome-wide patterns of repressive histone modifications (H3K9me3 132 and H3K27me3) in young and aged somatic tissues of C. elegans, we performed chromatin-133 immunoprecipitation followed by sequencing (ChIP-seq) in germlineless worms (S1 Table). We 134 used glp-1(e2141) mutants that lack germ cells at the non-permissive temperature (25°C) to 135 harvest post-mitotic somatic tissues and prepare whole worm extracts [35]. For both repressive 136 histone marks, we had two ChIP-seq biological replicates in each of the two time points (day 2 137 and day 12 of adulthood). We chose day 2 as the young time point when wild-type worms 138 normally reach peak reproduction, and day 12 as the old time point when 10-20% of the glp-139 1(e2141) population would have died [36]. To compare the similarity of the ChIP-seq replicates,

140
we performed pair-wise Pearson's correlation analysis on H3, H3K9me3 and H3K27me3 ChIP-141 seq tag coverage in 15kb sliding windows along the genome (S1 Fig). Biological replicates for 7 the two repressive marks and H3 control at the same age had genome-wide correlation 143 coefficients higher than 0.83, indicating high reproducibility.

145
Next, we identified H3K9me3-and H3K27me3-enriched peak regions normalized to the 146 control H3, using the broad peak parameters in MACS2 (v2.1.0). We combined all biological 147 replicates for peak calling and merged peaks identified at both time points (day 2 and day 12 148 adults) in order to generate common regions for downstream age-dependent analyses. We  Table), and 5,870 H3K27me3 peak regions (S3 Table). The

155
H3K9me3 peak regions were predominantly enriched in the repeat-rich distal arms of 156 autosomes (S2 Table and

171
To examine whether there were age-dependent changes in H3K9me3 and H3K27me3 172 enrichment at ChIP-seq peak regions in germlineless glp-1 animals, we first visualized the 173 similarity between young and old samples by performing multidimensional scaling (MDS) 174 analysis. The MDS plots showed that H3K9me3 ChIP-seq experiments were separated and 175 clustered by age, indicating there were substantial age-dependent differences in H3K9me3 176 enrichment ( Fig 1A). In comparison, H3K27me3 experiments in young and old adults showed 177 relatively minor age-dependent variations ( Fig 1B). Histone 3 ChIP-seq profiles were highly 178 similar between the two ages ( Fig 1A and 1B

181
Next, we examined how the average H3K9me3 and H3K27me3 signals at ChIP-seq 182 peak regions changed between young and old glp-1 mutants. We included our published 183 H3K4me3 and H3K36me3 ChIP-seq data in the analysis for comparisons [36,37]. At H3K9me3 184 peak regions, we found an overall reduced H3K9me3 enrichment at old age and concomitant 185 changes for the active histone marks H3K4me3 and H3K36me3 at the boundaries. In young 186 somatic tissues, H3K9me3 peak regions coincided with a sharp fall in the levels of these active 187 histone marks (H3K4me3 and H3K36me3), exhibiting clear boundaries ( Fig 1C). This distinctive 188 feature became blurred in aged somatic tissues, indicating invasion of active marks in 189 H3K9me3-marked constitutive heterochromatin during aging. At H3K27me3 peak regions, we 190 also observed decreased H3K27me3 enrichment but the more gradual decline of active histone 191 marks at the boundaries remained stable with age ( Fig 1D). It was known that the mutually 9 exclusive patterns of H3K27me3 and H3K36me3 were maintained by the antagonistic 193 relationship between Polycomb repressive complex 2 (MES-2/3/6) and H3K36 194 methyltransferase (MES-4). We found that the RNA expression of the four genes remained 195 stable during aging [36], which would be consistent with our observation that the H3K27me3 196 and H3K36me3 boundaries were preserved in aged glp-1 germlineless mutants.

198
Like other eukaryotes, the C. elegans genome is organized into transcriptionally active 199 and inactive chromatin compartments in the nucleus [40,41]. We next examined whether The H3K9me3 and H3K27me3 repressive histone marks are generally associated with 216 transcriptional repression. To confirm this, we used our previously published RNA-seq data of 217 10 day 2 and day 12 glp-1(e2141) germlineless adults [36] to assess the expression of the genes 218 whose coding sequences overlap with either H3K9me3 or H3K27me3 peak regions. As 219 expected, the majority of these genes showed no detectable RNA expression at both the young 220 and old time points (64.2% and 73.6%, respectively). Moreover, actively expressed genes 221 associated with the repressive peak regions showed significantly lower RNA levels compared to

226
To identify the peak regions with significant changes in H3K9me3 and H3K27me3 with 227 age in the glp-1 mutants, we performed differential analysis using edgeR. We found that 228 H3K9me3 peak regions had more statistically significant changes with age (595 peaks; 4.7% of 229 all peaks; S4 Table and  Next, we examined whether age-dependent changes in the histone marks correlate with 243 gene expression changes in the glp-1 mutants. As a control, we first analyzed our published 244 H3K4me3 peaks and RNA-seq data in the young and old glp-1 mutants [36]. The genes that 245 overlapped with significantly upregulated and downregulated H3K4me3 peaks showed net 246 positive and negative changes, respectively, in RNA expression levels comparing to genes in 247 non-peak regions (S6C, S6D and S6E Fig), consistent with our previous findings (2). Next, we 248 examined the expression of genes that overlapped with differential H3K9me3 and H3K27me3 249 peak regions. For genes associated with significantly upregulated H3K9me3 or H3K27me3 peak 250 regions, we found significant decrease in RNA expression changes at old age (Fig 2A and 2B).

251
On the contrary, genes overlapped with significantly downregulated H3K9me3 or H27me3 peak 252 regions showed significantly elevated RNA expression changes (Fig 2C and 2D). We then 253 examined the genes that were differentially expressed at old age (FDR < 0.01) and found a

257
It is however important to note that H3K9me3, and H3K27me3, might influence the expression 258 of genes whose coding sequences do not directly overlap with their peak regions, which would 259 be difficult to assess based on the available data.

261
New H3K9me3-marked repressive domains were formed in aged somatic tissues

262
To gain more insights into the genes associated with differential repressive peak 263 regions, we used WormCat to identify their enriched gene ontology (GO) terms [42] and 264 WormBase enrichment analysis tool to identify their tissue specificity and phenotype clusters 265 [43]. We found that there was no significant enrichment in tissue types and phenotypes for 266 genes associated with downregulated H3K9me3 and H3K27me3 peak regions. Hereafter, we 12 focused on the genes associated with the upregulated H3K9me3 and H3K27me3 peak regions.

268
In GO analysis, upregulated H3K9me3 peaks were associated with genes overrepresented for 269 the functional terms related to fundamental cellular processes, particularly "chromatin" (S6F

275
While the majority of genes associated with repressive H3K9me3 and H3K27me3 peak 276 regions were silent, it was puzzling that significantly upregulated repressive peak regions were 277 preferentially associated with actively expressed genes. One possibility was that the regions 278 that showed statistically significant upregulation of H3K9me3 and H3K27me3 in aged worms 279 had very low levels of the repressive marks in young worms, so the genes in those regions were 280 actively expressed at the young stage. Indeed, the average plots showed that the upregulated 281 peak regions were deprived of H3K9me3 and H3K27me3 marks at young age (Fig 3A and 3B).

282
At old age, these regions accumulated H3K9me3 and H3K27me3 at levels comparable to the 283 average of the non-differential peak regions (Fig 3A and 3B). This result indicated that certain ChIP-seq and ATAC-seq data. We included data for 26 histone marks, 268 293 transcription/chromatin factors, and one ATAC-seq analysis [31,36,37,[44][45][46][47][48][49]. We first used 294 correlation analysis to detect whether two peak profiles distributed and clustered in similar 295 regions in the genome. We compared the genome-wide distribution of the H3K9me3 or 296 H3K27me3 peaks that changed significantly with age with the peak profile of each of the 295 297 available datasets, by calculating the peak coverage in 50kb sliding windows across the entire 298 genome and performing pair-wise Person's correlation analysis (S7 Fig and S6 Table). Next, we 299 used permutation test to determine if two sets of peaks have more frequent overlap than 300 expected, and computed pair-wise z-scores (S4 Fig and S7 Table). We looked for top-ranked 301 peak profiles in both tests when compared with differential H3K9me3 and H3K27me3 peaks 302 during organismal aging. Among all the comparisons, we found that the peak regions with 303 significant gain in H3K9me3 strongly correlated with the peak profiles of H3K9me2 and 304 H3K36me1/2 in wild-type L3 larvae and mixed-stage populations (Pearson's correlation 305 coefficient > 0.45 and permutation z-score > 34). As expected, two H3K9me2-binding factors 306 (HPL-2, LIN-61) and the H3K9me1/2 methyltransferase (MET-2) also had high correlation 307 scores (Pearson's correlation coefficient > 0.39 and permutation z-score > 23). The remaining 308 pair-wise correlations were relatively weak in at least one of the two analysis and would not be 309 further discussed here. We also did not uncover any strong correlation for the differential 310 H3K27me3 peaks. It is important to note that the majority of the publicly available data were 311 from wild-type L3 larvae or mixed-stage embryos. Taken together, the results suggested that

354
only 19 of the 2,867 active repeats were silenced at young age, indicating repeat expression at 355 old age was not due to age-dependent derepression. Although the relative expression of repeat 356 elements increased (159%) with age ( Fig 4A; S8 Table), the overall abundance of repeat RNAs 357 remained relatively low (0.3~0.5% of the transcriptome) in the glp-1 worms (Fig 4A). Using 358 edgeR analysis, we identified 173 (0.28%) and 75 (0.12%) significantly upregulated and 359 downregulated repeats, respectively (S9A Fig and S8 Table). Taken together, these results to be distributed on all the chromosomes (Fig 4D and S9C Fig). Strikingly, the CELE45 family, in the glp-1 mutants. Since CELE45 is related to the SINE family of retrotransposons, which is 392 known to be transcribed by RNA polymerase III (pol III), we next examined whether the actively 393 expressed CELE45 copies were bound by RNA pol III in C. elegans, using the published ChIP-394 seq data of the RNA polymerase III subunit RPC-1 in wild-type L3 stage (modENCODE 6300).

395
We found most of the upregulated CELE45 copies (95%) were bound by Pol III whereas about 396 one-third of the silenced copies (36%) showed Pol III binding (S10E Fig  We were intrigued by the age-dependent overexpression of the CELE45 in aged somatic 411 tissues in glp-1 mutants. CELE45 is a putative short interspersed nuclear elements (SINEs). In 412 mammals, SINEs are rapidly induced in response to a variety of cellular stresses, including heat 413 shock and virus infection [32][33][34]. To this end, we examined whether CELE45 is induced after 414 heat shock in two previously published RNA-seq data in wild-type C. elegans [52,53]. Using the 415 same mapping and differential analyses strategies we applied to our glp-1 aging data, we found 416 18 that a large number of CELE45 copies were significantly upregulated after heat stress in the two 417 independent published data sets [52,53]. In one data set, we found 59 CELE45 copies were 418 rapidly upregulated after short-term heat shock (30 minutes at 33°C; S8

424
In this study, we used germlineless glp-1(e2141) worms raised at the non-permissive 425 temperature (25°C). We asked whether CELE45 induction was a result of lacking germline 426 and/or being cultured at 25C. To this end, we compared the locus-specific expression of 427 repetitive elements between wild-type worms raised at 20°C and glp-1(e2141) raised at 25°C 428 using published RNA-seq data [54]. We found that the expression of repetitive elements in 429 young adults at both conditions was surprisingly similar with very few differentially expressed 430 repeats (12 repeat loci in S8 Table), indicating that the lack of germline or a culturing 431 temperature between 20C and 25C did not significantly impact repeat expression. This 432 observation supported the notion that CELE45 upregulation was a specific response to aging.

435
In this study, we had investigated changes of H3K9me3 and H3K27me3 in aged somatic 436 tissues of C. elegans by ChIP-seq. We found that with aging, the overall H3K9me3 levels 437 increased at heterochromatic regions in the distal arms of chromosomes and decreased in 438 euchromatic central regions of autosomes (Fig 1E). H3K27me3 peak regions showed similar 439 patterns of age-dependent changes but at a lower magnitude (Fig 1E). Our data also revealed 440 specific regions with statistically significant age-dependent gains of H3K9me3, which resided in genomic regions not marked by H3K9me3 during the juvenile and young adult stage (Fig 3A), 442 thus representing aging-specific repressive regions (hereafter referred to as ASRRs). Newly 443 formed ASRRs at old age were enriched for H3K9me3 mark but not K3K27me3, whereas the 444 majority of heterochromatic regions formed during C. elegans development were dually marked 445 by H3K9me3 and H3K27me3. It is important to note that we used whole worm extracts of 446 germlineness glp-1 mutants for this study, and therefore our data could not distinguish whether 447 changes in H3K9me3 and H3K27me3 occurred in all or a specific subset of aged somatic 448 tissues. Nevertheless, the reproducible and robust gain of H3K9me3 in ASRRs detected by 449 ChIP-seq at old age suggested that it likely occurs in a substantial fraction of the aged somatic 450 tissues. Future studies investigating tissue-specific H3K9me3 profiles through a developmental 451 and aging time course will help to address this question.

452
Most healthy eukaryotic cells have either round-or oval-shaped nuclei. Abnormal 453 nuclear morphology, which can induce gene expression changes (1,2), occurs during natural 454 and premature aging in worms, flies, and human (3-5). Many previous studies have linked 455 depletion or defects of lamins, the major structural proteins at the nuclear periphery, to the 456 changes in nuclear morphology (6). However, a recent study provided evidence that altered 457 repressive histone modifications are sufficient to cause nuclear blebbing without detectable 458 changes in lamins (7). Therefore, it will be important to investigate whether the redistribution of 459 repressive marks observed in this study is linked to nuclear abnormalities during aging in future 460 studies.

461
Our analysis revealed that ASRRs occur in genomic regions marked by high levels of H3K9me2 462 and H3K36me2 at juvenile stages. H3K9me2 is a histone mark co-localized with 463 heterochromatin tethered to the nuclear lamina (8), whereas H3K36me2 is an active histone 464 mark associated with actively transcribed genes (9). H3K9me2-marked heterochromatin has 465 been reported to be transcriptionally permissive in yeast, worms and flies (10-13). Here, we 466 20 found that ASRRs preferentially overlap with gene bodies (92%, 529 of 578 ASRRs) and the 467 majority of the overlapped genes are transcriptionally active in somatic tissues (14). Importantly, 468 our data indicated that the H3K9me3 gain in ASRRs at old age results in transcriptional 469 repression (Fig 2B and 2F). Taken together, the reproducible aging-associated H3K9me3 gain   Table). In human and mouse, 508 SINE RNAs are upregulated by cellular stresses and implicated in regulating the expression of 509 stress-responsive genes (27-32). Although it is still unclear whether CELE45 participates in 510 stress responses, we found that CELE45 was induced after heat shock in the published RNA-511 seq data (33-35). Therefore, CELE45 over-expression at old age is likely a regulated process in 512 response to aging-associated cellular stresses. Going forward, it will be very interesting to

540
Two independent ChIP-seq experiments were performed. The first experiment (exp1) 541 generated two biological replicates of H3 and H3K27me3 ChIP-seq. The second (exp2) 542 generated two biological replicates of H3 and H3K9me3 ChIP-seq.

560
To compare the similarity between biological replicates, the number of mapped reads 561 were counted in 15kb sliding windows. Pearson's correlation coefficient was calculated between 562 the coverage tracks of individual replicates. To visualize ChIP-seq enrichment profiles in young 563 and old animals, mapped reads were first converted into bedGraph tracks of the log2 signal 564 ratio of treatments (H3K9me3 or H3K27me3) over the controls (H3) using deeptools 565 24 (bamCompare --bamfile1 H3K9me3orH3K27me3.bam --bamfile2 H3.bam --binSize 20 --566 operation log2 --scaleFactorsMethod readCount --smoothLength 60 --extendReads 200 -of 567 bedgraph). Next, the log2 signal ratio in bedGraph tracks was transformed into z-scores. Z-568 score bedGraph tracks were converted into bigWig tracks using UCSC bedGraphToBigWig 569 program. computeMatrix utility in deeptools was used to calculate a matrix of z-scores from the 570 bigWig tracks for generating metaplots and heatmaps.

572
Correlation analysis of peak profiles 573 To compare the similarity between peak profiles, we compared their genome-wide 574 distribution and frequency of overlap. We evaluated the genome-wide distribution by calculating 575 the peak area (bp) in 50kb sliding windows across the genome of individual peak profiles. Then,

576
we computed the pair-wise Pearson's correlation coefficients between peak coverage matrices.

577
To determine whether the pair-wise overlap frequency is significantly higher than expected by 578 chance, we used regioneR in the R environment to calculate z-score and p value by 579 permutation test [55]. We used the function "overlapPermTest" in regioneR with the following 580 arguments (A=peak1, B=peak2, ntimes=2000, alternative = "auto", genome=ce11, 581 non.overlapping=TRUE, mc.set.seed=FALSE).