Research ArticleBIOMOLECULES

Reorganization of chromosome architecture in replicative cellular senescence

See allHide authors and affiliations

Science Advances  05 Feb 2016:
Vol. 2, no. 2, e1500882
DOI: 10.1126/sciadv.1500882

Abstract

Replicative cellular senescence is a fundamental biological process characterized by an irreversible arrest of proliferation. Senescent cells accumulate a variety of epigenetic changes, but the three-dimensional (3D) organization of their chromatin is not known. We applied a combination of whole-genome chromosome conformation capture (Hi-C), fluorescence in situ hybridization, and in silico modeling methods to characterize the 3D architecture of interphase chromosomes in proliferating, quiescent, and senescent cells. Although the overall organization of the chromatin into active (A) and repressive (B) compartments and topologically associated domains (TADs) is conserved between the three conditions, a subset of TADs switches between compartments. On a global level, the Hi-C interaction matrices of senescent cells are characterized by a relative loss of long-range and gain of short-range interactions within chromosomes. Direct measurements of distances between genetic loci, chromosome volumes, and chromatin accessibility suggest that the Hi-C interaction changes are caused by a significant reduction of the volumes occupied by individual chromosome arms. In contrast, centromeres oppose this overall compaction trend and increase in volume. The structural model arising from our study provides a unique high-resolution view of the complex chromosomal architecture in senescent cells.

Keywords
  • Hi-C
  • genome organization
  • cellular senescence
  • chromatin
  • chromosome architecture
  • long-range interactions
  • chromosome conformation capture
  • DNA FISH
  • centromere
  • telomere

INTRODUCTION

Replicative cellular senescence is recognized as a fundamentally important biological process with roles in aging, tumor suppression, embryonic development, tissue repair, and wound healing (15). In many of these contexts, cellular senescence is believed to be triggered by genotoxic stresses, leading to irreparable DNA damage (for example, telomere dysfunction arising from replicative exhaustion, or replication fork collapse elicited by activation of oncogenes) (68). The accumulation of senescent cells in tissues has long been thought to contribute to age-related pathologies and functional decline, and recent evidence demonstrating beneficial effects of the removal of senescent cells supports this hypothesis (9).

The partitioning of chromosomes into accessible euchromatin and compact heterochromatin is central to the organization of the genome, and cellular senescence has been associated with numerous alterations to the chromatin landscape. Compacted chromatin domains, known as senescence-associated heterochromatin foci (SAHF), are formed in some but not all senescent cell nuclei (1014). DNA in heterochromatic regions becomes globally hypomethylated, whereas some CpG islands become hypermethylated (15). Expression of lamin B1 decreases in senescent cells and has been associated with a loss of peripheral heterochromatin, whereas centromeric regions become dramatically distended (1619). Using formaldehyde-assisted isolation of regulatory elements (FAIRE) (20), we found decreased chromatin accessibility in euchromatic regions, including gene promoters and enhancers, whereas heterochromatic regions became more accessible (16). It has been proposed that some of these chromatin changes contribute to the gene expression patterns characteristic of senescent cells and reinforce the irreversibility of the arrest (21).

These dramatic chromatin alterations are very likely to affect nuclear architecture, but the three-dimensional (3D) organization of chromosomes in senescent cells is not known. Recently developed high-throughput chromosome conformation capture methods, such as Hi-C and 5C, allow a high-resolution interrogation of 3D chromosome structure (22, 23). Hi-C experiments on interphase cells revealed a fundamental, large-scale organization of the genome into two types of compartments, designated as A and B (22, 24). These compartments are approximately 5 Mb in size and appear to be cell type–specific. Each chromosome is composed of multiple A and B compartments, which preferentially interact in cis (A with A, and B with B). The A compartments correlate with early replicating, euchromatic regions, whereas the B compartments correlate with heterochromatin (25). At higher resolution (<2 Mb), the genome is organized into topologically associated domains (TADs), which can be identified by the frequent interactions between loci located within these regions (2628). TADs can be located in either A or B compartments, and the activity of the genes in a TAD typically reflects its placement (that is, active in A and repressed in B). However, unlike compartments, TADs are highly conserved across cell types and even species, and are believed to constitute a fundamental organizational unit of the genome (26, 28). Supporting their organizational role, TAD boundaries are highly correlated with replication timing domains (29).

Several studies have explored changes in the 3D architecture of the genome during normal biological processes or under pathological conditions. For example, the fundamental organization of the interphase chromosomes breaks down in mitosis, where the genome adapts a homogeneous, locus-independent folding state (30). During differentiation, the basic interphase architecture is conserved, but frequent repositioning of TADs from A to B domains and vice versa is observed (31, 32). The repositioning between A and B compartments results in coordinated expression changes of the genes within the switching compartments and has been proposed to play a functional role in differentiation (31, 32). Here, we applied the genome-wide, high-resolution Hi-C approach, combined with direct measurements based on imaging methods, to investigate the 3D organization of chromosomes in senescent cell nuclei.

RESULTS

Reorganization of 3D chromosome interactions in senescence

Hi-C was performed on replicatively senescent, proliferating, and quiescent normal human diploid fibroblast (HDF) cells using three independent biological replicates. Quiescent cells were included as an additional control because they share features of cell cycle arrest with senescent cells (33). Sequencing reads were aligned to the genome by iterative mapping and filtered, binned, and corrected for biases using iterative correction (table S1) (25). Hi-C interaction matrices were normalized such that the sum of each row equals 1 and values in the Hi-C matrix represent the frequency of long-range interactions per bin. To allow comparisons across the experimental conditions, we resampled the Hi-C data sets at the fragment level such that an equal number of reads contributed to each Hi-C data matrix (table S2).

In proliferating cells, the Hi-C matrix displayed the checkerboard pattern characteristic of the previously described A and B compartments (Fig. 1A), and quiescent cells presented a very similar organization (Fig. 1B). In senescent cells, although the global A and B compartmentalization was qualitatively conserved, we observed a relative decrease in the long-range (>2 Mb) and an increase in the short-range (<2 Mb) interaction frequencies (Fig. 1C). These changes are particularly apparent in differential heatmaps (Fig. 1, D to F), where Hi-C frequency maps are subtracted between two conditions. The most pronounced differences were found between senescent and proliferating cells (Fig. 1E), followed by senescent versus quiescent cells (Fig. 1F), whereas quiescent and proliferating cells displayed the least differences (Fig. 1D). The same loss of off-diagonal (long-range) and gain of diagonal (short-range) signals in the Hi-C contact maps of senescent cells was observed in all chromosomes (figs. S1 to S4).

Fig. 1 Hi-C interactions matrices.

(A to C) Normalized heatmaps are shown for the q arm of chromosome 18 in proliferating, quiescent, and senescent cells. The color maps for relative interaction probability are displayed on the same scale for each heatmap. The PC1 signal was used to define the A and B compartments and is displayed below each heatmap (positive PC1 is shown in red and is designated as A compartments, and negative PC1 is shown in blue and is designated as B compartments). (D to F) Differential heatmaps are shown for the q arm of chromosome 18 for the indicated conditions. The color maps are displayed on the same scale for each comparison. Red is used to designate enrichment in the first condition (senescent or quiescent cells), and blue depletion.

To investigate this apparent redistribution of contacts in more detail, we calculated the mean intrachromosomal contact probability as a function of distance (24, 25). The mean contact probabilities of proliferating, quiescent, and senescent cells were all consistent with a power law (Fig. 2A), as previously observed (24), but notably, senescent cells displayed increased contact probabilities at ≤2 Mb and decreased contact probabilities at >2 Mb distances. This shift in contact probability can be quantified by computing the relative contact ratio of short-range versus long-range contacts, or SVL. We independently calculated the SVL ratio for each chromosome arm at a 40-kb binning resolution to better capture short-range contacts at <2 Mb, and found that the SVL was significantly higher in senescent cells when compared to either proliferating or quiescent cells (Fig. 2B).

Fig. 2 Quantitative comparison of short-range versus long-range contacts.

(A) Contact probability (Ps) was calculated as a function of genomic distance for interactions within individual chromosome arms across the whole genome. (B) The ratio of short-range (≤2 Mb) versus long-range (≥2 Mb) interactions (SVL) was calculated for each chromosome arm across the whole genome. The SVL is significantly higher in senescent cells compared to proliferating or quiescent cells (***P < 0.001).

We also examined in senescent cells the changes in mean contact probability as a function of distance at specific genomic features—gene promoters, lamin-associated domains (LADs), and regions with high GC content—using the approach of Zuin et al. (34). We observed an overall gain in short-range (<2 Mb) and loss of long-range (>2 Mb) interactions for all three genomic features (fig. S5, A to C). In addition, gene promoters and high GC content regions displayed a tendency to gain very short range (<400 kb) interactions, whereas LADs did not display this behavior (fig. S5, A to C). In aggregate, these analyses show that the most notable characteristic of senescent cells is a relative loss of long-range contacts and a gain in short-range contacts. A second independent Hi-C experiment performed with senescent cells showed that these changes are reproducible (as evidenced by hierarchical clustering of contact probabilities, A and B compartment profiles, and SVL ratios; fig. S6, A to C).

A and B compartment switching in senescence

We explored the definition of TADs using a high-resolution hierarchical method that captures sub-TAD interactions (35). We observed conservation of TAD boundaries across quiescent, senescent, and proliferating conditions (fig. S7A), which is consistent with the conservation of TADs across cell types (26). A and B compartments were defined using the first principal component (PC1) obtained by eigenvector decomposition of the Hi-C matrices, such that A is a positive PC1 and B is a negative PC1 (25). We then intersected ENCODE histone marks, replication timing domains, and LADs with the PC1 signal used to define A and B compartments. Using this approach, we found that constitutive heterochromatin marks globally overlap with B compartments in all conditions (fig. S7B). Euchromatin marks overlap with A compartments, and facultative heterochromatin marks overlap with a mixture of compartments (fig. S7C). Therefore, under all of our experimental conditions, we found a correlation between A and B compartments and active and repressed chromatin, respectively, in agreement with previous observations in other cell types and conditions (25).

We next examined whether the placement of TADs within the A and B compartments was changing. We defined TADs as being either in an A (active) or in a B (repressive) compartment by their mean PC1 signals, and found that overall, 636 of 5035 total TADs (12.6%) switched compartments in senescent or quiescent cells. A subset of TADs, 460 of 3196 (14.4%), defined as being A in proliferating cells, switched to B in senescent or quiescent cells (Fig. 3, A and C). Conversely, 176 of 1839 TADs (9.6%), defined as B in proliferating cells, displayed a switch in the opposite direction (Fig. 3, B and C). The number of TADs switching to an A compartment was more frequent in senescent cells than in quiescent cells.

Fig. 3 Switching of TADs between A and B compartments.

(A) Genome Browser view showing an example of a TAD that switches from an A compartment (red, positive PC1 signal) in proliferating and quiescent cells to a B compartment (blue, negative PC1 signal) in senescent cells. (B) Example of TADs switching from a B compartment in proliferating and quiescent cells to an A compartment in senescent cells. (C) Genome-wide summary of TAD switching using proliferating cells as a reference point. Venn diagrams show the number of TADs switching between proliferating, quiescent, and senescent conditions. Left: A-to-B compartment switches. Right: B-to-A compartment switches. (D) Same as (C) but displaying the number of genes that switch compartments. G1 to G6 designate gene sets containing the genes within each Venn diagram compartment. (E) GSEA analysis of a microarray expression data set from proliferating and senescent cells (33) using our gene set G2 (A to B switch in both quiescent and senescent cells). Significant overrepresentation of genes down-regulated in senescent/quiescent cells is evident [<0.001 false discovery rate (FDR)]. (F) The same analysis with gene set G3 (A to B switch only in senescent cells) showed overrepresentation of genes down-regulated in senescent cells (<0.001 FDR). (G) The same analysis with gene set G6 (B to A switch only in senescent cells) showed a significant overrepresentation of genes up-regulated in senescent cells (0.003 FDR).

The ENCODE project has used an integrative model of histone modifications and chromatin binding proteins to partition the genome into distinct regulatory states (36). To gain further insights into regulatory relationships, we intersected the TADs that switched between proliferating and senescent cells with the ENCODE-defined chromatin states. The switch from compartment A in proliferating cells to compartment B in senescent cells was significantly associated with several categories, including Polycomb repression, heterochromatin, insulator, weak transcription, weak enhancer, and weak promoter (table S3). A switch in the opposite direction (B in proliferating to A in senescent) was associated with Polycomb repression, heterochromatin, and insulator categories. Hence, switching TADs frequently overlap with regions of weak transcriptional and regulatory activity, and with facultative heterochromatin.

We next used the Hi-C PC1 signal to group genes into six gene sets according to their compartment switching behavior (Fig. 3D and table S4) and performed gene set enrichment analysis (GSEA) using a microarray data set from a published study that compared gene expression in proliferating, quiescent, and senescent HDFs (33). We found that genes switching from A to B compartments in both quiescent and senescent cells (gene set G2, Fig. 3D) or only in senescent cells (gene set G3, Fig. 3D) were significantly overrepresented for down-regulated genes in senescent cells (Fig. 3, E and F, and table S5). Conversely, for genes switching in the opposite direction, the corresponding gene set (gene set G6, Fig. 3D) was significantly overrepresented for genes up-regulated in senescent cells (Fig. 3G and table S5).

To investigate the contribution of compartment switching to the cellular phenotypes that characterize senescence, we examined the behavior of genes associated with the cell cycle, the senescence-associated secretory phenotype (SASP), the DNA damage response (DDR), p16 up-regulation, and the regulation of chromatin, using previously defined gene sets that encompass these pathways (3740). We found that most of the genes in these categories (63 to 80%) were in stable A compartments (table S6). However, senescence-specific A-to-B compartment switching (gene set G3) included 8 to 14% of the genes associated with cellular proliferation, p16 up-regulation, and chromatin (table S6). A similar proportion of genes in these categories were also present in the G2 gene set. In particular, the switching genes included key cell cycle regulators such as MCM6, RBL1, CENPF, CBX2, and CBX4 (fig. S8, A to D). We also observed overlap between B-to-A switching (gene set G6) and genes associated with senescence phenotypes (table S6), although to a lesser extent (1 to 4%). Two examples are the chromatin regulator TRIM24 and the SASP gene EGF (fig. S8, E and F).

Chromatin compaction in senescent cells

Hi-C does not provide measurements of physical distances between genomic regions nor can it address heterogeneity between cells. The preferential cis interactions between A and B domains (A with A, and B with B) should frequently position loci in different domains of the same type in closer physical proximity than indicated by the linear (genomic) distance between them, and fluorescence in situ hybridization (FISH) has been used to empirically verify the chromosome folding predictions of Hi-C (24, 41). To apply this validation method to our Hi-C data sets, we performed 3D DNA-FISH (42) to explore the actual separation between defined chromosome loci under our experimental conditions. We designed four FISH probes in one chromosome arm such that two probes (p1 and p3) were in two different A compartments, and the other two probes (p2 and p4) were in two different B compartments (Fig. 4A and table S7). We then verified that within each experimental condition, the distance between cis-compartment probes (p1-p3, p2-p4) was shorter in 3D space than that between intercompartment probes (p1-p2, p2-p3, p3-p4), despite intercompartment probes being linearly closer in the genome (Fig. 4, B to D, fig. S9, B to D, and table S8). These relationships are thus consistent with our assignments of genomic regions to the A or B compartment. The distances between intercompartment probes (A versus B) were, on average, significantly smaller in senescent cells than in quiescent cells (Fig. 4, C and D, fig. S9, C and D, and table S8). This was surprising because the distance between intercompartment probes was approximately 10 to 16 Mb, which is larger than the ≤2-Mb scale, where we observed increased contact probabilities in the Hi-C data. One possible explanation is that, globally, the chromosomes are becoming more compact in senescence.

Fig. 4 Physical distances between individual loci within a single chromosome arm.

(A) Four DNA-FISH probes (p1 to p4) were designed in the p arm of chromosome 4. Probes p1 and p3 are in nonadjacent A compartments, and p2 and p4 are in nonadjacent B compartments. Probes were chosen on the basis of strong A-A and B-B interactions in Hi-C data. (B) The schematic shows the 3D spatial relationship predicted by Hi-C for probes p2 and p4 (nonadjacent B compartments) and p3 (A compartment). (C) The distances between transcompartment (A-B) probes (p2-p3 and p3-p4) are significantly decreased in senescent cells (***P < 0.001). (D) Representative 3D DNA-FISH images of quiescent (upper panel) and senescent (lower panel) cells.

To test this hypothesis, we first investigated global chromatin accessibility in senescent cells using several complementary methods. The FAIRE method is based on the differential extraction of relatively open chromatin (20). We performed quantitative experiments to determine absolute yields of FAIRE-enriched DNA and obtained values that were significantly lower for senescent cells relative to proliferating or quiescent cells (Fig. 5A). We next performed deoxyribonuclease (DNase) I nuclease sensitivity assays on intact nuclei, a well-established method to assess chromatin compaction. We observed that compared to quiescent or proliferating cells, senescent cells were significantly more resistant to DNase I treatment (Fig. 5B and fig. S10). Hence, both the FAIRE and DNase I sensitivity results are indicative of globally decreased chromatin accessibility in senescent cells, which is consistent with chromosome compaction.

Fig. 5 Assessment of relative chromatin accessibility.

(A) FAIRE experiments were performed on three separate occasions using cells in the indicated conditions as starting material. Yields of FAIRE extracted (soluble) DNA are shown as picogram per cell. Relative to input DNA, these values represent 12, 10, and 4%, respectively. Error bars show SDs (**P < 0.01; *P < 0.05). (B) DNase I sensitivity of intact nuclei was visualized using the comet assay.

To directly assess the volume occupied by individual chromosomes, we performed chromosome painting. Using 3D-preserving fixation conditions (42), we found significant decreases of chromosome 4 and 18 volumes in senescent cells compared to quiescent cells (Fig. 6, A and B, and fig. S11). These data provide further direct support for the compaction of senescent chromosomes indicated by other methods. From these experiments, we also calculated the mean chromosome radii and combined these parameters with the Hi-C contact probabilities to build scaled 3D models of chromosomes. The resultant integrative 3D chromosome models (Fig. 6C and movies S1 and S2) provide a novel visualization of the structures of human chromosomes and their compaction in senescent cells. Consistent with what has been previously reported, the A and B compartments segregate on opposite hemispheres of the chromosome (43). In senescent cells, this global organization does not change, but the chromosome arms are significantly reduced in volume.

Fig. 6 Compaction of chromosomes in senescent cells.

(A) 3D chromosome painting of chromosomes 4 and 18 in quiescent and senescent cells. 3D renderings of Z stacks of images were generated with Imaris software. Representative nuclei are shown for each condition (quiescent, senescent) for chromosome 4 (top), chromosome 18 (middle), and both visualized in the same cells (bottom). (B) Chromosome volumes were calculated from the 3D renderings, and the resulting distributions of the volumes are shown as box plots. Chromosomes in senescent cells had a significantly smaller volume (***P < 0.001; **P < 0.01). (C) 3D modeling of chromosome 18 based on Hi-C contact probabilities and mean chromosome radii from chromosome painting as scaling factors. The colors designate A (red) and B (blue) compartment signals. (D) In the collapsing spring model, chromosome arms shrink in size as a consequence of an increased local compaction of the chromatin. Increased contact probability in TADs observed in senescent cells is consistent with a model in which the intra-TAD distances decrease more than the inter-TADs distances (d, intra-TAD distance; D, inter-TAD distance; q, quiescent; s, senescent). TADs are depicted here as spheres, shaded either red for A compartments or blue for B compartments.

Decompaction of satellite DNA in senescent cells

An increase in the volume of centromeric and pericentromeric regions in senescent cells, observed by FISH using satellite DNA probes, has been noted by several groups (15, 16, 19). However, the probes used for chromosome painting experiments are designed to target unique DNA sequences and exclude centromeric regions (as well as other repetitive sequences). Likewise, in Hi-C, repetitive, multimapping reads are excluded from the analysis. To obtain a direct assessment of centromeric volume changes during senescence, we performed 3D DNA-FISH using α-satellite DNA probes. Similar to previous observations, and unlike the chromosome arms, α-satellite DNA displayed an increase in volume in senescent cells, which we quantified as 1.72-fold (fig. S12). In contrast, we observed that telomeres decreased in volume by 1.73-fold. In aggregate, these data illustrate the opposing changes that affect chromosome architecture during cellular senescence: a decrease in the volume of chromosome arms and an increase in the volume of centromeres.

DISCUSSION

Here, we applied a combination of whole-genome Hi-C, biochemical, cytological, and computational methods to characterize the changes in 3D chromosome architecture during replicative senescence of human cells. This strategy is important because conformation capture methods, such as Hi-C, can only provide relative comparisons, whereas physical measurements made by the other methods are of lower resolution when used in isolation (41). The data provided by all the methods were in general agreement and suggest a model whose salient feature is a decrease in the volume and increase in the compaction of chromosomes during cellular senescence. The only regions that escape this trend are the centromeres, which, in contrast, increase significantly in volume.

The decrease in volume was indicated by chromosome painting experiments (Fig. 6, A and B, and fig. S11) and further substantiated by measurements of distances between specific loci made by 3D DNA-FISH (Fig. 4, C and D, and fig. S9, C and D). An apparent decrease in chromosome size evidenced by painting has been previously noted (13). The increase in compaction, which is inversely associated with accessibility, was further suggested by quantitative FAIRE yield and DNase I sensitivity experiments (Fig. 5, A and B). The opposite behavior of centromeric and pericentromeric regions was visualized by 3D DNA-FISH and has been previously reported by several groups (15, 16, 19).

This decrease in the volume occupied by chromosome arms stands in contrast to the known increase of nuclear volume during senescence (44, 45). Indeed, in the experiments reported here, the mean volume of senescent cell nuclei (663.2 μm3) was 2.23-fold larger relative to quiescent nuclei (297.2 μm3). This apparent discrepancy can be reconciled by noting that despite the extraordinary length of the DNA in the nucleus, chromatin is estimated to account for only some 15% of the nuclear volume (46). We have previously reported that total nuclear protein content increased 1.93-fold during senescence, whereas the DNA content remained relatively constant (45). Hence, although the volume of chromosomes appears to decrease during senescence, the fraction of total nuclear volume occupied by DNA is relatively minor and is overshadowed by a large increase in the amount of total nuclear protein.

The Hi-C data sets allowed the partitioning of the genome into A and B compartments at 200-kb resolution using conventional computational methods (25). Direct measurements by 3D DNA-FISH of distances between specific loci in A or B compartments were used to verify the validity of the assignments generated by our computational pipeline (Fig. 4, A to C). In a general sense, the 3D organization of chromosomes revealed by Hi-C was well conserved between proliferating, quiescent, and senescent cells (Fig. 1, A to C). Senescent cells were, however, uniquely characterized by a relative decrease in long-range interactions (>2 Mb) and an increase in short-range interactions (<2 Mb; Fig. 1, D to F).

At a higher resolution (40 kb), TADs were clearly visible and were well conserved between proliferating, quiescent, and senescent conditions. We found that 13% of all TADs switched compartments between senescent and quiescent cells with respect to proliferating cells. A similar TAD switching behavior has been observed during cellular differentiation of embryonic stem cells to a variety of lineages, encompassing 3.8 to 25% of the genome (31), or from pre–pro-B cells to pro-B cells, encompassing 12% of the genome (32). Because the transitions between the cellular states analyzed here are also characterized by an overall high conservation of TADs accompanied by frequent switching between compartments, they resemble cellular differentiation events.

The overlap between quiescence and senescence was 181 TADs (Fig. 3C). Senescence was associated with a much higher frequency (3.7-fold) of switching from B (inactive) to A (active) compartments. Hence, quiescence and senescence display unique and overlapping components, with senescence representing a significantly more altered state. Gene expression within the switching TADs was strongly correlated with the direction of the switch: genes that switched from A to B compartments were mostly repressed and vice versa (Fig. 3, E to G). Further examination revealed that many of the genes that are characteristic of the senescence phenotype were located in stable A compartments under all conditions examined (table S6). These are thus likely regulated independently of TAD localization, such as by transcription factors binding to their promoters or enhancers within permissive A compartments. For example, the transcription factor nuclear factor κB (NF-κB) has been identified as the master regulator of SASP genes (47), and these genes are mostly located in stable A compartments.

Genes switching from A to B compartments in senescent cells (gene set G3, Fig. 3D) comprised a significant proportion of genes in the cell proliferation (13.6%), chromatin (8.2%), and p16 up-regulation associated (13.3%) gene sets (table S6). A similar proportion of genes (13.6%) associated with cell proliferation switched in common in both quiescent and senescent cells (gene set G2), which is not unexpected because both states share the feature of cell cycle arrest. For example, MCM6, which is involved in DNA replication, switches to a repressed state in both quiescent and senescent cells (fig. S8A). Hence, the A-to-B compartment switching of cell cycle genes in senescent cells could play a role in regulating and maintaining their permanent cell cycle arrest.

Recently, Chandra et al. (48) also used Hi-C to analyze chromosome conformation changes during oncogene-induced senescence (OIS) to provide insights on the formation of SAHF. Their model, involving acute expression of a RAF1-ER oncogene in WI38 cells, was optimized for high (86%) efficiency of SAHF formation (48). It is well known that the frequency of SAHF formation is quite variable between different strains of HDFs, and SAHF formation is relatively infrequent during replicative senescence (49, 50). In the experiments reported here, we observed SAHFs in approximately 22% of the cells. In addition, entry into senescence in response to telomere shortening is asynchronous and occurs over a period of several weeks. Hence, Chandra et al. harvested their cells shortly after the induction of senescence (2 days), whereas we kept our cells in a senescent state for several weeks to allow a full development of the phenotype.

We processed the raw Hi-C data of Chandra et al. through our computational pipeline and compared it to our replicative senescence data set (fig. S13). The one aspect where the results differ is the pattern of long-range interactions: whereas we see a uniform loss of long-range and gain of short-range interactions along an entire chromosome, the contact matrices of Chandra et al. show interspersed regions at the TAD level, where short-range interactions are lost and long-range interactions are gained (fig. S14). They proposed that such focal increases in long-range interactions represent SAHF formation in these regions.

Several interpretations are possible when comparing these two studies. If the focal long-range interactions represent SAHFs, they could be missing from our data simply because our cells form SAHFs infrequently, and features associated exclusively with SAHF would generate a weak signal. Hence, different senescence states could trigger somewhat different genome structural endpoints. Alternatively, SAHF could represent an intermediate step that resolves at later times into the global shrinkage that we see. In agreement, the global increase in short-range interactions we detected was much more prominent in regions of high GC content (fig. S5), which Chandra et al. also linked to SAHF formation. If this were the case, then the genomic features we observe would be representative of mature, late-stage senescence.

An interesting additional comparison is a Hi-C study of Hutchinson-Gilford progeria syndrome (HGPS) fibroblasts (51). HGPS is characterized by accelerated aging phenotypes, which are caused by mutations in lamin A/C that disrupt the nuclear lamina (52). HGPS fibroblasts in culture display misshapen nuclei, detachment of LADs from the nuclear lamina, and accelerated senescence. Hi-C analysis of late-passage HGPS fibroblasts showed an almost complete loss of long-range interactions between A and B compartments (51). In contrast, although we observed that in senescent cells, chromosomes became more compact and there was a relative decrease in long-range contacts, for the most part compartment-level organization was preserved. The features observed in HGPS cells thus appear to be a more drastic version of the changes we found in replicatively senescent cells, and point to the importance of the interactions between chromatin and the lamina for maintaining chromosome structure.

In our data, the increase in short-range interactions in senescent cells was particularly evident in the differential matrix plots as an increase of signals close to the diagonal that were opposed by a relative loss of long-range interactions between TADs (Fig. 1, E and F). These contrasting behaviors, combined with the global reduction of chromosome volume, are consistent with a collapsing coil model (Fig. 6D). The model envisions that DNA regions harboring more frequent shorter-range interactions would experience a higher level of condensation than regions where interactions occur at much larger distances and are less frequent. This would cause chromosome arms to shrink nonhomogeneously at different scales, with sub-TAD regions shrinking more than inter-TAD distances. At this time, we can only speculate what leads to this new configuration. A global loss of lamin B interactions could promote compensatory short-range contacts. Other, hitherto poorly understood processes could also contribute: for example, histone H1 is replaced with the HMGA2 protein, which could have an important effect on chromatin architecture (53, 54).

In summary, although the underlying mechanisms remain to be elaborated in detail, our studies presented here show that senescent cells acquire a unique chromosome architecture characterized by an increase in short-range chromatin contacts accompanied by a genome-wide shrinkage of chromosome arms. More localized compartment switching is correlated with changes in gene expression, supporting the hypothesis that the 3D organization of chromosomes contributes an additional level of regulation to the cell cycle arrest (and other phenotypes) of senescent cells.

MATERIALS AND METHODS

Cell culture

The HDF cell strain LF1 (55) was grown at 37°C in an atmosphere of 5% CO2 and 2.5% O2 in F-10 nutrient medium (Thermo Scientific) supplemented with 15% fetal bovine serum (FBS), penicillin, streptomycin, and l-glutamine (7). Cells were serially passaged by trypsinization at 1:4 dilution after reaching 80 to 90% confluence. Proliferating cells were propagated under exponential growth conditions for two population doublings before harvesting. To obtain quiescent cells, proliferating cells were grown to 80% confluence, serum supplementation of the medium was changed to 0.1% FBS, and incubation was continued for 48 hours before harvest. To obtain senescent cells, cultures were propagated to replicative exhaustion and then further incubated for 4 months before harvesting. During the entire incubation under senescent conditions, complete medium was used (15% FBS supplementation) and was replaced every 4 days. Senescence was verified by staining for senescence-associated β-galactosidase (SA-β-Gal) (56). Senescent cultures contained a minimum of 93% of cells that stained positive for SA-β-Gal activity.

Construction of Hi-C libraries

Starting material was cells harvested under three different biological conditions: proliferating, quiescent, and senescent, as noted above. For each condition, three independent experiments were performed starting with a new inoculum of cells, for a total of nine independent harvests, which were then processed into nine libraries. Hi-C libraries were prepared using Hind III as described (57) with small changes to improve the final yield. Briefly, 1 million cells were fixed, lysed, and digested overnight with Hind III. After fill-in with biotin-dCTP (deoxycytidine 5′-triphosphate), the reactions were ligated, de–cross-linked, and purified using phenol chloroform extraction. Biotin was removed from unligated ends, and the DNA was sheared into 100– to 300–base pair (bp) fragments with a Covaris S2 instrument. The sonicated DNA was ethanol-precipitated in the presence of glycogen, resuspended, and end-repaired using the End-It DNA End-Repair Kit (Epicentre). The DNA was purified using Agencourt AMPure XP paramagnetic beads (Beckman Coulter) and eluted in molecular grade water. Deoxyadenosine triphosphate (dATP) was added to the DNA ends (58), and after another DNA purification, biotin-labeled DNA was pulled down with streptavidin beads (Dynabeads MyOne Streptavidin C1). Preannealed sequencing adapters were ligated to each sample (59). Thirteen cycles of PCR amplification were performed using Phusion High-Fidelity DNA Polymerase (New England Biolabs). A final purification step with AMPure XP beads yielded libraries in the 200- to 500-bp range that were sequenced on an Illumina HiSeq 2500 instrument in the Brown University Genomics Core. Digestions with Nhe I and Hind III to quantify the frequency of introducing new chimeric restriction sites in the fill-in reactions were performed as quality controls. Furthermore, aliquots of the final libraries were cloned using the Zero Blunt TOPO PCR Cloning Kit for Sequencing and Sanger-sequenced (50 clones per library) to evaluate the library content. The results of this validation are provided in table S9. The second independent set of Hi-C experiments for senescent cells used the protocols described above and two biological replicates.

Other molecular methods

For the preparation of FAIRE DNA, cells were grown in 15-cm culture dishes and collected at 80 to 90% confluence, as described in detail previously (16). Cells in culture dishes were washed twice with 5 ml of prechilled phosphate-buffered saline (PBS) and harvested by scraping into 1.5 ml of PBS containing 1 mM phenylmethylsulfonyl fluoride. While keeping the tubes on ice, harvested cells were counted with a Countess Automated Cell Counter (Invitrogen), 1% of the cells was set aside for preparing input DNA, and the remainder was immediately cross-linked by adding 37% formaldehyde to a final concentration of 1% for 4 min at room temperature with constant tube rotation. FAIRE DNA extraction from this point was performed as previously described (16). The final yield of FAIRE DNA, as well as of total (input) DNA, was quantified with a Qubit 2.0 instrument and the dsDNA High Sensitivity Assay Kit. FAIRE DNA yield was also normalized to the starting cell number. Cells were grown, passaged, harvested, and processed on three separate occasions.

For 3D DNA-FISH, four fluorescently labeled bacterial artificial chromosome (BAC) clone DNAs (RPCI-11-792J24, RPCI-11-49D20, RPCI-11-427M10, and RPCI-11-812M17) were purchased from Empire Genomics and designated as p1 to p4, respectively. Multicolor FISH on 3D-preserved nuclei was performed as described (42), with the exception that the hybridization buffer supplied by the probe manufacturer was used. To ensure full probe penetration into senescent nuclei, the optional pepsin treatment was included. Nuclei were counterstained with 4′,6-diamidino-2-phenylindole, and coverslips were mounted using ProLong Gold antifade medium (Life Technologies). Images were acquired using an LSM710 confocal microscope (Zeiss), visualizing three probes at one time (p1, p2, and p3, or p2, p3, and p4). Images were analyzed and distances were measured using the LSM Image Browser (Carl Zeiss) and Imaris software (Bitplane). Imaris software was also used to produce 3D renderings of the Z stacks of images in the 3D-DNA FISH experiments.

For whole chromosome painting, the same 3D-preserving protocol was used, substituting Kreatech FISH Probes (Leica Biosystems) for the BAC DNAs. Whole chromosome probes for chromosomes 4 and 18 were direct-labeled by the manufacturer with the PlatinumBright495 and PlatinumBright550 fluorophores, respectively. The hybridization buffer was supplied by the manufacturer. To visualize centromeres, Cy3-labeled α-satellite (CTTTTGATAGAGCAGTTTTGAAACACTCTTTTTGTAGAATCTGCAAGTGGATATTTGG) or Cy5.5-labeled satellite II (ATTCCATTCAGATTCCATTCGATC) probes were used (19). Telomeres were identified using a telomere-specific peptide nucleic acid (PNA) probe (7). Quantification of chromosome painting experiments was performed with Imaris software to measure chromosome volumes. Telomere and centromere volumes were measured in a similar manner. Imaris software was used to produce 3D renderings of the Z stacks of images in the 3D chromosome painting experiments.

For nuclease sensitivity experiments, cells were grown in 10-cm dishes, harvested by trypsinization, and washed twice with 10 ml of buffer A1 (60). Cells were collected by centrifugation at 2500 rpm for 5 min. The cell pellet was resuspended in 2.5 ml of buffer A2 (A1 + 1.0% NP-40), and cells were lysed by repeated pipetting with a 200-μl pipette tip for 10 min on ice. Nuclei were collected by centrifugation at 1000 rpm for 10 min. Nuclei were washed once with DNase I buffer (60) and resuspended to a final concentration of 5 × 106 nuclei/ml. Nuclei (1 × 105) were placed in an ice water bath, CaCl2 was added to a final concentration of 1 mM, and DNase I was added to a final concentration of 50 U/ml. The nuclei were incubated in the presence of DNase I on ice for 1 hour. The tubes were placed in a room temperature water bath for exactly 10 min and then placed back in the ice water bath for 30 s. The reactions were stopped by adding EDTA to 10 mM. Single-cell electrophoresis (comet assay) was used to analyze the extent of DNase I digestion as previously described (60). Images were analyzed with ImageJ open source software (http://rsbweb.nih.gov/ij/). The extent tail moment and the olive tail moment were calculated for each nucleus.

Computational analysis of Hi-C data

Sequencing reads were aligned, processed, and iteratively corrected as described (25), using the hiclib library for Python (commit b9a7db9, https://bitbucket.org/mirnylab/hiclib). Briefly, paired-end reads were mapped to the human reference genome (hg19) with bowtie2 using the “very-sensitive” option. For iterative mapping, reads were first truncated to 25 bp and mapped. Nonuniquely mapping reads were then extended by 5 bp for each iteration until the maximum length of 100 bp. Pairs of mapped reads were assigned to Hind III restriction fragments based on the 5′ alignment position. Restriction fragment–level data were filtered to remove artifacts from the Hi-C and sequencing protocols; specifically, single-sided reads, reads mapping within 5 bp of a restriction site, duplicate double-sided reads, reads from fragments smaller than 100 bp, reads from fragments larger than 100 kb, and reads from the top 0.5% of detected fragments were removed. After all alignment and filtering steps, senescent and quiescent samples had more sequencing reads than the proliferating samples. Therefore, each condition was resampled to the level of the proliferating sample. This process was optimized such that after fragment-level filtering, equal numbers of reads contributed to the final genome-wide interaction matrices for each condition and resampling was done independently three times (table S2). Fragment-level data were binned at a specific resolution by dividing the genome into fixed-length intervals. A 200-kb bin size was used for global examination of long-range contacts, and a 40-kb bin size was used for examination of local domain-level interactions. Iterative correction of the binned interaction matrices was done as described (25). Uninformative contacts on the diagonal and in neighboring bins, bins with less than 0.5*(binning resolution) sequenced counts, bins with the bottom 1% of counts, and bins with the top 0.05% of interchromosomal counts were removed. After filtering, iterative corrections were performed such that the rows of the interaction matrix approached a constant. Finally, to produce the normalized Hi-C interaction matrices, each matrix was divided by the mean row sum constant such that the sum of each row was equal to one.

Hi-C interaction heatmaps were produced using the normalized contact probabilities for each condition. Differential Hi-C interaction maps were calculated in a manner similar to Zuin et al. (34), subtracting the second condition from the first (using the 200-kb normalized interaction matrices). Hence, in the resultant differential matrix, positive values indicate higher and negative values indicate lower contact probability in the first condition. Contact probability as a function of distance (Ps) was calculated using the function included in the hiclib package (25). Only intrachromosome arm interactions were considered in the calculation of Ps. The short-range versus long-range, or SVL ratio, was calculated asEmbedded ImageOnly cis interactions within a single chromosome arm were considered for the calculation. To better capture the interactions in the <2-Mb range, we used the normalized interaction matrices at 40-kb resolution. Genomic feature analysis of contact probabilities was performed as described by Zuin et al. (34), using our high-resolution Hi-C interaction frequency data binned at 40 kb. The Hi-C interaction maps were used to calculate the mean change in contact probability as a function of genomic distance. We used only intrachromosomal arm contacts and separated contact pairs into three categories: 2×, for interactions where both bins contain the feature; 1×, for interactions where only one bin contains the feature; and None, for interactions where neither bin contains the feature. The features used in this analysis were RefSeq gene promoters (UCSC hg19, version update 6 October 2015), LADs from TIG3 fibroblasts (61), and high GC content (defined as >46% GC content using bedtools) (48, 62).

A/B compartment calling was performed as described (25) on the normalized interaction matrices. We used the hiclib package to calculate the eigenvector decomposition of the 40-kb interaction matrix, considering only intra-arm interactions. The first eigenvector (PC1) was used to compute the A and B compartments. Bins with PC1 values >0 were designated as A, and bins with values <0 were designated as B. The multiscale TAD calling approach was used to define TADs (35). TAD identification is dependent on the γ parameter for this algorithm. TADs were defined as the conserved set of nonoverlapping domains identified using a range of γ values. For additional analysis of the TADs, we used the TAD definition from the highest-resolution data set (quiescent cells). To identify TADs that switched between compartments, the mean first eigenvector of the PC1 signal for each defined TAD was calculated for each biological condition. TADs that converted from a positive mean PC1 (A compartment) to a negative mean PC1 (B compartment) in quiescent or senescent cells with respect to proliferating cells or the reverse (negative PC1 to positive PC1) were designated as TADs that switched compartments. To test for enrichment of switching TADs with ENCODE-defined regulatory states, we selected a comparable cell type or normal human lung fibroblasts and then used the Bits algorithm to test for significant overlap using randomization (36, 61, 63, 64).

To perform 3D modeling, we used a variant of the multidimensional scaling (MDS) optimization algorithm defined by Varoquaux et al. (65) as MDS2 with two notable changes. First, we added constraints to the MDS2 model that restrict the ratio of the inferred distances between FISH probe containing bins to be approximately equal to the ratios of the distances measured in the 3D DNA-FISH experiments. Second, we implemented a physical scaling factor that could be applied to the resultant 3D chromosome model, for which we used the mean radius calculated from 3D chromosome painting experiments. We implemented the visualization of the 3D chromosome structures using the python Mayavi package for 3D visualization.

Gene expression and GSEA

Gene expression analysis used published microarray data sets of Lackner et al. (33) available from the EBI ArrayExpress database (E-MTAB-2086). The probes were mapped using the hugene10sttranscriptcluster.db annotation package available in Bioconductor. The microarray data were then quantile-normalized using robust multichip analysis (RMA) to produce log2 expression scores, and differential expression analysis was performed using a linear model fit of the normalized data in the limma Bioconductor package (66). Compartment switching gene sets G1 to G6 were defined using transcription start sites from the UCSC RefSeq annotation (hg19 version update 6 October 2015). Genes on chromosomes not defined in the Hi-C experiments, including chromosome Y (LF1 cells are female) and unassigned contigs, were not considered. Genes were assigned to one of the gene sets using the 40-kb binned PC1 eigenvector values for the three conditions that were overlapped with the gene transcription start site. If one of the conditions had a first PC1 eigenvector value of zero, the gene was assigned to an undefined category. For 389 genes, the presence of multiple distinct isoforms resulted in gene placement in more than one category. GSEA (67) was then conducted on the Lackner et al. (33) data set using the gene sets G1 to G6 defined in Fig. 3D and table S4. GSEA was run using the java command line version of the preranked algorithm with 2000 permutations. To overlap compartment switching genes with molecular phenotypes characteristic of senescent cells, we used the following published senescence-associated gene annotations: the DDR response (37), cell proliferation (38), KEGG cell cycle, SASP (38, 39), genes whose expression changes in cells that have up-regulated p16 (38, 40), and genes involved in the regulation of chromatin (38).

Statistical analysis

Statistical significance of large sample size experiments including the SVL ratio calculations, DNA-FISH, 3D chromosome painting, and comet assays was assessed using nonparametric hypothesis testing with the Mann-Whitney U test. For the FAIRE yield assay, which had a low sample size (n = 3), significance was assessed using a parametric hypothesis test with the Welch’s t test. GSEA statistical significance was assessed using GSEA software that calculated FDR. To compare the positions of TAD boundary positions (fig. S7), Pearson’s correlation was performed. All statistical tests were implemented using the R programming language.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/2/e1500882/DC1

Table S1. Sequencing and iterative alignment statistics for Hi-C experiments.

Table S2. Resampling statistics for Hi-C experiments.

Table S3. Overlap of switching TADs with ENCODE regulatory elements.

Table S4. Compartment switching gene sets corresponding to Fig. 3D.

Table S5. GSEA summary statistics for compartment switching gene sets.

Table S6. Overlap of compartment switching genes with gene sets of senescence-associated phenotypes.

Table S7. FISH probes used for 3D DNA-FISH experiments.

Table S8. Summary statistics for 3D DNA-FISH experiments.

Table S9. Sanger sequencing validation of quiescent and senescent Hi-C libraries.

Fig. S1. Hi-C interaction matrices for the q arm of chromosome 2.

Fig. S2. Hi-C interaction matrices for the q arm of chromosome 3.

Fig. S3. Hi-C interaction matrices for the q arm of chromosome 4.

Fig. S4. Hi-C interaction matrices for the p arm of chromosome 4.

Fig. S5. Genomic feature analysis of contact probability.

Fig. S6. Comparison of first and second Hi-C experiments.

Fig. S7. Characteristics of TADs and A and B compartments.

Fig. S8. Representative genes that switch compartments.

Fig. S9. Physical distances between individual loci within a single chromosome arm.

Fig. S10. Quantification of comet assay images.

Fig. S11. Measurement of chromosome arm volumes.

Fig. S12. Measurement of centromere and telomere volumes in senescent cells.

Fig. S13. Comparison of Hi-C data between replicative senescence and oncogene-induced senescence.

Fig. S14. High-resolution comparison of Hi-C data between replicative senescence and oncogene-induced senescence.

Movie S1. Rotating movie of the 3D Hi-C model for chromosome 18 in quiescent (left structure) and senescent cells (right structure).

Movie S2. Rotating movie of the 3D Hi-C model for chromosome 4 quiescent (left structure) and senescent cells (right structure).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We would like to thank the Brown University Genomics Core Facility, Leduc Bioimaging Facility, and Center for Computation and Visualization for providing assistance. Funding: This work was supported in part by the following NIH grants: K25 AG028753 and K25 AG028753-03S1 to N.N., R37 AG016694 to J.M.S., R56 AG050582-01 to N.N. and J.M.S., F31AG050365 to S.W.C., and K01AG039410 to J.A.K. S.W.C. was also supported by the NIH Institutional Research Training Grant T32 GM007601. B.S. was supported in part by a Brown University Undergraduate Teaching and Research Award. The Brown Genomics Core Facility was supported by NIH grant P30 GM103410. Author contributions: Conceived the study: N.N. and J.M.S. Designed the experiments: N.N., J.M.S., S.W.C., M.D.C., and J.A.K. Performed the experiments and analysis: N.N., S.W.C., M.D.C., B.S., Y.Z., and J.A.K. Contributed to the supervision: N.N. and J.M.S. Wrote the manuscript: N.N., J.M.S., and S.W.C. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Hi-C sequencing data have been deposited under accession number SRP055421 in the National Center for Biotechnology Information short read sequencing archive.
View Abstract

Navigate This Article