Potential roles of DNA methylation in the initiation and establishment of replicative senescence revealed by array-based methylome and transcriptome analyses.

Cellular senescence is classified into two groups: replicative and premature senescence. Gene expression and epigenetic changes are reported to differ between these two groups and cell types. Normal human diploid fibroblast TIG-3 cells have often been used in cellular senescence research; however, their epigenetic profiles are still not fully understood. To elucidate how cellular senescence is epigenetically regulated in TIG-3 cells, we analyzed the gene expression and DNA methylation profiles of three types of senescent cells, namely, replicatively senescent, ras-induced senescent (RIS), and non-permissive temperature-induced senescent SVts8 cells, using gene expression and DNA methylation microarrays. The expression of genes involved in the cell cycle and immune response was commonly either down- or up-regulated in the three types of senescent cells, respectively. The altered DNA methylation patterns were observed in replicatively senescent cells, but not in prematurely senescent cells. Interestingly, hypomethylated CpG sites detected on non-CpG island regions ("open sea") were enriched in immune response-related genes that had non-CpG island promoters. The integrated analysis of gene expression and methylation in replicatively senescent cells demonstrated that differentially expressed 867 genes, including cell cycle- and immune response-related genes, were associated with DNA methylation changes in CpG sites close to the transcription start sites (TSSs). Furthermore, several miRNAs regulated in part through DNA methylation were found to affect the expression of their targeted genes. Taken together, these results indicate that the epigenetic changes of DNA methylation regulate the expression of a certain portion of genes and partly contribute to the introduction and establishment of replicative senescence.

Introduction

Cellular senescence is the irreversible cessation of cell proliferation [1] and is classified into two groups: replicative senescence and premature senescence [2]. Replicative senescence is caused by telomere shortening due to repeated DNA replication [3], while premature senescence is caused by stress, such as oncogene activation [4] and reactive oxygen species (ROS) [5], without apparent loss of telomere length and function. Although cellular senescence has been shown to be associated with tumor suppression in several cancers [2, 6, 7], it has been reportedly involved in cancer progression through the induction of epithelial-mesenchymal transitions and tumor invasion [8]. In addition, senescent cells secrete several factors associated with inflammation such as interleukin (IL)-6 and IL-8 [9], which are referred to as senescence-associated secretory phenotypes (SASP). Recently, SASP has been implicated in the pathogenesis of age-related diseases such as rheumatoid arthritis, periodontitis and Alzheimer’s disease [8, 10]. Therefore, elucidation of the mechanism contributing to induction and establishment of the senescent state will help overcome such age-related diseases.

Epigenetic regulation such as histone modification, DNA methylation and interference with micro RNA (miRNA) is one of the mechanisms that modulate gene expression. Alteration of the chromatin structure occurs in human senescent fibroblasts, where epigenetic regulation contributes to the establishment of the senescent state partly through the p16INK4A / retinoblastoma (RB) protein pathway [11]. Modification of histones, trimethylated histone H3 at lysine9 (H3K9me3) and trimethylated histone H3 at lysine 27 (H3K27me3) is enriched in the senescent-associated heterochromatic foci (SAHF), although the spreading of repressive histone marks is not necessary for SAHF formation [12]. So far as known, histone modification is reportedly involved in expression of key molecules of senescence, such as p16INK4A, p14ARF and p53. Loss of H3K27me3 is involved in expression of p16INK4A and p14ARF, whereas H4 acetylation and trimethylation of histone H3 at lysine 4 (H3K4me3) are involved in expression of p53 [13-16]. Furthermore, Takahashi, A. et al. [17] showed that the expression levels of SASP factors, IL-6 and IL-8, are regulated by demethylation of H3K9 through the APC/CCdh1-G9a/GLP pathway.

DNA methylation changes have also been observed during cellular senescence and individual aging. Many studies have elucidated the role of epigenetics in senescence using various approaches such as pyrosequencing, array-based methylome and combined bisulfite restriction analysis (COBRA), as well as using various types of cells and different tissues derived from genetically different individuals [18-21]. These studies have shown that DNA methylation profiles are tissue- and cell-type specific, and the epigenetic control of gene expression seems to promote in part tissue- and cell-type specific differentiation in addition to cellular senescence. For instance, the DNA methylation profile of fibroblast cells is different to that of mesenchymal stromal cells [22]. Moreover, principal component analysis showed that the methylation pattern of fibroblast cells from one dermal region is different to that from other dermal regions [23]. Christensen et al. also showed interindividual variation in methylation profiles among 11 tissues, including blood and brain [24]. In contrast to the studies using genetically different tissues and cells derived from individuals, age-related changes in DNA methylation between monozygotic twins were reported to arise with chronological time, indicating that the differences in genetically identical individuals are driven by different cellular responses to environmental changes [25]. Therefore, we hypothesized that a specific type of cultured cell leading to senescent states induced by different methods and demonstrating differences in methylation profiles will allow us to characterize in detail the role of epigenetic regulation in response to cellular senescence.

The effects of DNA methylation on gene transcription have been extensively studied in relation to the states of CpG islands (CGIs) near the transcription start sites (TSS). According to recent studies examining the relationship between gene body methylation and gene expression, hypermethylation correlated to not only high gene expression but also low gene expression [26], suggesting a more complex regulation. Varley et. al. also reported that the relationship between methylation and gene expression is context-dependent, although the current models reported by several groups indicate that methylation in the promoter regions is associated with gene silencing, and gene body methylation is associated with expression [20].

Furthermore, some senescence-associated genes are regulated in part by miRNA [27-31]. MiR-34 has been well known as a tumor suppressor and its targeted genes encode the cell cycle regulators including E2F, c-Myc, cyclin D1, cyclin E2, cdk4, and cdk6 [32-34]. Recent studies on cancers have revealed that the expression of some tumor suppressive miRNAs is regulated by DNA methylation in the miRNA promoter regions [35-37]. Despite these recent advances, much of the relationship between DNA methylation and miRNA expression and how this relates to senescence genes remains unclear.

In this study, we examined characteristic features of DNA methylation during cellular senescence using TIG-3 cells established from fetal lung fibroblasts. Although such TIG-3 cells have been a common focus in senescence research, much of the epigenetics remains unexplored. To test the hypothesis that the differences in genetically identical cells are driven by different cellular responses to senescence, we examined array-based gene expression and DNA methylation profiles using genetically identical TIG-3 cells which had been induced to a senescent state by three different methods, namely, replicatively senescent, ras-induced senescent (RIS), and senescent SVts8 cells. Originally derived from TIG-3 cells, SVts8 cells can induce senescence under non-permissive temperature by inactivation of the temperature-sensitive SV40 large T antigen. We then searched for the positional trend of methylation changes and the possible effects of DNA methylation changes on gene expression using integrated analysis.

Materials and methods

Cells and cell culture

Normal human diploid fibroblast TIG-3 cells (obtained from the Health Science Research Resources Bank, Japan) were cultured in DMEM+GlutaMAX-I (GIBCO) supplemented with 10% fetal bovine serum (FBS) (HyClone) and 1% penicillin/streptomycin (Nacalai Tesque) at 37°C under a 5% CO2 atmosphere. TIG-3 cells were cultured until they senesced at population doubling level (PDL) 85. The TIG-3 cells were harvested for DNA extraction at PDL 36, 49, 69 and 85, for RNA extraction to analyze gene expression at PDL 36 and 84, and for RNA extraction to analyze miRNA expression at PDL 44, 60, 78 and 80.

To prepare RIS cells, retroviral infection was performed as reported previously [38]. Briefly, Phoenix-Eco cells (obtained from Dr. G. P. Nolan, Stanford University, CA, USA) were transfected with the pBabe-puro-H-Ras-V12 or pBabe-puro plasmid by the Chen-Okayama method [39]. Viral supernatants were prepared from the cells after transfection, passed through a 0.45-μm-pore-size syringe filter, and pooled. The supernatant and 8 μg/mL hexadimethrine bromide (Sigma-Aldrich) were added to TIG-3 cells expressing an ecotropic receptor at a proliferating phase. After infection, the cells were selected with growth medium containing 300 μg/mL G418 and 2 μg/mL puromycin for 9 days before being harvested.

SVts8 cells (obtained from the Health Science Research Resources Bank, Japan) [40] continued to proliferate at a permissive temperature (33.5°C), because of suppression of RB and p53 through induction of a temperature-sensitive mutant of the simian virus (SV) 40 large T antigen to TIG-3 cells and the high ability of telomere maintenance. Senescent SVts8 cells were obtained by culturing SVts8 cells at a non-permissive temperature (38°C) for 6 days.

Microarray assays

Gene expression

Total RNA was isolated with a ReliaPrep RNA Cell MiniPrep System (Promega) according to the manufacturer’s instructions. Starting with 200 ng of the isolated RNA for each sample, double stranded cDNA and cyanine 3 labeled cRNA were synthesized using a low input quick amp labeling kit (one-color) and RNA spike-in kit (Agilent). The labeled cRNA was purified with an RNeasy mini kit (Qiagen), and hybridized to a SurePrint G3 Human GE microarray 8×60K Ver. 2.0 (Agilent). After washing the microarray to remove unhybridized cRNA, the microarray was scanned with an Agilent DNA microarray scanner G2505B, and then feature extraction was performed using the GE1_QCMT_Sep09 protocol.

DNA methylation

Genomic DNA was isolated using DNeasy Blood & Tissue (Qiagen) according to the manufacturer’s instructions. Genomic DNA (1.5 μg) was bisulfite-converted using the EpiTec Plus DNA Bisulfite Kit (Qiagen). From each sample, 300 ng of bisulfite-treated DNA was subjected to DNA methylation profiling using an Infinium HumanMethylation27 or HumanMethylation450 BeadChip array (Illumina) according to the manufacturer’s standard protocol. The array slides were scanned with an iScan system (Illumina).

MicroRNA (miRNA) expression

Total RNA including miRNA was isolated with a ReliaPrep miRNA Cell and Tissue MiniPrep System (Promega) according to the manufacturer’s instructions. Approximately 100 ng of the isolated RNA for each sample was labeled with cyanine 3 using miRNA Complete Labelling and Hyb Kit (Agilent) according to the manufacturer’s protocols. The labeled RNA was hybridized to a SurePrint G3 Human miRNA Microarray (Release 21.0) (Agilent). After washing, the microarray was immediately scanned using one color scan setting for 8×60k array slides with an Agilent DNA microarray scanner G2505C. The images were extracted with Feature Extraction Software 10.7.3.1 (Agilent) using default parameters.

Data analysis

The SurePrint G3 Human GE microarray data were analyzed using the Subio Platform (Ver. 1.18.4625). The data were normalized as to low signal cutoff (cutoff 1.0 and replace), log transformation (base 2), and global normalization (percentile 75), and then the ratios to those of the control sample (mean) were obtained. In this study, > = 2-fold and = < 0.5-fold changes were regarded as up- and down-regulated gene expression, respectively. In order to calculate fold-change differences and construct a heatmap, each gene expression level in the senescent cells was compared with that in the control cells, namely replicative senescent cells versus proliferating cells, RIS cells versus cells infected with the empty vector, and senescent SVts8 cells versus proliferating SVts8 cells. The heatmap with sample clustering was drawn using Subio Platform with Uncentered Correlation. The microarray data for gene expression have been deposited in the GEO (GSE81798).

The Infinium HumanMethylation27 BeadChip includes probes for 27,578 CpG sites [41]. The Infinium HumanMethylation450 BeadChip includes probes for 485,577 CpG sites covering 21,231 RefSeq genes (99%), and 26,658 CGIs (96%), and 3,091 probes for non-CpG loci [42]. The image data obtained by the iScan system were subjected to background subtraction and control normalization using GenomeStudio V2011.1 (Illumina). Methylation levels were calculated as β values (= intensity of the methylation allele / [intensity of the unmethylated allele + intensity of the methylated allele +100]), which ranged from 0 (0% methylation) to 1 (100% methylation). Probes for CpG sites with a detection p-value of > 0.05 or a blank β value were eliminated from further analysis. The β values of control cells were subtracted from those of senescent cells (Δβ). Δβ > = 0.2 and Δβ = < −0.2 were regarded as hyper- and hypomethylation in this study. TIG-3 cells at PDL 36 were regarded as control cells for those at higher PDLs (49, 69, and 85). An Infinium HumanMethylation27 BeadChip was used to obtain methylation profiles for all samples (TIG-3, RIS, and SVts8 cells). An Infinium HumanMethylation450 BeadChip was used to obtain more comprehensive methylation profiles of TIG-3 cells at PDL 36 and PDL 85, and RIS cells. The BeadChip data have been deposited in the GEO (GSE81788 and GSE81797). Scatter plots for β values were drawn using Genome Studio. Hyper- and hypo-methylated CpG sites were classified into six CpG subcategories depending on the location relative to the CpG island (CpG island, N_Shore, S_Shore, N_Shelf, S_Shelf, and the open sea) and into seven gene feature subcategories (− 1500 to − 200 bp upstream of the TSS (TSS1500), − 200 bp to 0 bp upstream of TSS (TSS200), 5’ untranslated region (5’UTR), first exon (1st exon), gene body, 3’ untranslated region (3’UTR), and intergenic region), according to the probe annotation (HumanMethylation450_15017482_v.1.1.csv) provided by Illumina [43].

Integrated analysis of DNA methylation and gene expression

Normalized gene expression and DNA methylation array data were integrated based on the genomic locations of the RefSeq genes’ TSS composed in the expression array and the genomic location of the CpG sites placed in the BeadChip array using custom perl scripts. RefSeq information was retrieved from http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz. We assessed Δβ values for the CpG sites located within 8 kb distance from the closest TSS of the RefSeq genes and the fold-change value of the corresponding gene expression. The gene expression data were obtained from probes annotated by the RefSeq ID for mRNA. When multiple expression probes existed in the same RefSeq ID, the mean of multiple intensities was calculated for the RefSeq ID. The distance from the methylation site to the TSS was calculated with a computer using the location information on the methylation probes and RefSeq (as described above). Promoters registered in RefSeq were classified into two classes, CGI and non-CGI promoters, using the following criteria: > = 50% GC content and observed-to-expected ratio of CpG > = 0.6.

Gene ontology (GO) analysis

GO analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 (http://david.abcc.ncifcrf.gov/) [15-17] with GOTERM_BP_ALL and functional annotation clustering. The genes for GO analysis were extracted according to each cut-off value given in the methylation and integrated analysis section. For GO analysis using only gene expression data, genes exhibiting a more than 3-fold change of expression were used due to the limited gene numbers loaded onto the annotation tool. The top 10 represented GO terms are shown in some tables. Enrichment scores of more than 1.3 gave p-values of less than 0.05.

To conduct GO analysis on gene feature- and CpG site-categories, the genes related to differentially methylated CpG sites on the gene-coding region based on the probe annotation were analyzed with GOTERM_BP_ALL and a functional annotation chart due to the small numbers of extracted genes. The ratios of the number of genes with similar functions were calculated for all gene feature- and CpG site-subcategories.

miRNA expression regulated by DNA methylation and its targeted genes

To measure the expression level of miRNAs during replicative senescence, the Human miRNA microarray data were analyzed using GeneSpring (Ver. 12.5). The data were normalized by a 90-percentile shift, and then the ratios to those of the control sample (PDL 44) were obtained. The cut-off for miRNA expression change used in this study was a ±1.5-fold change. The microarray data for miRNA expression have been deposited in the GEO (GSE90942). miRNAs exhibiting methylation changes were obtained from HumanMethylation450 BeadChip data using the “MIR” keyword in the “UCSC_REFGENE_NAME” column provided by Illumina. We searched for miRNAs with hypermethylated promoter regions and down-regulated expression, and vice versa, by comparing methylation changes (Δβ cut-off was ±0.2) with expression changes (fold-change cut-off, ±1.5). Genes targeted by those miRNAs for which expression seemed to be regulated by DNA methylation were picked up by TargetScanHuman Release 7.0. The resulting genes were then further filtered by more than two miRNAs and the gene expression levels based on our microarray data.

Results

Similar biological outcomes, but different gene expression patterns were detected in three types of senescent cells

We prepared three types of senescent cells, namely, replicatively senescent, RIS and senescent SVts8 cells. The senescence state in SVts8 cells was rapidly induced by the inactivation of the temperature-sensitive mutant of the SV40 large T antigen under non-permissive temperature. Senescence was confirmed by the growth arrest, appearance (a large flat morphology), senescence-associated beta-galactosidase (SA-β-Gal) activity, and protein and mRNA levels of p16 INK4A and p21Cip1/Waf1 (S1 Fig). Gene expression data for senescent cells obtained with SurePrint G3 Human GE microarrays (Agilent) showed that less than 10% of the probes tended to show similar changes in the three types of senescent cells (up-regulated: 358 of 4,355 probes in replicatively senescent cells, 3,927 probes in RIS cells and 5,124 probes in senescent SVts8 cells; down-regulated: 278 of 4,679 probes in replicatively senescent cells, 4,557 probes in RIS cells and 3,823 probes in senescent SVts8 cells), whereas more than 50% of the probes showed changes specific to each type of senescent cell (up-regulated: 2,396 of 4,355 probes in replicatively senescent cells, 2,384 of 3,927 probes in RIS cells, and 3,408 of 5,124 probes in senescent SVts8 cells; down-regulated: 2,561 of 4,679 probes in replicatively senescent cells, 2,750 of 4,557 probes in RIS cells, and 2,840 of 3,823 probes in senescent SVts8 cells) (Fig 1). A heatmap also showed different patterns of gene expression among the three types of senescent cells, even though sample clustering indicated that replicatively senescent cells and RIS cells were closer than senescent SVts8 cells (Fig 2).

Fig 1

The number of probes indicating up- or down-regulated genes in three types of senescent cells.

The number of up-regulated (A) or down-regulated probes (B) are shown in three types of senescent cells. The probes exhibiting a more than 2-fold change were counted. Purple circles, replicatively senescent cells; yellow circles, RIS cells; green circles, senescent SVts8 cells.

Fig 2

Heatmap of gene expression in three types of senescent cells.

Gene expression patterns in three types of senescent cells are shown as a heatmap with sample clustering using uncentered correlation.

According to the GO analysis of more than 3-fold changed genes, the down-regulated genes in the replicatively senescent, RIS and senescent SVts8 cells were mostly related to the “cell cycle” (S1 Table). On the other hand, the up-regulated genes were related to “immune response”, “locomotion”, and “cell migration” in all three types of senescent cells that we examined. Both up-regulated and down-regulated genes included “developmental process” genes. The results suggest that different genes with similar functions contributed to the senescent process, although the biological outcomes were similar among the three types of senescent cells.

Replicatively senescent cells showed DNA methylation changes

The DNA methylation profiles of senescent and proliferating cells were obtained using the Infinium HumanMethylation27 BeadChip. Among the three types of senescent cells examined, only the replicatively senescent cells (TIG-3 at PDL 85) were notably differentially methylated compared to the proliferating control cells (TIG-3 at PDL 36), whereas prematurely senescent cells (RIS and SVts8) were not (Fig 3, upper panels). Among the 629 and 366 CpG sites that were hyper- and hypo-methylated, respectively, in TIG-3 at PDL 85 (Fig 3, lower table), 565 (89.8%) and 310 (84.7%) sites showed a stepwise increase and decrease, respectively, of DNA methylation along with the progression of PDLs (36, 49, 69, and 85) (Table 1). Genes hosting the differentially methylated CpG sites detected in TIG-3 at PDL 85 were subjected to GO analysis using DAVID. Genes associated with hypomethylated CpG sites were found to be most enriched in the term “immune response” (enrichment score 7.86) and its related terms (S2 Table). On the other hand, genes associated with hypermethylated CpG sites were enriched in a wider variety of terms such as regulation of biological and developmental processes with lower enrichment scores (4.24 or lower).

Fig 3

Comparison of DNA methylation profiles in three types of senescent cells.

DNA methylation β values of senescent (x-axis) and proliferating control (y-axis) cells in the three types of senescent models are shown in scatter plots. The data were obtained by Infinium HumanMethylation27 (upper panels, A–C) and HumanMethylation450 BeadChip (middle panels, D–E). The correlation coefficient (r2) is shown in each plot. The lower table (F) shows the numbers of CpG sites that were hyper- (Δβ > = 0.2) and hypo-methylated (Δβ = < −0.2) upon senescence.

Table 1

Sequential changes of DNA methylation during replicative senescence.

		PDL 49 vs PDL 36	PDL 69 vs PDL 36	PDL 85 vs PDL36
Hypomethylated CpG sites	Total number	64	165	366
	Number of sites with sequential change	–	156 (94.5%)	310(84.7%)
Hypermethylated CpG sites	Total number	76	277	629
	Number of sites with sequential change	–	273 (98.6%)	565 (89.8%)

To more comprehensively examine the alteration of DNA methylation upon senescence, we obtained the DNA methylation profiles of TIG-3 cells at PDL 36 and PDL 85, and RIS cells using the Infinium HumanMethylation450 BeadChip. Again, replicatively senescent TIG-3 cells were differentially methylated compared to the proliferating control cells, whereas RIS cells were not (Fig 3D and 3E). Among 484,662 probes that passed quality control procedures, 14,214 and 12,727 probes were hyper- and hypo-methylated, respectively, in TIG-3 cells at PDL 85 as compared with TIG-3 cells at PDL 36 (Fig 3F). The ratios of differentially methylated probes were determined for each of the seven gene features and six CpG subcategories (Fig 4). The frequency of hypermethylated CpG sites was higher than that of hypomethylated CpG sites in TSS1500 (1.6 fold), TSS200 (2.1 fold), 5’UTR (1.5 fold), 1st exon (2.0 fold) and gene body (1.2 fold) subcategories (Fig 4A). As shown in Fig 4B, the frequency of hypermethylated CpG sites was higher in CGI (16.0 fold), N_Shore (2.2 fold), and S_shore (2.4 fold) subcategories, but was lower in N_Shelf (0.44 fold), S_Shelf (0.39 fold), and the open sea (0.51 fold) subcategories. We applied GO analyses to genes hosting the differentially methylated CpG sites to characterize the features of the genes that were supposed to be regulated in part by DNA methylation in each gene feature- and CpG site-subcategories. A total of 8,114 genes hosting the differentially methylated CpG sites were classified by GO terms into the subcategories. The GO terms obtained from functional annotation charts using DAVID were further categorized into seven groups: “immune response”, “metabolic process”, “transport”, “cell adhesion”, “development”, “signal transduction”, and “transcription”, plus an additional group, others (Figs 5A–5H and 6A–6H and S3 Table). The “immune response”-related genes hosting hypomethylated CpG sites were enriched in all of the gene feature subcategories, whereas a small portion of those hosting hypermethylated sites were located in the TSS1500, 5’UTR, 1st exon and gene body subcategories (Fig 5A). Interestingly, in the region from − 200 bp to 0 bp upstream of TSS (TSS200), all of the classified 63 genes related to “immune response” consisted of the hypomethylated CpG sites. Consistent with the results of GO analyses using HumanMethylation27 data (S2 Table), the genes hosting hypomethylated CpG sites were found to be markedly enriched with the genes involved in “immune response”. In sharp contrast, the “transcription”-related genes hosting hypermethylated CpG sites were enriched in all of the gene feature subcategories except for “1st exon”, whereas no “transcription”-related genes hosting the hypomethylated ones were located in the gene feature subcategories except “gene body” (Fig 5G). The results of the GO analyses conducted on the CpG site subcategories are shown in Fig 6. Genes related to “immune response” were enriched only in the genes hosting hypomethylated CpG sites in the “open sea”, but not those in other subcategories (islands, shores, and shelves) (Fig 6A). In addition, genes related to “transcript” were enriched in the genes hosting hypomethylated CpG sites in the “CpG islands” (Fig 6G). However, genes related to the other six groups of GO terms were enriched in the genes hosting hypomethylated CpGs in several CpG subcategories (Fig 6B–6F and 6H). When we focused on the genes hosting hypermethylated CpG sites, genes related to “immune response” were enriched in “N_shore” rather than in the “open sea”. In addition, genes related to other groups were enriched in the genes hosting hypermethlated CpG sites regardless of the CpG site locations (islands, shores, shelves, and the open sea) (Fig 6B–6G). Thus, a certain portion of “immune response”-related genes might be regulated in part by DNA methylation via CpG sites being outside of the CpG islands.

Fig 4

Rates of differentially methylated CpG sites in replicatively senescent TIG-3 cells in seven gene- feature (A) and six CpG- site (B) subcategories.

CpG sites showing Δβ > = 0.2 and Δβ = < −0.2 in senescent TIG-3 cells (at PDL 85) compared to those at PDL 36 were regarded as hyper- and hypo-methylated, respectively. Seven gene-feature subcategories: TSS1500, TSS200, 5’UTR, 1st exon, gene body, 3’UTR, and intergenic region. Six CpG subcategories: CpG island, N_Shore, S_Shore, N_Shelf, S_Shelf, and the open sea.

Fig 5

The characterization of genes hosting differentially methylated CpG sites in the gene feature subcategories.

GO analyses for gene hosting differentially methylated CpG sites were performed using DAVID (functional annotation chart and GOTERM_BP_ALL). The genes hosting hypermethylated (red) or hypomethylated (blue) CpG sites were classified into six gene feature categories (TSS1500, TSS200, 5’UTR, 1st exon, gene body and 3’UTR). Using the GO terms detected to be enriched (p-value = < 0.05) by DAVID, genes with similar functions were classified into seven groups: immune response (A), metabolic process (B), transport (C), cell adhesion (D), development (E), signal transduction (F), and transcription (G), plus an additional group, others (H). Histograms (A–H) show the ratio of the number of genes classified by GO terms in gene features subcategories (S3 Table shows the number of genes analyzed and classified).

Fig 6

The characterization of genes hosting differentially methylated CpG sites in the CpG site subcategories.

The same type of GO analyses shown in Fig 5 was applied to six CpG subcategories (N_Shelf, N_Shore, CpG island, S_Shore, S_Shelf, and the open sea). S3 Table shows the number of genes analyzed and classified.

Hypomethylation observed in the open sea was frequently associated with the up-regulation of genes related to immune response

Next, we performed an integrated analysis of HumanMethylation450 BeadChip and gene expression profiles to investigate potential functional effects of DNA methylation on gene expression in replicative senescence. Most methylation changes were unlikely to have a significant impact on gene expression (S2 Fig). Out of 212,885 probes located within 8 kb distance from a TSS of RefSeq genes, 1,596 CpG sites were differentially methylated along with the altered expression of the nearest genes to these CpG sites: 1,101 hyper- and 495 hypo-methylated CpG sites for which the nearest genes were down- and up-regulated, respectively, upon senescence (S4 Table). The genes nearest to these 1,596 differentially methylated CpG sites were subjected to GO analysis using DAVID. Those genes with up-regulated expression that were in a hypomethylated state in close proximity (within 8 kb) to their TSS were involved in “immune response” and “cell death” as ranked in the top 10 terms (Table 2). The immune response-related genes included MHC II, and the cell death-related genes included FAS and oxidized low density lipoprotein receptor 1 (OLR1). In contrast, those genes with down-regulated expression that were in a hypermethylated state in close proximity to their TSS were categorized into several categories such as “development” and “cell cycle”.

Table 2

GO analysis of genes showing the relationship between methylation and gene expression level in replicative senescence.

Rank	Hypomethylation & up-regulated gene expression(n = 351)		Hypermethylation & down-regulated gene expression(n = 516)
	GO term categories	Enrichment Score	GO term categories	Enrichment Score
1	Response to stimulus	3.20	Developmental process	7.08
2	Immune response	3.19	Organ development	5.14
3	Inflammatory response	2.73	Lung alveolus development	3.69
4	Response to hormone stimulus, Response to corticosteroid stimulus	2.50	Embryonic development	3.64
5	Response to nutrient levels, Response to extracellular stimulus	2.30	Cell cycle	3.16
6	Regulation of transport	1.69	Negative regulation of transcription	3.11
7	Developmental process	1.62	Epithelial tube morphogenesis	2.70
8	Cell death	1.37	Mesenchymal cell development	2.47
9	Actin cytoskeleton organization	1.26	Ear development	2.46
10	Negative/positive regulation of kinase activity	1.25	Positive regulation of transcription, DNA-dependent, positive regulation of metabolic process	2.34

* More than 1.3 of enrichment scores gave less than 0.05 p-value.

We also assessed the distribution of the distances between each of the 1,596 differentially methylated CpG sites and the nearest TSS. In this analysis, we classified genes into two subcategories based on promoter types, CGI and non-CGI promoters. The number of probes for all CpG sites within CGI and non-CGI promoters on the array were 204,829 and 94,432 respectively (Fig 7A). The number of probes located within 1 kb distance to the TSS was 186,998 (62%). The 1,101 hypermethylated CpG sites for which the nearest gene was down-regulated upon replicative senescence consisted of 744 and 357 CpG sites within CGI and non-CGI promoters, respectively. The hypermethylated CpG sites tended to be located more frequently in CGI promoters than in non-CGI promoters, and located close to the TSS: 630 out of 1,101 (57%) within 1 kb distance (Fig 7B). The distribution patterns were similar between all CpG sites (Fig 7A) and 1,101 hypermethylated CpG sites (Fig 7B). In sharp contrast, 495 hypomethylated CpG sites for which the nearest gene was up-regulated were mainly located in non-CGI promoters (Fig 7C). Within 1 kb distance of the TSS, 202 hypomethylated CpG sites, for which the nearest gene was up-regulated upon replicative senescence, consisted of 36 and 166 CpG sites within CGI and non-CGI promoters, respectively (Fig 7C). This distribution pattern was similar to that of a subset of hypomethylated CpG sites, for which the nearest gene was up-regulated and related to “immune response” (Fig 7D). This result suggests the possibility that DNA methylation on the promoter modulates the expression of a subset of “immune response” genes in replicative senescence.

Fig 7

Distribution patterns of the distances between CpG sites and the nearest TSS.

Distribution patterns for all 299,261 CpG probes on the HumanMethylation450 BeadChip (A), 1,101 hypermethylated CpG sites for which the nearest gene was down-regulated upon replicative senescence (B), 495 hypomethylated CpG sites for which the nearest gene was up-regulated upon replicative senescence (C), and 90 hypomethylated CpG sites for which the nearest gene was up-regulated upon replicative senescence and related to immune response (D) are shown.

DNA methylation affected the expression of miRNAs and the target genes during replicative senescence

A HumanMethylation450 BeadChip includes 3,436 probes related to miRNA. In replicatively senescent TIG-3 cells, 98 CpG sites corresponding to 66 miRNAs were hypomethylated, whereas 102 CpG sites corresponding to 62 miRNAs were hypermethylated (S5 Table). The expression levels of miRNAs were examined using Human miRNA microarray (Release 21.0), where 2,549 human miRNAs were represented. During replicative senescence, 178 miRNAs (fold-change > = ±1.5) showed up- or down-regulated (data not shown). To select miRNAs for which expression seemed to be regulated by DNA methylation, we compared the methylation changes of CpG sites related to miRNA with the miRNA expression changes. Among 18 miRNAs selected, seven miRNAs showed a decrease in expression accompanied with a hypermethylated CpG site in the promoters (S6 Table). This result was consistent with a previous report using IMR 90 cells, where the expression levels of six of the seven miRNAs were down-regulated during replicative senescence [27]. In contrast, no miRNAs showed an increase in expression along with a hypomethylated CpG site in the promoters. We next searched for candidate genes targeted by the seven miRNAs using TargetScan Human Release 7.0, and these genes were further selected by more than two miRNAs and the targeted gene expression levels based on our microarray data. As a result, we identified 27 genes for which expression seemed to be indirectly regulated by DNA methylation on the promoters of targeting miRNAs (Table 3). The genes encoding IL-6 signal transducer (IL6ST) and Zinc finger matrin-type 3 (ZMAT3) were included in the targeted genes.

Table 3

Targeted genes mediated in part by miRNAs, with miRNA expression regulated via methylation in the promoter regions.

Targeted gene symbol	Gene expression (fold-change)	miRNAs of regulating targeted gene
EIF4EBP2	3.0	hsa-miR-7-5p	hsa-miR-193a-5p	hsa-miR-335-5p
HECW2	4.2	hsa-miR-7-5p	hsa-miR-25-3p	hsa-miR-505-3p
CADM1	8.2	hsa-miR-7-5p	hsa-miR-505-3p
CNNM2	2.9	hsa-miR-335-5p	hsa-miR-505-3p
CPEB2	2.6	hsa-miR-7-5p	hsa-miR-505-3p
CRY2	2.2	hsa-miR-7-5p	hsa-miR-17-5p
DAAM1	2.5	hsa-miR-130b-3p	hsa-miR-335-5p
DGKH	2.5	hsa-miR-130b-3p	hsa-miR-505-3p
DOK6	2.2	hsa-miR-17-5p	hsa-miR-335-5p
EGR2	2.8	hsa-miR-17-5p	hsa-miR-25-3p
FAM134C	2.2	hsa-miR-17-5p	hsa-miR-335-5p
GRIN2A	2.0	hsa-miR-7-5p	hsa-miR-130b-3p
HPCAL4	7.6	hsa-miR-7-5p	hsa-miR-335-5p
IL6ST	2.2	hsa-miR-130b-3p	hsa-miR-505-3p
KLHL28	2.4	hsa-miR-7-5p	hsa-miR-335-5p
MEF2D	2.7	hsa-miR-335-5p	hsa-miR-505-3p
MYO1D	5.9	hsa-miR-193a-5p	hsa-miR-335-5p
NR4A3	2.1	hsa-miR-7-5p	hsa-miR-335-5p
OXR1	2.9	hsa-miR-7-5p	hsa-miR-17-5p
PGM2L1	4.8	hsa-miR-17-5p	hsa-miR-130b-3p
PIP4K2C	2.1	hsa-miR-25-3p	hsa-miR-505-3p
PPARGC1B	2.6	hsa-miR-7-5p	hsa-miR-505-3p
RAB11FIP5	2.8	hsa-miR-7-5p	hsa-miR-17-5p
RNF141	2.4	hsa-miR-7-5p	hsa-miR-335-5p
SEMA6D	5.9	hsa-miR-7-5p	hsa-miR-193a-5p
ZDHHC8	5.2	hsa-miR-17-5p	hsa-miR-335-5p
ZMAT3	3.4	hsa-miR-7-5p	hsa-miR-130b-3p

Discussion

In this study, we examined DNA methylation levels in three types of senescent cells and found that only replicatively senescent cells were differentially methylated, whereas prematurely senescent ones (RIS and SVts8) were not (Fig 3 and Table 1). These results were in good agreement with previous studies showing DNA methylation changes in aged tissues and cells [16, 21, 24, 44]. There are several reasons why the DNA methylation profile is strongly modified during replicative senescence and not during premature senescence. Firstly, errors may accumulate due to repeated cell division. Laird et al. evaluated the fidelity of transmission of the DNA methylation state in the CpG island of the FMR1 gene in normal human lymphocytes using hairpin-bisulfite PCR [45]. Although the high fidelity of inheritance of the methylated state of cytosine was estimated, the results clearly showed that errors in maintaining DNA methylation occurred to some extent in every DNA replication. When culturing TIG-3 cells from PDL 36 to PDL 85, the cells would be divided approximately 2^50 times. In contrast, as RIS and senescent SVts8 cells were rapidly induced to a senescent state, prematurely senescent cells would not have enough rounds of cell division to accumulate any errors. Secondly, reduced activity of DNA methyltransferase 1 (DNMT1), which primarily maintains DNA methylation patterns during replication, may increase the incidence of errors. The expression profiles in this study exhibited a marked decrease in DNMT1 in replicatively senescent cells (0.23-fold reduction, data not shown), whereas the expression levels of DNMT1 in RIS (0.70-fold reduction, data not shown) and senescent SVts8 (0.80-fold reduction, data not shown) were slightly reduced. This is consistent with the results reported by Kaneda et al., who reported that DNA methylation was not altered in RIS using methylated DNA immunoprecipitation (MeDIP) sequencing and bisulfite sequencing, and the Dnmt1 expression level was not altered in RIS, or during 3 passages (passage 2 to 5) using mouse embryonic fibroblasts [16]. Thirdly, a marked decrease in Ten-eleven translocation 1 (TET1) expression could alter DNA methylation patterns. Recent studies have suggested active DNA demethylation is mediated by TET, the enzyme that converts 5-methylcytosine (5mC) to 5-hydroxymethylcytosine (5hmC), 5-formylcytosine (5fc), and 5-carboxylcytosine (5caC) [46-49]. There are three TET proteins: TET1, TET2 and TET3. The expression levels of TET1 and TET3 genes, but not TET2, reportedly decrease along with aging, and TET3 expression is important for decreasing genomic 5hmC during aging in human T cells [50]. In replicatively senescent TIG-3 cells, the expression levels of the TET1 gene were drastically decreased, although we have no expression data for the TET3 gene due to a lack of TET3 probes in the expression array. A recent study showed that Tet1/Tet2 double-knockout mouse embryonic fibroblasts (MEFs) had defects in maintaining hypomethylation and resulted in hypermethylation of DNA methylation canyons where developmental genes are associated [51]. Taken togather, methylation changes during senescence would be required for many rounds of cell division and the decreased activities of DNA methyltransferase and methylcytsine deoxygenase, such as DNMT1 and TET1, might produce the altered methylation patterns within long-term culture.

We found that hypomethylation observed in the open sea was frequently associated with the up-regulation of genes related to immune response. When the genes hosting hypomethylated CpG sites were classified into six CpG subcategories using annotated GO terms, genes related to “immune response” were enriched only in the genes hosting hypomethylated CpG sites in the “open sea” (Fig 6A). Consistent with the results, the promoter types of genes categorized as “immune response” were mainly non-CGI promoters (Fig 7D). Nevertheless, the probes in the HumanMethylation450 BeadChip covered much more CpG island promoter genes than non-CpG island promoter genes (Fig 7A). In contrast, genes categorized into other groups, namely “metabolic process”, “transport”, “cell adhesion”, “development”, “signal transduction”, and “transcription”, exhibited hypermethylation rather than hypomethylation in the CpG sites close to the TSS (S3 Table). Furthermore, integrated analyses showed that “immune response” was ranked in GO terms of genes up-regulated and concomitantly hypomethylated in close proximity to their TSS (Table 2). These results suggest that hypomethylation in the open sea may increase the expression of the “immune response”-related genes during replicative senescence. For example, MHC II has non-CGI promoters exhibiting hypomethylated CpG sites and its expression was up-regulated during replicative senescence. Several studies also indicated the effects of DNA methylation in the promoter regions on the expression of inflammatory genes [52-54].

Currently we cannot explain the mechanism by which hypomethylation in the non-CpG promoter increases the expression of a certain portion of the “immune response”-related genes during replicative senescence. We speculate that TET2 may play a role in this regulation. As mentioned above, TET proteins contribute to active DNA demethylation. In our data, the expression levels of one of the TET2 probes showed a 2-fold increase in replicatively senescent cells and senescent SVts8 cells, whereas no increase in TET2 expression was detected in the RIS cells. Unlike TET1 and TET3, TET2 does not have a CXXC domain, which is required for binding to the CpG site [55, 56]. These results suggest that collaboration between TET2 and its associated factor(s) may be required for DNA binding and demethylation. In fact, TET2 was reported to need a cofactor for binding to DNA [49]. During differentiation of helper T (Th) cells, TET2 induced DNA demethylation at the loci of key cytokine genes in a lineage-specific transcription factor-dependent manner and promoted signature cytokine expression in Th1 and Th17 cells [57]. Depending on the cofactor, TET2 may bind to non-CGI promoters in senescence-associated genes, and increase gene expression. Further studies are needed to investigate the effects of TET2 on hypomethylation at the specific CpG sites and its interacting factor(s) during senescence.

Recently, senescence-associated miRNAs (SA-miRNAs) have been reported to regulate genes associated with senescence [27-31]. We examined the effects of DNA methylation in the promoter regions of miRNAs on miRNA expression using human miRNA expression microarray and further selected miRNAs that likely affected the predicted gene expression. As a result, we identified seven miRNAs and 27 targeted genes (Table 3). These genes included eukaryotic translation initiation factor 4E binding protein 2 (EIF4EBP2), which inhibits translation initiation [58], and cryptochrome circadian clock 2 (CRY2), which is a key component in regulating circadian rhythm [59]. Among the 27 genes, the IL-6 signal transducer (IL6ST) and zinc finger matrin-type 3 (ZMAT3) are involved in immune response. IL6ST, which is supposed to be up-regulated by decreased levels of has-miR-130b-3p and has-miR-505-3p, encodes glycoprotein 130 (gp130) [60]. The protein gp130 is a signal transducer shared by many cytokines including IL-6, one of the SASP factors [61], IL-11, IL-27, and oncostatin-M [62-64]. In addition, ZMAT3, a predicted gene regulated by has-miR-7-5p and has-miR-130b-3p, encodes a double-stranded-RNA-binding zinc finger protein Wig-1 (for wild-type p53-induced gene 1). Wig-1, a transcriptional target of p53, stabilizes p53 by binding to the 3’ UTR of p53 mRNA and protecting it from deadenylation [65]. A high level of p53 triggers cell cycle arrest, senescence and apoptosis, and efficiently inhibits tumor development [7, 66, 67]. In addition to the p53 transcription factor, miRNAs with expression regulated by DNA methylation via their promoter regions may also contribute to ZMAT3 expression. MiR-34 is one of the SA-miRNAs and is up-regulated by p53 [27, 28, 32, 33, 68, 69]. Although increased levels of miRNA-34 were detected upon replicative senescence, the methylation changes (Δβ cut-off was ±0.2) in the promoter regions of miRNA-34 were not included in our data (data not shown).

In this study, we investigated the possibility of regulation by DNA methylation during senescence. We found that hypomethylation in the open sea may contribute to the up-regulation of genes related to immune response. Several miRNAs targeting genes associated with a senescent state seem to regulate expression by DNA methylation in the promoter regions. However, we have to consider the possibility that DNA methylation results from alteration of gene expression. To investigate this possibility, we need to collect more data and explore the mechanism of methylation and demethylation change. Moreover, in order to reveal the whole mechanism of DNA methylation, we also need to focus on the methylation profile in hypomethylation and outside of the TSS and CpG islands.

Conclusion

Three types of senescent TIG-3 cells showed similar biological outcomes, but the regulatory mechanisms were different. Replicatively senescent cells showed sequential DNA methylation changes, but prematurely senescent RIS and SVts8 cells did not. In replicative senescence, hypomethylation with up-regulated gene expression often occurred in the open sea. Moreover, hypomethylation was observed in non-CGI promoters of genes related to the immune response. These results suggested that hypomethylation in the open sea regulates the expression of a certain portion of immune-related genes in replicative senescence. In addition, several miRNAs that seemed to have expression levels regulated in part by DNA methylation may also contribute to the expression of senescence-associated genes.

Confirmation of senescence.

Senescent cells were subjected to senescence-associated beta-galactosidase (SA-β-Gal) staining, qRT-PCR and immunoblotting. SA-β-Gal staining of the control A), C), E), replicatively senescent B), RIS D), and senescent SVts8 cells F). The percentages of SA-β-Gal-positive cells are shown at the bottom of each picture. Bar, 200 μm. Objective, ×10. G) The expression levels of p16INK4A and p21Cip1/Waf1 obtained with SurePrint G3 Human GE microarrays and qRT-PCR. H) Representative western blotting of p16INK4A and p21Cip1/Waf1, Ras, and loading control (actin). Images of p16INK4A are shown separately shown due to different exposure time.

(TIF)

Click here for additional data file.

Integrated analysis of methylation and gene expression in replicative senescence.

Each plot represents the values obtained from a single gene using integrated analyses. DNA methylation β values and gene expression levels are plotted along the abscissa and the ordinate, respectively. Red lines show the cut-off border. For methylation, Δβ for hypermethylation is > = 0.2, Δβ for hypomethylation is = < −0.2. For gene expression, the cut-off for increased/decreased expression was a ±2-fold change.

(TIF)

Click here for additional data file.

GO terms for up- or down-regulated genes in three types of senescent cells.

(XLSX)

Click here for additional data file.

GO terms for hypo- or hyper-methylated genes in replicatively senescent cells.

(XLSX)

Click here for additional data file.

The number of genes hosting differentially methylated CpG sites in the gene feature- and the CpG site- subcategories.

(XLSX)

Click here for additional data file.

Integrated analysis of gene expression and methylation changes in replicatively senescent cells.

(ZIP)

Click here for additional data file.

Hypo- or hyper-methylated miRNAs in the promoter regions of replicatively senescent cells.

(XLSX)

Click here for additional data file.

miRNA expression regulated by DNA methylation in the promoter region.

(XLSX)

Click here for additional data file.