Nucleolus association of chromosomal domains is largely maintained in cellular senescence despite massive nuclear reorganisation

Mammalian chromosomes are organized in structural and functional domains of 0.1–10 Mb, which are characterized by high self-association frequencies in the nuclear space and different contact probabilities with nuclear sub-compartments. They exhibit distinct chromatin modification patterns, gene expression levels and replication timing. Recently, nucleolus-associated chromosomal domains (NADs) have been discovered, yet their precise genomic organization and dynamics are still largely unknown. Here, we use nucleolus genomics and single-cell experiments to address these questions in human embryonic fibroblasts during replicative senescence. Genome-wide mapping reveals 1,646 NADs in proliferating cells, which cover about 38% of the annotated human genome. They are mainly heterochromatic and correlate with late replicating loci. Using Hi-C data analysis, we show that interactions of NADs dominate interphase chromosome contacts in the 10–50 Mb distance range. Interestingly, only minute changes in nucleolar association are observed upon senescence. These spatial rearrangements in subdomains smaller than 100 kb are accompanied with local transcriptional changes. In contrast, large centromeric and pericentromeric satellite repeat clusters extensively dissociate from nucleoli in senescent cells. Accordingly, H3K9me3-marked heterochromatin gets remodelled at the perinucleolar space as revealed by immunofluorescence analyses. Collectively, this study identifies connections between the nucleolus, 3D genome structure, and cellular aging at the level of interphase chromosome organization.


Introduction
The spatiotemporal regulation of genomes correlates with transcription, replication, recombination and repair. In the post-genomic era a new model of human genome organization emerged, which is largely based on high-throughput genomics analyses. A key concept of this model is that chromosomal domains, megabase-ranged functional units of the chromosomes, represent an essential operational level of genome regulation (reviewed in [1][2][3][4][5][6][7][8]). However, several questions remained open about the domain organization of human chromosomes. These include the dynamics and functional consequences of chromosomal domain interactions with nuclear sub-compartments in different cell types and under various physiological conditions. The nucleolus is a paramount example for the functional organization of the genome in space and time. Active nucleolar organizer regions (NORs) from different chromosomes build in all cell types a microscopically visible, dynamic nuclear compartment after mitosis, which disassembles again during the next cell division. Consequently, nucleolusassociated chromatin [9] also undergoes cyclic changes in proliferating cells. While the involvement of the nucleolus in ribosome biogenesis, furthermore in facilitating cell cycle progression, stress sensing and RNP function is well characterized [10,11], we just begin to uncover the molecular characteristics of the nucleolus-associated chromatin and to understand the role of the nucleolus in genome organization and function [12][13][14]. Nucleolusassociated chromosomal domains (NADs) represent the mappable genomic fraction of the nucleolus-associated chromatin. NADs were first identified in HeLa cervical carcinoma and HT1080 fibrosarcoma cells using the combination of high-throughput genomics and immuno-FISH analyses [15,16]. These studies provided us with a snapshot of global genome organization in and around the nucleolus and revealed that NADs represent mainly, but not exclusively a specific heterochromatin compartment. In addition, the mechanistic role of several cis-and trans-acting factors in chromosomal domain-nucleolus interactions has also been initially addressed [17][18][19][20]. However, comprehensive high-resolution maps of NADs and the characterization of the nucleolus-associated genome of normal diploid cells have not been determined. Moreover, the involvement of NADs in the formation of interphase chromosome structure and its dynamics during diverse cellular processes such as cellular aging, differentiation or cell cycle, remained largely unknown.
Cellular senescence is a stable arrest of cell proliferation, whose common form, replicative senescence, is induced by telomere attrition and chromosomal instability. In human and other multicellular organisms, cellular senescence plays a pivotal role in several physiological processes, namely tumour suppression, tissue repair, embryonic development and organismal aging [21,22]. Notably, senescence-related, genome-wide reprogramming of gene expression is accompanied by massive structural reorganization of chromatin [23][24][25][26][27][28]. Changes in the spatial organization of chromosomes include the formation of senescence-associated distension of satellites (SADS), an early and consistent marker of cellular senescence [29], which can be followed by the development of senescence-associated heterochromatin foci (SAHF) under certain conditions [30,31]. Remarkably, SAHF-like chromosome condensation can be induced also in a senescence-independent manner [32]. The nucleolar hallmarks of senescent cells include increased size and fusion of nucleoli in mammals, furthermore instability of ribosomal DNA (rDNA) in yeast (reviewed in [33,34]). Notably, the instability of the rDNA cluster is considered as a key control element of genome maintenance and inducer of cellular senescence in yeast [35,36]. While recent works highlighted the involvement of various histone modifications, lamins and DNA methylation in shaping the epigenetic landscape during cellular senescence in mammals (reviewed in [37][38][39][40]), the molecular details of senescence-related nucleolar remodelling are still elusive.
Here we present high-resolution NAD maps of young and replicative senescent IMR90 human embryonic fibroblasts. Our findings identify a central role for the nucleolus in the nuclear arrangement of specific transcriptionally silent chromosomal domains in proliferating and senescent cells, and suggest molecular mechanisms involved in the regulation of nucleolus-associated chromatin and spatial genome organization.

Mapping of NADs in primary human cells
In two independent experiments, nucleoli of each young and proliferating IMR90 fibroblasts were isolated by biophysical disruption of cells and subsequent fractionation by differential centrifugation. To assess the quality of the purification, the enrichment of the nucleolar transcription factor UBTF and ribosomal DNA (rDNA) were monitored in immunoblot and quantitative PCR experiments, respectively, and the nucleolar fraction was examined also by microscopy (S1 Fig). The DNA of isolated nucleoli was extracted and subjected to comparative genomic hybridization experiments on whole genome tiling microarrays to identify the nonrepetitive DNA content of the nucleolus-associated genome and to address its chromosomal organization. The bimodal nature of hybridization signals provided the basis for genome-wide mapping of NADs with the use of a two-state hidden Markov model (HMM) analysis. By using this approach 1,646 autosomal NADs of young, proliferating IMR90 cells were discovered ( Fig 1A, S2 Fig, and S1 Table). We refer here to this comprehensive list of IMR90 NADs simply as 'NADs', and to the previously identified HeLa NADs [15], which were mapped using an 85%-threshold-based method, as 'top NADs'. NADs cover 1.2 billion bp, approximately 38% of the annotated human genome, and their median sequence length (361 kb) resides in the typical size range of the higher order structural and functional chromosomal domains of mammalian genomes (Fig 1B).
In order to validate the nucleolus association of selected chromosomal regions, 3D immuno-FISH experiments were performed. Nucleoli were labelled by immunofluorescence staining of nucleophosmin (NPM1) and different genomic regions by fluorescence hybridization. The frequency of nucleolus interaction was determined for several genomic regions and plotted against the NAD score value, which was calculated from the microarray data by averaging the log-enrichment values of two replicate experiments. Additionally, our previously collected 3D immuno-FISH data [15] was also integrated in the analysis (Fig 1C and 1D, and S2 Table). The results of population-based and single-cell analyses correlate well, as chromosomal regions with high NAD scores showed typically more frequent association with nucleoli in 3D immuno-FISH than regions with low NAD scores.

Comparative epigenomics of NADs and inter-NAD regions (iNADs) reveals specific heterochromatic features of NADs
In order to uncover the chromatin features of NADs at the genome level, the 15-state chromatin HMM (ChromHMM) map of IMR90 cells [41] was quantitatively analysed. The size of NAD and iNAD genomic regions occupied by the various chromatin states was identified and illustrated on a bar graph. Next, the values of each chromatin state were calculated for the entire population of NADs and iNADs, and their log2 ratios were plotted (Fig 2A). The results provide compelling evidence that NADs are depleted in active chromatin features and they can be mainly characterized by the 'heterochromatin' and 'quiescent/low' ChromHMM states. The main features of these chromatin states are high levels of DNA methylation, low DNaseI accessibility, and low incidence of genes, which are mainly repressed. Further investigations of the overall DNaseI accessibility and transcriptional activity in NADs and iNADs revealed large differences with clearly lower values in NADs (Fig 2B and 2C), supporting the results of the ChromHMM analysis.
To identify the replication timing profile around NAD borders, NADs were aligned at their 5' and 3' ends (p to q chromosome orientation), and normalized Repli-seq signals were averaged in 500 kb windows spanning the borders. Replication signals were measured recently at five different time points of the S/G2 phases of the cell cycle [42], and their distributions appear as five waves around NAD borders (Fig 2D). The waves show a regular spatiotemporal pattern within NADs, in which early-to-late replication timing correlates with the distance to  (13, 14, 15, 21 and 22), centromeres and some pericentromeric regions were not analysed because they are not present in the current human genome assembly. (B) Histogram of NAD sizes. Median = 361kb, a total of 1,646 NADs were identified. (C) 3D immuno-FISH analysis of NAD and inter-NAD regions (iNADs) in IMR90 cells. Nucleolus association of a chromosomal domain is illustrated by showing the Z-projection of an IMR90 nucleus on the left and the corresponding single light optical sections with the associated and non-associated allele on the right. BAC hybridization signals are shown in green, nucleolar staining in red and DAPI counterstain in blue (scale bar: 1.6 μm). (D) Hybridization signals (percentage of nucleolus-associated alleles) are plotted against the according microarray signals (average log2-fold difference of the nucleolar signal over the background). Red and grey circles indicate genomic regions that reside in NADs and iNADs, respectively (see S2  Table for   Bar graphs show total and relative amounts of ChromHMM states in NADs and iNADs on the left and right, respectively. The chromatin states and their colour code correspond to the Primary Core Marks segmentation the border. Interestingly, only the depletion of the latest replicating fraction (G2) is notable outside of NADs in the 0 to -500kb distance range, whereas the signal intensities of the other four fractions (S1 to S4) get indistinguishable with increasing distance. Altogether, the result is in good agreement with the late-replicating nature of heterochromatin, the characteristic chromatin state of NADs.
Comparative genomics of NADs, iNADs, lamina-associated domains (LADs) and inter-LADs (iLADs) uncovers NAD-specific genomic features As described above, several lines of experimental evidence point to a considerable enrichment of heterochromatin in NADs. The constitutive heterochromatin of cultured mammalian cells predominantly localizes to the nuclear periphery, pericentromeric bodies and the perinucleolar region, and these nuclear sub-compartments are to large extent functionally overlapping (recently reviewed in [43]). Genomic analyses of the nuclear periphery identified the maps of lamina-associated domains (LADs) in proliferating Tig3 human fibroblasts by using the DamID method [44], and later also in IMR90 cells by a Lamin B1 ChIP-seq approach [24,25]. To statistically evaluate the similarity of LADs and NADs of human embryonic fibroblasts, the three LAD datasets were compared to the NAD list of IMR90 cells and to each other in a quantitative manner. As illustrated by the large intersections of the Venn diagrams, the IMR90 NADs show a substantial overlap with the LADs of human diploid fibroblasts ( Fig 3A). However, about one-third of the total NAD-and LAD-covered genomic regions are non-overlapping and chromosome-specific differences in the patterns are clearly recognizable (S3 Fig). To identify differences between genomic features of NADs and inter-NAD regions (iNADs), furthermore LADs and inter-LADs (iLADs) of the annotated human genome iNADs, the incidence of various sequence features was addressed in these chromosomal domains of human primary fibroblasts (Fig 3B-3E). The results revealed that defensin (DEF) and olfactory receptor (OR) genes are enriched in otherwise gene-poor NADs, whereas they are depleted in generich iNADs and show intermediate levels in LADs and iLADs. Although zinc finger (ZNF) genes are depleted in NADs, the extent of their depletion is less than the overall depletion of protein-coding genes in NADs. A similar ZNF gene distribution can be observed between LADs and iLADs ( Fig 3C). Non-coding RNA genes are also less prevalent in NADs and LADs than in iNADs and iLADs, however 5S RNA genes are slightly enriched in NADs and LADs ( Fig 3D). Repetitive DNA analyses showed that satellite repeats and LTR elements are enriched in NADs, and to a far lesser extent also low complexity repeats, simple repeats, LINEs and DNA repeats. In contrast, SINEs are depleted in NADs due to the lower incidence of Alu repeats in these genomic regions. Similar distributions of repetitive elements can be observed between LADs and iLADs with the clear exception of satellite repeats, which are more frequent Average percentage-normalized signals of different replication domains within a distance of 500 kb from aligned 5' and 3' NAD borders are shown. Repli-Seq signals from five different time points of the cell cycle (S1-S4 and G2) were taken from Pope et al., [42] and averaged over 1 kb windows. Only NADs with a width >500 kb were considered (n = 652). The arrows drawn at the bottom show away from the aligned NAD/iNAD borders.
in iLADs than in LADs (Fig 3E). Taken together, satellite repeats, furthermore OR and DEF gene-coding regions represent the most specific, characteristic genomic elements of NADs. While the prevalence of specific sequence features in NADs is similar to that of top NADs [15], the enrichment of individual features is less pronounced.
Weak nucleolus association is characteristic for escaper regions on the X chromosome As IMR90 cells possess one active and one inactive X chromosome having substantially different chromatin conformations and no Y chromosome, the sex chromosomes were not included in the aforementioned genome-scale NAD analyses. In an extended mapping of NADs a few regions of the Y chromosome emerged as nucleolus associated, due to cross-hybridization or false assignment on the microarrays. In contrast, almost the entire X chromosome appeared as NAD in this analysis. The heterochromatic nature of the inactive X chromosome and the frequent association of the entire chromosome with the perinucleolar heterochromatin might explain that irrespective from the sub-nuclear localization of the active X, most regions of the X chromosome are classified as NAD (S4 Fig). Nevertheless, certain regions of the X chromosome were identified as iNADs. Since some of them clearly overlapped with the pseudo-autosomal regions that escape X inactivation (Fig 4A), the correlation between escaper genes and iNADs of the X chromosome was further addressed. The 'escaper', 'heterogeneous' and 'inactive' states of genes were classified according to Carrel and Willard ([45] and S3 Table) and NAD score values were computed for all three groups. A boxplot representation of the results demonstrates that the escaper group is less frequently associated with nucleoli than the two others ( Fig 4B). The visualization of microarray data along with escaper gene positions provides a more detailed picture. Besides the fact that nearly all iNAD regions contain escaper genes, the remaining escapers also show a strong coincidence with weaker nucleolus association within NADs ( Fig 4A). These observations strongly suggest that local spatial distension, for instance looping of active genes leads to their efficient separation from the more compact heterochromatin of the inactive X chromosome during nucleolus isolation.
Integration of NAD and Hi-C data reveals that NAD-NAD interactions dominate in the 10-50 Mb distance range over iNAD-iNAD and heterotypic NAD-iNAD interactions Chromosome Conformation Capture (3C) provides a way to divide chromosomes into domains by measuring contact probabilities between chromosomal segments. Hi-C analyses with increasing resolution led to the genome-wide determination of A and B compartments [46], topological domains [47], and finally contact domains [48] of human chromosomes. However, this methodology does not deliver information about the nuclear position of chromosomal domains and their interaction with nuclear bodies. To integrate complementary information about the spatial organization of the human genome in IMR90 cells, a comparative genomic analysis was performed using the NAD data and the to date highest resolution Hi-C datasets of this cell line [48]. Venn diagrams show the overlap between NADs and the six genomic subcompartments, which were determined by Hi-C-based segregation of the human genome ( Fig 5A). First, the results demonstrate that there is only little overlap between NADs and the euchromatic, gene dense A1-A2 subcompartments. Second, nearly half of the B1 subcompartment, which represents facultative heterochromatin, coincides with NADs. Third, approximately three-quarter of the B2 and two-third of the B3 subcompartments, primarily constitutive heterochromatin, correspond to NADs. Importantly, these are the largest Hi-C subcompartments. Of note the strongest overlap between B2 and the NADs, which was also detected with the top NADs [48], and also the increased B3 -NAD overlap. Finally, NADs coincide to a greater extent with B4 than with A subcompartments. The interpretation of this overlap is ambiguous due to the minute size of B4. B4 contains KRAB-ZNF gene clusters and displays a specific chromatin pattern, and it is the far smallest Hi-C subcompartment. Altogether, the results of the comparative analysis of Hi-C subcompartments and NADs are largely consistent with the results of the genomic analyses shown in Figs 2 and 3, and uncover that 74% of NADs reside in B2/B3-type constitutive heterochromatic chromosomal regions.
In the nucleus of mammalian cells the chromosomes occupy distinct territories [1,49,50], and their ultrastructural organization is influenced by intra-chromosomal association of chromosomal domains. To compare the contribution of NADs and iNADs to intra-chromosomal interactions at different size scales, frequencies of NAD-NAD, iNAD-iNAD and mixed NAD-iNAD contacts were calculated for each chromosome from IMR90 Hi-C datasets and plotted against linear sequence distance ( Fig 5B and S5 Fig). Intra-and interarm interactions were visualized separately to consider a possible influence of centromeres on intrachromosomal contact probabilities. Visual inspection of the plots revealed that NAD-NAD contacts dominate at the 10-50 Mb distance range over iNAD-iNAD and mixed contacts. Over larger distances iNAD-iNAD contacts became the most frequent, whereas NAD-iNAD contact probabilities are clearly the lowest ones at all distance ranges. Moreover, a strong increase of contact frequencies can be observed in the <10 Mb distance range, which corresponds well to the size range of topologically associating domains (TADs). TADs are considered as the basic units of chromosome folding and they can be defined by measuring genomic interaction frequency changes along the chromosomes: the interaction frequencies of two loci are high within TADs but sharply drop at the boundary between neighbour TADs [51]. In order to visualize NAD-NAD and iNAD-iNAD contacts at >10Mb distances in a locus-specific manner, the Hi-C contact frequencies of individual chromosomes were displayed also on heat maps (S6 Fig). To aid the identification of homotypic contacts, NADs or iNADs were masked in additional individual heat maps. Strikingly, NAD/iNAD and Hi-C-contact-based TAD segmentations of the chromosomes revealed highly similar patterns. In addition, the heat maps showed again that NAD-NAD contacts are enriched at 10-50 Mb and depleted at >50 Mb distances.
Next we calculated the average interaction intensity within each NAD-NAD, iNAD-iNAD, and NAD-iNAD contact to address the involvement of NADs in chromosomal contacts at different distance ranges. To exclude intra-NAD contacts from the analyses, contact frequencies at <10Mb distances were omitted. The contact frequencies of different contact classes were calculated at two different distance ranges, at 10-50Mb versus >100Mb and shown in a boxplot diagram (Fig 5C). Due to the <100Mb size of human chromosomes 16 to 22 and the fact that the p arms of the acrocentric chromosomes 13 to 15 are not annotated, only chromosomes 1 to 12 were investigated here. The boxplots show that homotypic interactions of chromosomal domains are the prominent ones in the organization of chromosome territories at >10Mb distances. Importantly, the results also reveal that NAD-NAD contacts are the dominant ones at the 10-50Mb distance range. Next, the contribution of the different contact classes to the overall contact frequencies within and between chromosome arms was visualized (S7 Fig). It is clearly recognizable that NAD-NAD interactions are the dominant ones within intraarm contacts over 10Mb distances, whereas NAD-NAD and iNAD-iNAD interactions show similar frequencies within interarm contacts. The prevalence of NAD-NAD interactions in the intraarm contacts compared to interarm contacts suggests that the centromere may reduce interarm contact probabilities, however the linear distance seems to play the major role in determining intrachromosomal NAD-NAD contact probabilities.

Replicative senescence causes fusion of nucleoli, but only local, transcription-dependent changes in NADs
Next, we addressed the questions how nucleolar morphology, NAD organization and global gene expression alter during cellular aging. IMR90 fibroblasts were grown under standard cell culture conditions and the aging of the cell population was monitored (S8 Fig). In order to evaluate senescence-associated changes in nucleolar number and morphology, proliferating and senescent IMR90 cells were stained for NPM1 by immunofluorescence. Counting of NPM1 signals showed that proliferating cells exhibit typically 2-4 nucleoli and senescent cells a single, large nucleolus ( Fig 6A).
In order to discover genome-wide changes of nucleolus-associated chromatin at high resolution, NADs of senescent cells were mapped and compared to NADs of young, proliferating cells. Surprisingly, the two NAD maps were highly similar. Altogether 39 Mb sequence was specific for young NADs and 70 Mb for senescent NADs, which respectively represent 3.3% and 5.9% of the total NAD size (1.2 billion bp) determined in young IMR90 fibroblasts (S4 Table). Importantly, chromosomal fragments that can be determined as young-only or senescent-only NAD regions have a median size of less than 20 kb. Accordingly, most of the changes involved just parts of individual NADs (Fig 6B and S9 Fig). Systematic analysis of sequence features in young-only and senescent-only NAD regions revealed that they are particularly enriched in protein-coding genes (Fig 6C and S4 Table). To identify whether the changes in nucleolus association correlate with transcriptional dynamics gene expression microarray experiments were performed. A targeted bioinformatics analysis revealed that the loss of nucleolus association correlates with an overall significant increase in transcript levels, while the gain of nucleolus association with a significant decrease of it ( Fig 6D). Additionally, hierarchical clustering analysis of gene expression microarrays was performed using the data from this and two other recent studies [25,52], in which the same experimental system was investigated. The result showed co-clustering of young and senescent samples providing an additional, robust quality control of the proliferating and senescent states (S10 Fig). Nucleolus-Satellite repeat cluster interactions are impaired in senescent cells, whereas rDNA and telomeres do not display remarkable changes in nucleolus association frequency To uncover the role of cellular aging on the nuclear position of selected genomic loci 3D immuno-FISH analyses were performed. The comparison of the association frequencies of five chromosomal loci with nucleoli showed only little differences between young and senescent cells (Fig 7). This result is in agreement with the largely unaltered maps of NADs in senescence. In the next experiments the spatial dynamics of tandem repeat arrays was addressed, which build the core of the nucleolus-associated chromatin, but are not present on NAD maps because they reside in non-annotated genomic regions. The clusters of rDNA and the distal junctions of NORs displayed basically no difference between young and senescent cells. There was no indication for complete inactivation of NORs, that is, complete separation of strong hybridization signals from the nucleolus. The distal junctions appeared as discrete loci at the nucleolar periphery, whereas rDNA signals were distributed within the nucleoli. Similar to these loci, the global distribution of telomeres did not show detectable aging-related difference at the nucleolar and peri-nucleolar space, however the signal intensities were severely reduced in replicative senescent cells indicating telomere attrition. In contrast, senescence-associated distension of satellite repeats led to their substantial disassociation from nucleoli (Fig 7). Collectively, the 3D immuno-FISH and genomics analyses showed that mainly satellite repeat arrays and individual genes within NADs are the subjects of nucleolus-associated chromatin reorganization in senescence.
The senescence-dependent spatial reorganization of H3K9me3-marked chromatin differs at the perinucleolar space and the nuclear periphery Next we decided to uncover senescence-dependent alterations of H3K9me3-marked constitutive heterochromatin in our model system for the following reasons: First, the role of H3K9 methylation has been demonstrated in the organization of nucleolus-associated chromatin in Drosophila [53]. Second, while genomic redistributions of H3K4me3 and H3K27me3 histone modifications, markers for active and facultative heterochromatin respectively, were investigated in replicative senescence in IMR90 cells [26], H3K9me3-marked chromatin dynamics was analysed only in oncogene-induced senescence [24,54]. Third, H3K9me3-marked centromeric and pericentromeric satellite repeat clusters are essential components of nucleolus-associated chromatin [15], and they are subject to senescence-associated distension (see Swanson et al. [29] and Fig 7). Semi-quantitative immunoblot analyses were performed to monitor the overall level of H3K9me3, and the results showed that there is a strong decrease in senescent cells compared to young, proliferating ones. However, this change goes hand in hand with a somewhat less decrease in global histone H3 level. In addition, robust decrease in Lamin B1 and no detectable alteration in Lamin A/C levels were observed (Fig 8A and S11 Fig). Next, the co-localization of the strongest H3K9me3 signals with pericentromeric HSATII repeat clusters was demonstrated in 3D immuno-FISH experiments in IMR90 cells (Fig 8B). In order to reveal senescence-associated changes in the nuclear distribution of constitutive heterochromatin, quantitative immunofluorescence experiments were performed [55]. The analysis focused on the perinucleolar space and the nuclear periphery (Fig 8C), the two main sites of heterochromatin accumulation in the human nucleus. First, the fluorescence intensity of the H3K9me3 staining was measured in the respective areas and calculated as per cent of the total fluorescence intensity. The results revealed that H3K9me3 signals decrease remarkably, from 5.2% to 2.2%, at the perinucleolar space in replicative senescence (Fig 8D). In contrast, there was only a little difference between the relative fluorescence intensities of proliferating and senescent cells at the nuclear periphery (9.5% vs. 8.4%). Second, the coefficient of H3K9me3 signal variation (C.V. = standard deviation/mean of fluorescence intensity) was calculated, which gives an indication of the heterogeneity of the staining in the respective areas (Fig 8E). The results showed a more heterogeneous staining in the nucleus of senescent IMR90 cells with 0.677 vs. 0.573 C.V. values in senescent and proliferating cells, respectively. Interestingly, at the perinucleolar space the C.V. values did not differ significantly in proliferating and senescent cells (0.576 compared to 0.555), but a significant increase was detected at the periphery in senescent compared to proliferating cells (0.709 vs. 0.632). Last, the abundance of the most heterochromatic regions in the perinucleolar space and at the nuclear periphery was evaluated. Therefore, the image containing the H3K9me3 signals was segmented by thresholding the fluorescence signals at 90% of the maximal intensity. While the abundance of the 10% brightest pixels, which represent the most heterochromatic regions, did not change significantly at the nuclear periphery (15.7% in proliferating vs. 13.7% in senescent cells), it strongly decreased at the perinucleolar space (5.6% in proliferating vs. 1.3% in senescent cells) (Fig 8F). We suppose that the global loss of H3K9me3 at the perinucleolar space (Fig 8D) is mainly due to the senescence-associated distension of satellites, which results also in the selective loss of intensely stained areas (Fig 8F). In contrast, the signal heterogeneity remains largely unaltered at the perinucleolar space (Fig 8E), which could be explained with the stable nucleolus association of NADs in senescence (Fig 6 and S9 Fig). Considering that the H3K9me3 signal distribution gets heterogeneous in other nuclear regions (Fig 8E), we conclude that the nucleolus may play a functionally important role in maintaining the 3D organization of constitutive heterochromatin in replicative senescence.
In order to collect initial insights about the epigenetic regulatory network underlying the senescence-dependent alterations of nuclear architecture, gene set enrichment analyses (GSEA) were performed (S12 and S13 Figs). Comprehensive gene sets were compiled that include genes encoding ATPase and regulatory subunits of chromatin remodelling complexes (CRC), histone chaperones regulating chromatin assembly and disassembly, DNA modifying enzymes and proteins that bind modified DNA, furthermore histone modifying enzymes (S5 Table). The results revealed a global reduction of the mRNA levels of epigenetic regulators in senescence, which correlates well with the non-proliferative status of the cells. However, the transcript levels of a small number of genes, including ERCC6, BAZ2A, KAT2B, SETDB2, EZH1, KDM3A and KDM7A, were significantly increased. This suggests that their gene products could be actively involved in genome regulation in replicative senescence. Notably, we also observed significant up-or down-regulation of various histone modifying enzymes that use H3K9 as substrate (S12 Fig). The result suggests that this residue may act as target point of different senescence-dependent cellular signals. To validate the results of gene expression array experiments, quantitative RT-PCR was performed (S12 Fig and S6 Table). The levels of nineteen mRNA species, which showed significant alterations in the microarray analysis, were measured in proliferating and senescent cells, and two control mRNAs (TFRC and SMARCA1) were used for normalization. Seventeen targets (ERCC6, BAZ2A, CHD4, ATRX, CHD7, BAZ1B, NPM1, NCL, NPM3, CHAF1A, DNMT1, DNMT3B, PCAF, HDAC2, PRMT5, SUV39H1, KDM2B) were selected based on their known nucleolar function [56]. The remaining two, PSG4 and LMNB1, were selected because they show strong up-and down-regulation in senescence, respectively. The qRT-PCR results largely supported the results of the microarray analysis, as only the results for ATRX showed clear difference in the two independent assays. Taken together, the gene expression analyses lay the foundation for the identification of chromatin regulatory mechanisms that are involved in senescencedependent nuclear remodelling.

Discussion
In this study we used cell fractionation and genomics analyses to discover the features of NADs in proliferating and replicative senescent human embryonic fibroblasts and to untangle intensity projections of Z-sections of representative confocal microscopy images are shown to aid the visualization of all repeat clusters within a single nucleus. Nucleolar staining is in red, DAPI counterstain in blue and FISH signals are in green on all images (scale bars: 1.6 μm). From top to bottom: The distal junction regions of NORs and rDNA repeat clusters remain associated with nucleoli in senescent cells. The signal intensity of telomeres is reduced in senescence, but their overall association with nucleoli does not show specific changes. HSATII and alpha-satellite repeat clusters display senescence-associated distension and reduced nucleolus association.
https://doi.org/10.1371/journal.pone.0178821.g007 Our results demonstrate that NADs are predominantly heterochromatic chromosomal domains, which remain largely associated with nucleoli also in senescent cells. We show that, despite microscopically observable nuclear and nucleolar re-organization, only sub-NADs change their association with nucleoli during cellular aging. These changes clearly correlate with altered transcriptional activity. Our data suggests that transcriptional activitydependent disassociation is characteristic also for X escaper genes [45] located on the nucleolus-associated inactive X chromosome [57,58]. The largest senescence-dependent rearrangements of the nucleolus-associated chromatin are, however, not visible on the genomic maps of NADs, because they are caused by the distension of satellite repeat clusters that are not present in genome assemblies. We show altered nucleolar association of pericentromeric and centromeric satellite repeats, which corresponds to the SADS phenotype. Notably, SADS was shown to coincide with increased transcription of satellite repeats [59]. To summarize these findings we propose a model in which active NORs, bearing active and potentially also inactive rDNA repeats [60] and additional NOR-specific sequence elements around the rDNA clusters [61], safeguard the maintenance of the NAD network in replicative senescence. At the same time, transcriptionally activated regions of specific NADs and satellite repeat clusters get disassociated from nucleoli and consequently alter the composition of nucleolus-associated chromatin. Regarding the nuclear remodelling, we demonstrate global loss and local spatial rearrangements of H3K9me3-marked heterochromatin in replicative senescence. We screened for potential epigenetic regulators of senescence-dependent nuclear reorganization also by using a GSEA approach and observed several H3K9 modifying enzymes on the list of genes with significantly altered expression level. These results support each other, and suggest a role for this histone modification in the senescence-dependent dynamics of nucleolus-associated chromatin. Remarkably, we found also BAZ2A among the handful of chromatin regulator genes that show elevated expression level in senescence. The BAZ2A protein is thought to be involved in silencing and constitutive heterochromatin formation at rDNA, centromeres and telomeres [62,63], and it has been recently reported that depletion of BAZ2A promotes escape from senescence [64]. This suggests a role for BAZ2A in the regulation of a subset of H3K9me3-marked chromatin domains in senescence.
Although senescence-dependent chromatin alterations have been investigated in recent studies [24][25][26][27][28]54], several aspects of the nuclear changes remained unclear. One of them, the dynamics of nucleolus-associated chromatin, has been addressed here. Still, the picture is far from complete, and differences in the experimental systems make direct comparisons between on a light optical section of a representative confocal microscopy image. The lamina-and nucleolusassociated areas label 240 nm distances from the edges of the DAPI and nucleolus staining, respectively. (D) Bee swarm plots of relative fluorescence intensities show senescence-dependent small decrease in H3K9me3 levels at the nuclear periphery and strong reduction at the perinucleolar space. Proliferating and senescent IMR90 cells were stained for H3K9me3 and the relative immunofluorescence intensities were measured at the nuclear periphery (lamina) and at the perinucleolar space ( NADs in cycling and senescent human cells the different studies often difficult. For instance, SAHF formation was not prominent in our replicative senescence experimental setup, which distinguishes this study from comprehensive chromatin and nuclear architecture analyses of oncogene-induced senescence (OIS) with a characteristic SAHF phenotype [24,27,54]. Nevertheless, the changes observed in constitutive heterochromatin rearrangements in replicative senescence in our work are partially consistent with the results of OIS and other replicative senescence studies. Most importantly, the alterations of H3K9me3 patterns at the nuclear periphery are reminiscent of its previously reported reorganization in LADs [24]. This reorganization accompanies senescent-dependent Lamin B1 depletion, which was demonstrated also here. Notably, we find about 70% overlap between the genomic maps of NADs and LADs in human fibroblasts, which is in agreement with the stochastic association of several chromosomal domains with either the lamina or the nucleolus [65,66], and in general with the probabilistic model of the spatial organization of the human genome [3]. Indeed, when nucleoli are in close proximity to the lamina, several chromosomal domains may also be simultaneously associated with both nuclear compartments. We predict that NADs on chromosomes bearing an active NOR show higher contact frequencies with the nucleolus than with the lamina. Our chromosome-based NAD-LAD comparison, and a pioneering single-cell genomic analysis of LADs in KBM7 cells seem to support this idea [67]. However, this issue should be addressed by single-cell nucleolus genomics in future studies, in which the dissimilar senescence-dependent dynamics of LADs and NADs could also be more precisely addressed. We hypothesize that globally lowered contact frequencies with the nuclear periphery may accompany with increased nucleolar contact frequencies. This scenario would be consistent with the senescence-dependent destruction of the lamina and the simultaneous enlargement of nucleoli. Notably, the latter one correlates with increased ribosomal RNA precursor levels, which is due to delayed processing and acts as a senescence inducer [68]. In a next step, in nuclei undergoing SAHF formation the association probability of chromosomal domains with the nucleolus might also be reduced.
Concerning the senescence-independent role of NADs in chromosome organization, our bioinformatics analyses of high-resolution IMR90 Hi-C data [48] revealed that NAD-NAD interactions are the dominant ones in the 10-50 Mb genomic distance range and in intra-chromosome-arm contacts, whereas iNAD-iNAD interactions are more frequent over larger genomic distances. The high interaction frequency of NADs in the 10-50 Mb distance range correlates well with the current view that a compact, transcriptionally inactive nuclear compartment builds the core of chromatin domain clusters [69]. Importantly, the pattern of NADs and iNADs largely resembles that of the TADs, which connects the spatial organization of the interphase genome to the nucleolus, the centre of ribosome biogenesis. As ribosomal RNA synthesis determines nucleolar assembly and represents the primary response site of cell growth regulation, we speculate that active NORs may physically link the cell's metabolic activity to 3D genome organization. Notably, the spatial organization of the genome can be considered as moderator of chromosomal communication [70]. The presented study and several previous observations about the role of the nucleolus in shaping genome architecture from yeast to human (reviewed in [12][13][14]71,72]) lead us to postulate that the nucleolus is an important factor in moderating this communication.

Cell culture and nucleolus isolation
Human IMR90 embryonic fibroblasts were obtained from Coriell Repositories (Cat. No. I90-79) and cultivated in DMEM (Gibco Cat. No. 21885-025 supplemented with 10% v/v Foetal Calf Serum, 100 U/mL Penicillin, 100 μg/mL Streptomycin) at 37˚C in humidified, 5% CO 2 atmosphere and regularly tested for mycoplasma contamination. Proliferating young and senescent cells were cross-linked with 1% formaldehyde and nucleoli were isolated as described [73]. Nucleolus-associated DNA was prepared from two independent experiments, both from young and senescent cells, for subsequent microarray analysis. Total genomic DNA of the same IMR90 cells served as control for one young and one senescent sample during comparative hybridization, whereas the other young and senescent samples were hybridized against nonnucleolar DNA. This DNA was collected during the nucleolus preparation and only the DNA of the last, nucleolar fraction was excluded from it. Since the depletion of nucleolus-associated DNA in the control did not alter markedly the NAD patterns, the results obtained from the two independent experiments were combined both for the young and senescent samples (S2 Fig). Quality controls of nucleolus preparations were performed as described [15].

Cellular senescence assays
Cell populations were kept in culture for two weeks after they stopped growing. During this time the growth medium was replaced every second day. The senescence status of the IMR90 cell populations was monitored by senescence-associated beta-galactosidase staining of fixed cells (BioVision Senescence Detecion Kit), and immunofluorescence staining of the MKI67 proliferation marker protein in fixed cells by using a rabbit polyclonal antibody (Santa Cruz sc-15402).

Computational analyses
If not stated otherwise, all computational analyses were performed in R/bioconductor using default parameters (R-project.org/bioconductor.org, version 3.3; R Core Team (2016). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/). The hg19 version of the human genome served as reference. We used hg19, because most of the comparative genomics analyses were performed with datasets available in this format.

Gene expression microarray experiments
Gene expression microarray analyses were performed using TRIzol-extracted total RNA from young, proliferating and senescent IMR90 human diploid fibroblasts. Two biological replicate experiments were performed, and the Affymetrix Human Gene 1.0 ST microarray platform was used. Labelling of the samples and hybridizations were carried out at Source BioScience. We calculated expression values using the Robust Multichip Average (RMA) algorithm. Many-to-one probesets to gene relationships were resolved by retaining the probeset with the highest variance across all arrays. In order to compare expression levels with the ones of

Gene set enrichment analysis
Gene sets of chromatin remodelling enzymes, histone chaperones, DNA modification enzymes and binding proteins, furthermore histone modifying enzymes were compiled based on literature search. Gene Set Enrichment Analyses (GSEA) were performed according to the instructions described on the GSEA homepage http://software.broadinstitute.org/gsea/index. jsp and in [74].

Quantitative RT-PCR
To validate the results of the gene expression microarray experiments, RNA of young and senescent cells was prepared from three additional, independent experiments. 200 ng RNA of each sample was reverse transcribed by using 200U of MMLV-RT and random hexamers. The resulting cDNA was amplified in real-time quantitative PCR experiments using SybrGreenintercalation-based quantification. Primer pairs used for amplification were selected using the qPrimerDepot (http://primerdepot.nci.nih.gov). Data were collected with a Rotor-Gene Q system (Qiagen) and analysed using the comparative quantitation module of the system software. The mean and standard deviation values are derived from three independent experiments analysed in triplicate quantitative PCR reactions.

Mapping and genomics of NADs
Raw tiling array signals were subjected to quantile normalization. Nucleolar enrichment was defined as the log2-fold difference of the nucleolar signal over the background (genomic input or supernatant, respectively). Enrichment signals were smoothed per sample by sliding medians in 100 kb windows as in our previous study [15]. Subsequently, enrichments were averaged for proliferating and senescent cells, respectively. Given the bimodal nature of signal distribution, a two-state hidden Markov model (library 'tileHMM') was employed to classify nucleolus-associated domains with a minimum length of 10 kb (S2 Fig). The average enrichment values within each domain served as 'NAD-scores'. Genomic features of NADs and iNADs were determined by using the UCSC Table Browser (http://genome.ucsc.edu/cgi-bin/ hgTables). Jaccard similarity indices that reflect the ratio of the number of intersecting base pairs between two sets to the number of base pairs in the union of the two sets and thus serve as a similarity statistic for comparing the distribution of regions of different genomic features were calculated with bedtools' 'jaccard' function (version 2.24).

Hi-C data analysis
Raw inter-and intracromosomal contact matrices (10 kb windows, MAPQ> = 30) for the IMR90 map presented in [48] (GSE63525) were normalized using the provided normalization vectors. Each interaction was classified for interaction distance, overlap with NAD or iNAD, and chromosomal arm localization.

Immunofluorescence
Cells were rinsed twice in 1xPBS and fixed in 4% paraformaldehyde in 1xPBS for 10 min. All steps were carried out at room temperature. During the last minute of fixation, few drops of 0.5% Triton X-100 in PBS were added. Cells were then washed three times in 1xPBS/0.01% Triton X-100 for 3 min, followed by 5 min in 1xPBS/0.5% Triton X-100. Finally, cells were washed twice with PBS-T (1xPBS/0.1% Tween20) for each 5 min before antibody staining. Cells were incubated with the primary antibodies (anti-NPM1: Santa Cruz sc-6013R, anti-H3K9me3: CMA318 [76], anti-H3K27me3: Active Motif 61017, anti-LMNB1: Santa Cruz sc-6216) in 4% bovine serum albumin (BSA) in PBS-T for 1h in a humidified chamber, washed in PBS-T three times for 3 min and incubated with fluorescence dye-conjugated secondary antibodies in 4% BSA in PBS-T for 1 h. After washing with PBS-T twice for 5 min, the DNA was counterstained with 50 ng/ml of DAPI in PBS-T for 5 min. Slides were rinsed in PBS-T and mounted in Vectashield (Vector). Middle optical sections of 200 nm z-stacks were selected from 8-bit CLSM (Leica TCS SP8) images based on the quality of nuclear morphology. The xy pixel size was 80.25x80.25 nm. A custom-made ImageJ script was used for semi-automated quantitative analysis of the immunofluorescence signals [55]. The following parameters were measured: i) The percentage of total nuclear fluorescence intensities of H3K9me3 within the lamina-and nucleolus-associated areas (regions of interests-ROIs), which label 240 nm distance from the edge of the DAPI and nucleolus staining, respectively, ii) The heterogeneity of the staining in the ROIs as coefficient of variation (C.V. = standard deviation/mean of fluorescence intensity [77]), iii) The percentage of the 10% brightest nuclear pixels in the two ROIs. Nucleolar number and volume measurements were performed on z-stack images acquired with a Zeiss Axiovert200 microscope, using a 'Plan-Apochromat 63x/1.40 Oil' objective. Whole nuclei were recorded with a z-step size of 500 nm. Nucleoli were counted and assigned to cells manually. Voxel sizes were obtained by calibrating the objective, and values were saved and transferred to FIJI/ImageJ and the volume in μm 3 was calculated by measuring the number of object voxels using the '3D Objects Counter' plug-in. The measurements were performed for both the DAPI and nucleolus staining to determine nuclear and nucleolar volumes, respectively. The resulting values were transferred to an Excel sheet and the average nucleolar volume, the total nucleolar volume per cell and the ratio of nucleolar/nuclear volume per cell was calculated. The results were displayed as bee swarm plots, on which the data points represent individual cells.
Supporting information S1 Fig. Scheme and controls of nucleolus isolation. Nucleolus isolation and quality control experiments were performed as described [15,73]. Differential interference contrast (DIC) micrographs show an intact IMR90 cell ('C') and isolated nucleoli ('No'). The nucleus is marked with a dashed line and blue arrows point to nucleoli. Immunoblots show robust enrichment of the nucleolar transcription factor UBTF and depletion of tubulin in the nucleolar fraction. 0.5% of each fraction of the two-step purification was loaded. 'C'-Input (whole cell extract), 'S1' and 'S2'-supernatants, 'No'-nucleolar fraction. Quantitative PCR measurements indicate the enrichment of ribosomal DNA in nucleolus-associated DNA compared to genomic DNA ('No'/'C'). Error bars represent the standard deviation of two nucleolus isolation experiments from young (Y1 and Y2) and replicative senescent (S1 and S2) IMR90 cells, each of which was analysed in nine technical replicate measurements (triplicate quantitative PCR reactions in three PCR runs). (TIF) were used for the analysis, because they showed the highest similarity to NADs (see Fig 2G). The sizes of overlapping and non-overlapping regions (in 10Mb) are shown above and Jaccard coefficients below the diagrams. Chromosome numbers are shown in bold in the middle of each Venn diagram. (B) The proportions of chromosomal NAD regions that do not overlap with LADs (%NAD-not-LAD = (NAD-LAD size)/ (chromosome size)x100) are shown in the diagram. In the case of the acrocentrics the size of the annotated q-arm was considered as 'chromosome size'. Note, that four of the NOR-bearing, acrocentric chromosomes (13,14,21,22) are NAD-dominated. Since chromosome 15 appears as LAD-dominated in this analysis, we speculate that its NOR contains mostly inactive rDNA repeats in the IMR90 cells used in this study. (C) Integrative Genome Browser snapshots of chromosomes 17 and 22. Chromosomal regions that are covered only by LADs (LAD-NAD) or NADs (NAD-LAD) are shown above the LAD and NAD tracks as indicated.
Chromosome 17 is clearly LAD-dominated, whereas chromosome 22 is NAD-dominated. Collectively, the results shown here highlight the nucleolus-targeting potential of NORs, while showing also differences between NORs of different chromosomes. In addition, HSA19 and HSA4 are also particularly enriched in NADs, which is in good agreement with our previous observations [15]. HSA19, a chromosome localized usually internal in the nucleus, contains a large number NADs which do not appear as LADs. The 'top NADs' of HSA4 include the centromeric and both subtelomeric regions, which suggests tethering or anchoring of the chromosome by these regions to the nucleolus. Overexpression of senescence-associated beta-galactosidase was inspected using an enzyme activity assay as shown next to the growth curve. IMR90 cells were fixed with formaldehyde and beta-galactosidase activity assay was performed. Intense blue staining marks senescent cells. The level of the proliferation marker protein MKI67 was monitored by immunofluorescence staining shown in the bottom panel. Young and senescent IMR90 cells were fixed with formaldehyde and stained for MKI67 and DAPI. Young, proliferating cells show positive MKI67 staining, whereas senescent cells are negative for MKI67. MKI67 staining is labelled in red, and DAPI staining in blue on the merged image.  S5 Table for the full list) reveals frequent decrease (blue bars) in their mRNA levels in senescence. The few genes with increased mRNA levels, which are indicated in red, can be considered as active regulators of senescence. Bar graphs show log2 fold changes in mRNA levels of senescent vs. young cells. Epigenetic regulators with known nucleolar, rDNA-related activity are shown in bold, and regulatory subunits of chromatin remodelling complexes (CRC) in italic. Histone modifying enzymes that act certainly (n = 6) or possibly (n = 5) on H3K9 are labelled with dark and light grey background, respectively. (B) Quantitative RT-PCR validation of selected microarray data. Relative amounts of specific mRNA species in total RNA preparations from young and senescent cells were determined by quantitative RT-PCR with primer pairs listed in S6 Table. The bar graphs depict log2 fold changes in mRNA levels of senescent vs. young cells as determined using two different calibrator mRNA species in qRT-PCR experiments. The log2 fold changes measured in microarray experiments are shown next to the corresponding qRT-PCR data (see legend of the graph). Error bars represent the standard deviation of three independent biological replicate experiments, each of which was analysed in triplicate quantitative PCR reactions.  Table. The genes on the GSEA output images are ordered according to decreasing senescent vs. young (S/Y) expression ratios. High and low mRNA levels are marked with red and blue, respectively. Note the occasional differences between the results of the two studies (e.g. NAP1L3), indicating that the S/Y expression ratios of certain genes might be influenced by factors other than the senescence status. The Enrichment in Phenotype sections of the analysis reports are shown in tabular format on the right bottom side, below the 'DNA MODIFICATION' analysis data. (TIF) S1 Table. Genomic coordinates (hg19 genome build) of IMR90 NADs (this study) and HeLa top NADs [15]. (XLSX) S2 Table. Single cell measurements of nucleolus association frequencies of twelve selected genomic regions support the results of high-throughput genomics experiments. 3D immuno-FISH analyses of NAD and inter-NAD regions (iNADs) were performed in young IMR90 cells in this and our previous study [15], and nucleolus association frequencies (% allele) were determined in 50 and 40-44 cells, respectively. The NAD score values of the regions of interest were calculated from the microarray experiments as average log2-fold difference of the nucleolar signal over the background. BAC names and their genomic coordinates (hg19 genome build) are indicated. (XLSX) S3 Table. Expression features of 630 protein-coding genes on the inactive X chromosome as determined by Carrel and Willard [45]. Genomic coordinates (hg17 and hg19 genome builds) of the genes are also shown.  Table. qRT-PCR primer sequences, genomic locations, GC%, Tm, and estimated cDNA and genomic amplicon sizes. (XLSX)