Two contrasting classes of nucleolus-associated domains in mouse fibroblast heterochromatin

  1. Paul D. Kaufman1
  1. 1Department of Molecular, Cellular and Cancer Biology, University of Massachusetts Medical School, Worcester, Massachusetts 01605, USA;
  2. 2Program in Computational Biology, Fred Hutchinson Cancer Research Center, Seattle, Washington 98109-1024, USA;
  3. 3Department of Molecular, Cellular and Cancer Biology, Program in Bioinformatics and Integrative Biology, Program in Molecular Medicine, University of Massachusetts Medical School, Worcester, Massachusetts 01605, USA
  • Corresponding authors: Julie.zhu{at}umassmed.edu, paul.kaufman1{at}umassmed.edu
  • Abstract

    In interphase eukaryotic cells, almost all heterochromatin is located adjacent to the nucleolus or to the nuclear lamina, thus defining nucleolus-associated domains (NADs) and lamina-associated domains (LADs), respectively. Here, we determined the first genome-scale map of murine NADs in mouse embryonic fibroblasts (MEFs) via deep sequencing of chromatin associated with purified nucleoli. We developed a Bioconductor package called NADfinder and demonstrated that it identifies NADs more accurately than other peak-calling tools, owing to its critical feature of chromosome-level local baseline correction. We detected two distinct classes of NADs. Type I NADs associate frequently with both the nucleolar periphery and the nuclear lamina, and generally display characteristics of constitutive heterochromatin, including late DNA replication, enrichment of H3K9me3, and little gene expression. In contrast, Type II NADs associate with nucleoli but do not overlap with LADs. Type II NADs tend to replicate earlier, display greater gene expression, and are more often enriched in H3K27me3 than Type I NADs. The nucleolar associations of both classes of NADs were confirmed via DNA-FISH, which also detected Type I but not Type II probes enriched at the nuclear lamina. Type II NADs are enriched in distinct gene classes, including factors important for differentiation and development. In keeping with this, we observed that a Type II NAD is developmentally regulated, and present in MEFs but not in undifferentiated embryonic stem (ES) cells.

    As the site of ribosome biogenesis, nucleoli are organized around rDNA repeats (McStay and Grummt 2008; Pederson 2011) and consist of liquid-phase separated and functionally distinct layers of proteins and RNAs (Brangwynne et al. 2011; Feric et al. 2016). Nucleoli are surrounded by perinucleolar chromatin that is the subject of this work. This perinucleolar localization is a functionally relevant example of heterochromatin-mediated silencing, because transgene proximity to mammalian nucleoli is correlated with gene repression (Fedoriw et al. 2012b). To date, the DNA associated with nucleoli has been analyzed in human cells and the plant Arabidopsis thaliana (Németh et al. 2010; van Koningsbruggen et al. 2010; Pontvianne et al. 2016; Dillinger et al. 2017), and these studies have detected specific nucleolus-associated domains (NADs). Consistent with their heterochromatic status, NADs are enriched in the H3K9me3 histone modification characteristic of constitutive heterochromatin, frequently include large multigene clusters, and are correlated with low gene expression and density (Matheson and Kaufman 2016).

    In eukaryotic cell nuclei, genomic loci are distributed in three-dimensional space in a nonrandom manner. In particular, heterochromatic regions frequently reside at the nuclear lamina (NL) and/or in the perinucleolar region (Németh et al. 2010; van Koningsbruggen et al. 2010; Kind et al. 2013, 2015; Ragoczy et al. 2014). A major question in eukaryotic cell biology is how nonrandom localizations of loci are achieved and how they foster functional sequestration or, conversely, how interactions among accessible regions are coregulated (Schoenfelder et al. 2010). Multiple studies suggest that positioning of human cell heterochromatin to either the lamina or perinucleolar region appears to be stochastic and reshuffled at each mitosis when chromosomes are condensed and then decondensed in the subsequent interphase (van Koningsbruggen et al. 2010; Kind et al. 2013). Heterochromatin interactions are also regulated during cellular differentiation (Peric-Hupkes et al. 2010; Meuleman et al. 2013). These long-range chromosome interactions in mammalian cells are of great interest because they can regulate the developmental timing (Vernimmen et al. 2007) or the variegation of gene expression (Noordermeer et al. 2011). However, mechanisms governing these dynamics remain poorly understood.

    The architecture of nucleolar-associated satellite DNA repeats changes during mouse development (Aguirre-Lavin et al. 2012), but systematic characterization of murine NADs had not previously been done. To lay the groundwork for studying the biological roles of NADs in a system that allows analysis of mammalian developmental processes, we have mapped the nucleolus-associated domains (NADs) in mouse embryonic fibroblasts (MEFs). The methods and data presented here demonstrate that in these cells, there is a distinct, previously uncharacterized subset of perinucleolar heterochromatin that is not subject to stochastic laminar positioning. Further, this work provides a platform for future investigation of the dynamics of the nucleolar heterochromatic compartment.

    Results

    Isolation of nucleoli from cross-linked and non-cross-linked MEFs

    For our initial analyses of murine NADs, we compared nucleoli prepared with or without formaldehyde cross-linking (Fig. 1A; Supplemental Fig. S1A). Both methods yielded highly purified nucleoli based on phase microscopy images (Supplemental Fig. S1B; see also the “NAD-seq” QC files at https://data.4dnucleome.org). Quantitative PCR suggested that rDNA sequences were enriched relative to the input total genomic DNA in both the cross-linked and non-cross-linked preparations (Supplemental Fig. S1C). Immunoblot analyses showed enrichment of nucleolar protein fibrillarin relative to lamin A (Supplemental Fig. S1D). Additionally, we characterized the small RNAs present in the non-cross-linked preparation, observing that the nucleoli were highly enriched for the nucleolar U3 RNA (Supplemental Fig. S1E). We have not attempted to reverse cross-links for the recovery small RNAs from cross-linked samples.

    Figure 1.

    Identification of mouse nucleolus-associated domains (NADs). (A) A schematic workflow of nucleolar isolation. Details that differ between the cross-linked and non-cross-linked versions of the preparation are illustrated in Supplemental Figure S1. (B) Analysis of raw NAD sequencing data and comparison to heterochromatin features. All data shown are from MEF cells. The strongly nucleolar-associated Chromosome 18 is shown in its entirety. From the top, tracks shown are DNA replication timing (RT) data (Hiratani et al. 2010) (early replicating regions have positive values in cyan; late replicating regions have negative values in red) and LADs (red) (Peric-Hupkes et al. 2010). Below those are raw sequence data from cross-linked cells (blue, experiment #26): read counts for total genomic DNA (genomic), nucleolar-associated DNA (nucleolar), and the nucleolar/genomic ratio, followed by the analogous data from non-cross-linked (nonXL) experiment #29. (C) Analysis of NAD-seq by multiple bioinformatic packages, including NADfinder. As in B, MEF cell data for Chr 18. Below the RT and LAD data are the nucleolar/genomic ratios from cross-linked (XL) experiment #26 and the peaks called by NADfinder (default cutoff q < 0.05), followed by the analogous data from non-cross-linked (nonXL) experiment #29. Next are shown comparisons of peaks called by the bioinformatic packages MACS2 (broad peaks setting), EDD (50-kb-sized bins), hiddenDomains (50-kb-sized bins), and normR (50-kb-sized bins, v1 < 0.01, q < 10 × 10−3). All software tested except NADfinder called almost the entire chromosome as peaks. NADs called by NADfinder generally correlate with LADs and late RT. (D) As in C, showing Chromosome 5. Note the gradual decrease in the ratio going from left to right (away from the centromere on the left of this acrocentric chromosome). Unlike MACS2 and EDD, NADfinder was able to call peaks on regions distal from the centromere. (E) A shorter chromosomal region, ∼1.3 Mb from Chr 11qA1, analyzed with the same tracks as C and D, plus MEF RNA-seq data at the bottom (GSM2453367, GSM2453368). This region encompasses a gap (red box) between two NADs identified by NADfinder, but was predicted to be part of a single large NAD region by EDD, hiddenDomains, and normR. Note the transcribed genes within the boxed gap. (F) As in C, showing a distal telomere-proximal 2.5-Mb region of Chr 7, including the Kcnq1 locus. This locus has previously been shown to associate with nucleoli (Pandey et al. 2008; Fedoriw et al. 2012a). Of the bioinformatics programs used, only NADfinder correctly called this locus as a peak. (G) As in F, showing a 2-Mb region of Chr XqE. Note the flat level of nucleolar enrichment across this region. The red bar at the top is probe pPK831, a BAC that was observed to be negative for nucleolar association in DNA-FISH experiments (see Fig. 3). This region was correctly called negative by NADfinder, but incorrectly identified as positive by EDD, hiddenDomains, and normR.

    Bioinformatic analysis of nucleolar-associated DNA

    We performed three biological replicate experiments for both the cross-linked and non-cross-linked protocols. For each experiment, we purified DNA from the isolated nucleoli, alongside genomic DNA from a sample of whole, unfractionated cells from the same population for normalization of the nucleolar sequencing data. These DNAs were used to generate sequencing libraries via PCR-free protocols, yielding genome-scale maps. Visual inspection of these “NAD-seq” data showed that the total genomic DNA was largely evenly distributed across the genome (Fig. 1B), occasionally interspersed with short, poorly represented regions (Supplemental Fig. S3A). The nucleolar DNA was distributed with distinct peaks and valleys. The ratio of nucleolar/genomic read counts provides a raw metric of enrichment of nucleolar association.

    We expected mouse nucleolar-associated regions to be highly enriched in heterochromatin, as observed in human cells (Németh et al. 2010; van Koningsbruggen et al. 2010; Dillinger et al. 2017). Therefore, we compared our raw data to data sets that distinguish mouse euchromatin and heterochromatin. For example, late DNA replication (Hiratani et al. 2008, 2010) and association with the NL (Peric-Hupkes et al. 2010) are properties generally associated with heterochromatin (these features are colored red throughout the figures). We noticed that peaks in our raw nucleolar-associated DNA sequence data frequently overlapped these features. Additionally, we observed that the patterns of nucleolar enrichment often appeared similar in each of the triplicate experiments performed with the cross-linked (XL) and non-cross-linked (nonXL) protocols (Supplemental Fig. S2). Consistent with those profiles, Pearson's correlation analysis of the nucleolar/genomic ratios from each of the six experiments displayed a high degree of similarity between cross-linked and non-cross-linked data (Supplemental Fig. S2). Together, these data suggested that both of our biochemical preparation methods reproducibly yielded mouse nucleolar heterochromatin.

    To facilitate comparisons to other genomic features, we then sought to call peaks in our data sets and tested several established bioinformatic packages developed for analysis of LADs or ChIP-seq data (Fig. 1C,D; Supplemental Fig. S3A). We tested MACS2 in the broad mode, which is frequently used for the analysis of ChIP-seq data (Zhang et al. 2008; Feng et al. 2012). Because of the large size of heterochromatin regions like NADs and LADs compared to most transcription factor binding sites, we also tested enriched domain detector (EDD) (Lund et al. 2014), hiddenDomains (Starmer and Magnuson 2016), and normR (Helmuth et al. 2016; Kinkley et al. 2016), packages suited for detecting enrichments over large genomic regions. We tested each of these using at least two different settings (Supplemental Fig. S3A).

    Small chromosomes often contain nucleolar organizing regions (NORs, that is, clusters of rDNA repeats), which tether chromosomes to the nucleoli (Strongin et al. 2014). Small chromosomes without NORs are also generally more frequently associated with nucleoli, as observed in the recent genome-wide SPRITE data, which can detect clustering of DNA loci with the ribosomal RNA (Quinodoz et al. 2018). We observed that small chromosomes were often annotated as almost entirely nucleolus-associated using existing tools (e.g., Chr 18) (Fig. 1C), without distinguishing most peaks and valleys observed in the nucleolar DNA sequencing data. Conversely, we observed that for longer chromosomes, most of these packages identified very few peaks on the chromosome end distal from the centromere (pictured on the right side of these acrocentric mouse chromosomes, e.g., Chr 5 in Fig. 1D). We noted that the nucleolar/genomic ratios displayed a notable slope across large chromosomes, with higher average ratio values near the centromere-containing end of the acrocentric mouse chromosomes. This likely reflects the frequent nucleolar association of pericentric heterochromatin (Ragoczy et al. 2014). However, even in regions far from the centromeres, there are still local peaks with smaller nucleolar/genomic ratios and that are surrounded by deeper valleys. These local peaks often coincide with heterochromatic features. These observations suggested that local baseline correction would improve the accuracy of peak identification in regions that are not near centromeres, while also facilitating the discrimination between peaks and valleys on small chromosomes.

    Therefore, we developed a new tool to analyze our data, a Bioconductor package termed NADfinder that includes local baseline correction (Supplemental Fig. S4A–C; Methods). We note that local baseline correction was required for the NADfinder program to identify peaks distal to the centromere. That is, without local baseline correction, and similar to other tools (e.g., Fig. 1D), only a few strong peaks near the centromere-proximal end of the chromosome were identified (Supplemental Fig. S4B,C, green tracks). In contrast, inclusion of the baseline correction step generated peak profiles that are distributed across the chromosomes (Supplemental Fig. S4B,C, blue tracks). We observed that these baseline-corrected peaks frequently overlap with lamin-associated regions and late replicating regions (Supplemental Fig. S4B,C, red), consistent with previous results suggesting that NADs are largely heterochromatic. In addition to the chromosome-wide views shown in Figure 1B–D, we also analyzed the peaks called by NADfinder on smaller scales. First, we performed genome-wide analysis of the NAD boundaries determined by NADfinder. The analysis results show a sharp upward transition from called non-NADs to NADs in the aggregated plot of the read coverage (Supplemental Fig. S3E). The sharp transition extends over 10 kb centered at just outside of the peak boundaries and indicates a good correspondence between the signal and the peak boundaries. Next, we paid particular attention to the gaps between the peaks identified by NADfinder (e.g., the boxed regions in Fig. 1E; Supplemental Fig. S3B–D), in cases when these fall inside large contiguous peaks identified by the other software (Fig. 1C). We observed that these gaps coincided with diminution of the nucleolar DNA signal. It also appeared that these gaps frequently overlapped with highly expressed genes, suggesting that these gaps are indeed non-heterochromatic regions. This idea will be explored more thoroughly in Figure 4 (see below).

    We also illustrated two examples in which NADfinder was important for proper peak calling. First, we note that NADfinder, but not the other software examined, identified a peak containing the known nucleolar-associated locus Kcnq1 (Fig. 1F; Pandey et al. 2008; Fedoriw et al. 2012a). Second, a locus that was not nucleolar associated in DNA-FISH experiments (Fig. 3, see below) was placed in peaks by several programs, but not by NADfinder (Fig. 1G). Together, these observations suggested that NADfinder would be particularly useful for analyzing our data.

    In the first biological replicate experiment for each purification method, we sequenced approximately 200 million reads. To determine whether this amount was necessary, we performed preliminary peak calling using NADfinder and then performed subsampling analyses (Supplemental Fig. S4D,E). Subsampling analysis results showed that the percentage of peaks detected and overlapped with those from the whole data set initially increases quickly with increased depth of sequencing and plateaus at a subsampling depth of ∼25% (i.e., 50 million reads). We therefore performed 50 million reads for samples in the subsequent second and third replicate experiments for both experimental protocols: cross-linked (XL) and non-cross-linked (nonXL). We updated the NADfinder package with functionality for analyzing data sets with biological replicates and for generating P-values adjusted for multiple hypothesis testing. We used this feature to generate the preliminary NADfinder peak tracks shown in Figure 1C,D and Supplemental Figure S3A.

    Distinct subsets of nucleolus-associated domains

    In addition to MEF LADs (Guelen et al. 2008; Peric-Hupkes et al. 2010), we analyzed the “constitutive interLADs” (ciLADs), which are sequences that are not lamin-associated at multiple stages during neural differentiation of mouse ES cells, first into neural precursor cells and then astrocytes (Peric-Hupkes et al. 2010; Meuleman et al. 2013). That is, ciLADs generally represent euchromatic regions that display much higher levels of gene activity than do LAD regions. To analyze the overlap among NADs and ciLADs, we separated the nucleotides within each NAD peak into subsets of those overlapping with one of the other data sets, or neither. For example, for Chr 17 (Fig. 2A), NADs (blue bars) and LADs (red) displayed large regions of overlap (magenta), in both the XL and nonXL data sets. Analysis of the overlaps in a genome-wide fashion (Fig. 2B) indicated that the majority but not all NAD peak sequences overlap with LADs, regardless of the nucleolus purification method. For example, ciLADs (Fig. 2A, cyan) overlap with a subset of NADs (green). The ciLAD-overlapped subset was smaller than the LAD-associated regions, but occurred in similar regions in samples from the two experimental protocols (XL and nonXL). Thus, we had subdivided NADs into LAD-overlapped and ciLAD-overlapped subsets. For convenience, we will refer to these as Type I and Type II, respectively. Additionally, a subset of NAD sequences was found in neither LADs nor ciLADs (Fig. 2B).

    Figure 2.

    Distinct classes of mouse NADs. (A) Subdivision of MEF NADs based on overlap with LAD or ciLAD regions. Chr 17 is shown. Features associated with heterochromatin (H3K9me2, lamin association, and late RT) are in red, and features associated with euchromatin (early RT and ciLAD regions) are in cyan. At the top is shown DNA RT (Hiratani et al. 2010), followed by H3K9me2 (from MEF cells lacking methyltransferase EHMT2 (also known as G9a), making heterochromatic H3K9me2 more easily visualized) (GSM2803100) (Chen et al. 2018). Next are constitutive LAD and MEF LAD peaks (Peric-Hupkes et al. 2010; Meuleman et al. 2013). Following those are eight tracks of NAD data, the top four from a cross-linked preparation (XL experiment #26), and the bottom four from a non-cross-linked preparation (nonXL experiment #29). At the top of each set is the log2 ratio of nucleolar/genomic sequencing reads (Nol./Genom. Ratio), followed by the peaks determined by NADfinder (blue). These are followed by a pair of tracks in which the NAD peaks were separated at the nucleotide level into regions that overlap the MEF LADs (magenta) and regions that overlap ciLADs (light green). The bottom two tracks are ciLADs in cyan (Peric-Hupkes et al. 2010; Meuleman et al. 2013), and H3K27me3 in olive green (GSM2417097) (Chronis et al. 2017). (B) Venn diagram illustrating the overlap among NADs, LADs, and ciLAD regions. Data from non-cross-linked experiments are on the left, and cross-linked ones are on the right. Numbers show the size of the indicated regions (Mb). (C) Venn diagram illustrating overlap of XL and nonXL NAD peaks. Numbers indicate size of each region in Mb, which are drawn proportional to size. (D) Box plot distributions of NAD lengths, showing median and quartiles in the box.

    We next compared NADs in our XL and nonXL data sets. Genome-scale analysis showed that almost all peaks in the nonXL data set were also found in the XL data set. Additionally, cross-linking resulted in a greater proportion of the genome identified as NAD peaks (Fig. 2C). Specifically, the NAD peaks comprised a total of 1137 Mb in the XL and 811 Mb in the nonXL data sets, representing 41% and 30% of the nonrepetitive mouse genome, respectively. For comparison, LADs comprise ∼40% of the mouse genome (Guelen et al. 2008; Peric-Hupkes et al. 2010), and NADs in human diploid IMR-90 fibroblasts cover 38% of the genome (Dillinger et al. 2017). Also, XL and nonXL murine NADs had similar median lengths (Fig. 2D), comparable to MEF (534 kb) and human fibroblast LADs (553 kb) (Guelen et al. 2008; Peric-Hupkes et al. 2010). Therefore, the overall characteristics of mouse NADs are similar to those of other large heterochromatic domains. However, as shown below, these genome averages mask distinct subsets within the population of mouse NADs.

    Immuno-DNA-FISH validation of genomic data

    To validate peaks called by NADfinder, we performed DNA fluorescent in situ hybridization (DNA-FISH) experiments. We used bacterial artificial chromosome (BAC) probes to detect specific genomic loci, and anti-fibrillarin or anti-nucleophosmin antibodies to visualize nucleoli. We performed at least three biological replicate experiments per probe, scoring the percentage of alleles (n > 100) associated with nucleoli in a population of cells (Fig. 3A–C). First, we chose a negative control BAC (termed pPK871) that was identified by NADfinder as not containing nucleolus-enriched sequences in either the cross-linked or non-cross-linked data (Supplemental Fig. S5A). As expected, we observed a low frequency of overlap between pPK871 and nucleoli (median value of 14% of alleles scored) (Fig. 3D), similar to the values observed for non-nucleolar-associated probes in human cell studies (Németh et al. 2010; van Koningsbruggen et al. 2010; Smith et al. 2014; Dillinger et al. 2017). Then, we analyzed multiple additional probes to characterize NAD associations in MEFs (e.g., Fig. 3B,C; Supplemental Fig. S5) and observed that probes from both Type I and Type II NAD peaks displayed significantly greater frequencies of association with nucleoli than did pPK871. This was evident in the measurement of the percentage of total alleles associated for the individual probes analyzed (Fig. 3D) and when data for the probes of each group were combined (Fig. 3E). Additionally, we compared the observed percentage of nucleolar association of FISH probes to the adjusted P-value (Q-value) of the corresponding peaks calculated by NADfinder (Fig. 3F,G). For both the nonXL and XL data sets, we observed high Q-values for “negative” probes which displayed nucleolar association levels statistically indistinguishable from those of negative control pPK871 (Fig. 3F,G, in blue). In contrast, the Type I (green) and II (red) NAD probes displayed lower Q-values, clearly distinct from the negative probes (note that the x-axis in these graphs is logarithmic). We conclude that NADfinder identified bona fide NADs. Additionally, we compared how NADfinder and other software predicted the nucleolar association of the DNA-FISH probes (Supplemental Table S4). NADfinder outperformed the other software at these specific regions; a genome-scale analysis is presented as part of Figure 4 below.

    Figure 3.

    Three-dimensional (3D) DNA-FISH experiments confirm nucleolar association of NADs. (A) Maximum projection of confocal microscopy images from a 3D immuno-FISH experiment: (red) pPK871 BAC DNA probe; (green) fibrillarin; (blue) DAPI; (scale bar) 10 μm. (B) As in A, with pPK828 as the DNA probe. (C) As in A, with pPK914 as the DNA probe. (D) Graph of percentage of alleles that are associated with nucleolus (mean ± SEM for n ≥ 3 biological replicate experiments). Negative probes (blue) display nucleolar association that is not significantly different than the negative probe pPK871 (P > 0.05, Welch's t-test). In contrast, all Type I (red) and Type II (green) probes display association significantly different from pPK871 (P < 0.03). (E) Data for nucleolar association from D were grouped (see coloration). P-values for the differences between the different groups are shown. (F) Graph comparing the percentage of nucleolar-associated alleles from DNA-FISH (y-axis) versus the adjusted P-value (x-axis) of the nonXL NAD peaks identified by NADfinder that overlap with the indicated BAC probe. Note that the “negative” probes (blue) have an adjusted P-value = 1 = 10 × 100. (G) As in F except that the x-axis is the adjusted P-value of the XL NAD peaks. Note that the “negative” probes (blue) have an adjusted P-value > 10 × 10−3. (H) Single optical sections of confocal microscopy images of a 3D immuno-FISH experiment. pPK834 BAC DNA probe hybridizing to a Type I NAD (green); (red) fibrillarin; (blue) nuclear pore protein NUP62; (scale bar) 10 μm. (I) As in H, with pPK828 hybridizing to a Type II NAD. (J) Graph of percentage of alleles that are associated with nuclear periphery (mean ± SEM for n ≥ 3 biological replicate experiments). Only Type I NAD probe pPK834 (*) displays laminar association that is significantly different than the negative control pPK871 (P = 0.042, Welch's t-test). As expected, three Type II NAD probes (pPK828, 914, and 915) do not display significant laminar association.

    Figure 4.

    Differential gene density, expression, and chromatin state in different classes of NADs. (A) Gene densities of the indicated genomic regions are shown, ranked left to right. (B) Distribution of gene expression levels from RNA-seq data, expressed as log10(FPKM + 1) for the same genomic regions analyzed in A. (C) Jaccard similarity coefficients were clustered to visualize similarities among the indicated genomic regions. Type I NADs cluster with LADs, Type II NADs cluster with ciLADs. (D) As in C, in this case showing comparison of the indicated regions with MEF early and late RT data (Hiratani et al. 2008, 2010). The scale of this figure is different than C because the RT data is from microarray experiments with much sparser representation than deep sequencing data. (E) As in B, in this case including the NAD-splitting regions (NSRs). NSRs are regions between NADs determined by NADfinder, which fall within large contiguous peaks determined by the indicated two other software programs (for examples, see Fig. 1E; Supplemental Fig. S3B–D). That is, NSRs are uniquely identified by NADfinder as non-NADs in these comparisons. NSRs from the XL and nonXL NAD data sets were analyzed separately. In all cases analyzed, NSRs were found to have FPKM distributions most similar to the non-NAD regions of the genome. (F) Venn diagram analysis of the overlap among NADs (nonXL), LADs, and H3K9me3 peaks. Numbers indicate size of each region (Mb), which are drawn proportional to size. (G) As in F, except that the overlaps with LADs and H3K27me3 peaks are shown. (H) As in F, except that the overlaps with ciLADs and H3K9me3 are shown. (I) As in F, except that the overlaps with ciLADs and H3K27me3 are shown.

    Testing association of NADs with the nuclear lamina

    Designation of Type II NADs was based on the overlap with ciLAD regions, which do not associate with the NL in MEFs or during neural differentiation of ES cells (Peric-Hupkes et al. 2010; Meuleman et al. 2013). However, we considered the possibility that some regions could have been classified as Type II NADs if their apparent lack of laminar association was a false negative. We therefore directly tested a subset of our BAC probes for associations with the NL in additional immuno-DNA-FISH experiments. In this case, we tested for colocalization with the edge of DAPI-stained material, which can be highlighted by staining with antibodies against a nuclear pore protein (Fig. 3H,I). We observed that our negative control BAC, pPK871, displayed on average ∼20% association with the lamina, whereas a Type I probe that overlaps with a LAD, pPK834, displays significantly higher association frequencies (Fig. 3J). In contrast, three different Type II probes displayed frequencies of lamina association that were not significantly different from those of the negative control. These data indicate that Type II NADs represent chromatin that preferentially associates with nucleoli rather than the NL in MEF cells.

    Biological properties of the NAD subsets

    To explore NAD biological properties, we first compared gene densities of different NAD subsets to the rest of the genome (Fig. 4A). Consistent with their heterochromatic character, the complete sets of NADs identified by either the XL and nonXL protocols (red bars) display much lower average gene density than the genome average (brown). Furthermore, LAD-overlapping Type I NADs (pink bars), display even lower gene density levels. In contrast, the ciLAD-overlapped Type II NADs (purple bars) display greater gene densities than do Type I NADs, NADs that overlap with neither LADs nor ciLADs, or NADs as a whole. Next, we examined steady-state mRNA levels in the different classes, expressed as FPKM in RNA-seq data (Fig. 4B). We observed that NADs are much less frequent sites of gene expression than the genome-wide average. However, the Type II NADs are distinct in that they display elevated gene expression compared to the other subclasses, consistent with the euchromatic character of the ciLAD subset of the genome (Meuleman et al. 2013). In summary, these data indicate that Type II NADs frequently contain actively expressing genes.

    Our observation of distinct distributions of gene expression levels for the different NAD subsets provided us an opportunity to use genome-wide RNA-seq data to help evaluate the NAD peak calling by NADfinder. Specifically, we explored our observation that the EDD, hiddenDomains, and normR software often called very large, contiguous NAD peaks on small chromosomes, whereas NADfinder split these large regions into smaller peaks that better correlated with the boundaries of lamin B enrichment and late DNA RT (e.g., Fig. 1C). We reasoned that if such NAD-splitting regions (NSRs) between the peaks called by NADfinder were not NADs, they should display a gene expression profile more similar to the non-NAD portion of the genome than to NADs. We identified NSRs (Supplemental Materials; Methods), and found that they are of a similar size range as NADs themselves (Supplemental Fig. S3F, cf. to Fig. 2D). In comparisons to mRNA data sets, we indeed observed that genes in NSRs displayed high gene expression levels similar to the non-NADs; these levels are higher than those observed in NADfinder-identified NADs, including the more highly expressed Type II NAD regions (Fig. 4E). This analysis result demonstrates that NADfinder correctly splits large regions into NADs and non-NADs, whereas other tools aggregated them as contiguous NADs. Therefore, the NSR analysis confirmed our DNA-FISH results on a genome-wide scale; that is, NADfinder distinguished NADs from non-NADs more effectively than other software.

    We next compared LADs and NADs regarding their DNA replication timing (RT), a genomic feature that largely distinguishes euchromatin and heterochromatin (Hiratani et al. 2008, 2010). To do this, we calculated Jaccard similarity coefficients between regions of the genome we designated as NADs and those previously designated as early or late replicating (Hiratani et al. 2008, 2010), or those designated as LADs or ciLADs (Guelen et al. 2008; Peric-Hupkes et al. 2010; Meuleman et al. 2013), and we performed cluster analysis to visualize the similarities. As expected, we observed that LAD regions frequently overlapped with NADs, and more frequently overlapped with late than early replicating regions (Fig. 4C). We also analyzed the Type I and Type II NAD subsets separately. We observed that Type I NADs were more similar to the LADs and NADs as a whole, whereas Type II NADs comprised a distinct subset that clustered most strongly with ciLAD regions (Fig. 4C). Type I NADs were even more heavily enriched in late RT regions than NADs as a whole, whereas Type II NADs were more evenly distributed between early and late RT regions (Fig. 4D). Together, these analyses indicate that Type II NADs comprise a distinct subset of nucleolar-associated chromatin.

    We then analyzed the distribution of heterochromatic histone modifications, namely H3K9me3, a hallmark of constitutive heterochromatin that is generally inaccessible to transcription factors, and H3K27me3, a hallmark of facultative heterochromatin that is frequently regulated in a tissue-specific manner and is accessible to transcription factors and associated with paused RNA polymerase (Becker et al. 2016). We observed that almost half of the NAD regions overlapping with LADs (Type I NADs) also overlapped with H3K9me3 peaks (261 Mb of 542 Mb) (Fig. 4F), whereas a smaller proportion (185 of 542 Mb) overlapped with H3K27me3 peaks (Fig. 4G). These data are consistent with the extensive H3K9me3 marking across LADs (Sadaie et al. 2013), which also feature less broad H3K27me3 peaks at their boundaries (Guelen et al. 2008; Harr et al. 2015). In contrast, the NAD-ciLAD overlapped regions (Type II NADs) displayed different enrichments (Fig. 4H,I). Only 37 Mb of these 134 Mb overlapped H3K9me3 peaks (Fig. 4H), but twice as much (74 Mb) overlapped with H3K27me3 peaks (Fig. 4I). Thus, compared to Type I NADs, Type II NADs lack lamin association, display greater gene density, gene expression, and H3K27me3 marks, and are more often early replicating. Thus, the Type II subset of NADs shares many features with facultative heterochromatin.

    Role of H3K27 methylation

    A hallmark of facultative heterochromatin is methylation of histone H3K27 by the EZH2 subunit of the Polycomb Repressive Complex 2 (PRC2) (Margueron and Reinberg 2011). Previous studies implicated H3K27me3 in the localization of LADs to the nuclear periphery (Harr et al. 2015), so we tested whether the same would be true for nucleolar localization. Using two commercially available inhibitors of the EZH2 methyltransferase, we observed that we could efficiently reduce cellular H3K27me3 levels, detected either via immunofluorescence or immunoblotting (Fig. 5A–E). Both these compounds significantly reduced both the nucleolar and the laminar association of the Type I NAD detected with BAC probe pPK834 (Fig. 5F–J). Analysis of a Type II NAD probe (pPK828) showed a significant loss of nucleolar association in the presence of EPZ6438; the inhibitor GSK126 reduced association but not to a statistically significant level (Fig. 5K–O). As expected, Type II NAD probe pPK828 did not associate with the NL above background levels, and neither drug altered that pattern. Together, these data indicate that H3K27 methylation plays an essential role in the localizations of LADs and NADs to the nuclear and the nucleolar periphery, respectively.

    Figure 5.

    Inhibition of H3K27 methylation blocks NAD localization (A) Immunofluorescence measurement of H3K27me3 levels in control cells. (B) As in A, except cells were treated with GSK126. (C) As in A, except cells were treated with EPZ6438. (D) Quantitation of IF data from AC. GSK126: (****) P < 0.0001, Welch's t-test; EPZ6438: (****) P < 0.0001, Welch's t-test. (E) Immunoblot analyses of cells treated as in AC. (F) 3D immuno-FISH localization of Type I probe pPK834 (green), anti-nucleophosmin antibodies (NPM; ab10530 [Abcam], red), and Anti-Lamin A antibodies (blue). Untreated control cells. (G) As in F, except that cells were treated with GSK126. (H) As in F, except that cells were treated with EPZ6438. (I) Quantitation of pPK834 nucleolar association measured by DNA-FISH as in Figure 3. GSK126: (*) P = 0.040, Welch's t-test; EPZ6438: (*) P = 0.023, Welch's t-test. (J) Quantitation of pPK834 laminar association measured by DNA-FISH as in Figure 3. GSK126: (**) P = 0.0029, Welch's t-test; EPZ6438: (*) P = 0.021, Welch's t-test. (K) 3D immuno-FISH localization of Type II probe pPK828 (green), with other markers as in A. Untreated control cells. (L) As in K, except that cells were treated with GSK126. (M) As in K, except that cells were treated with EPZ6438. (N) Quantitation of pPK828 nucleolar association. GSK126: P = 0.083, Welch's t-test; EPZ6438: P = 0.018, Welch's t-test. (O) Quantitation of pPK828 laminar association. GSK126: P = 0.33, Welch's t-test; EPZ6438: P = 0.63, Welch's t-test.

    Nucleolar associations are sensitive to hexanediol treatment

    Recent data suggest that higher-order chromosomal interactions can be driven by phase separation, involving liquid demixing of heterochromatin components (Larson et al. 2017; Strom et al. 2017). Phase separation also contributes to the self-assembly of nucleoli themselves (Feric et al. 2016). One diagnostic used to study these compartmentalized interactions is the sensitivity of these to the aliphatic alcohol 1,6-hexanediol, which is thought to perturb liquidlike condensates via the disruption of weak hydrophobic interactions (Ribbeck and Görlich 2002). We therefore tested whether features of MEF heterochromatin would also be sensitive to this compound. First, we demonstrated that incubation of cells in media containing 1% hexanediol for 2 h did not affect cell viability (Fig. 6H). We note that the 1% hexanediol concentration tested here is lower than that used to disrupt other liquid demixed entities, such as the nematode synaptonemal complex (Rog et al. 2017), Drosophila HP1 clusters (Strom et al. 2017), yeast or human cell stress granules (Kroschwald et al. 2015, 2017), or nuclear pore complexes (Ribbeck and Görlich 2002). We next tested whether NAD interactions with nucleoli were affected by the treatment with 1% hexanediol for 2 h. We observed different responses depending on which BAC probe was analyzed. The frequency of nucleolar associations of a Type I NAD probe (pPK834) (Fig. 6C) were reduced by the treatment, although this change did not achieve statistical significance (P = 0.11, Welch's t-test) (Fig. 6A). Also, the laminar interactions of this probe were not reduced (Fig. 6B). In contrast, hexanediol significantly reduced the nucleolar associations of two different Type II probes (Fig. 6D–G). For one of the Type II probes (pPK828), we tested for a dose-dependent response and observed a significant loss of nucleolar association of NADs at 1% but not 0.5% hexanediol (Fig. 6D). In sum, these data indicate that different chromatin localizations can have distinct sensitivities to hexanediol.

    Figure 6.

    Response of NADs to hexanediol treatment. (A) Quantitation of Type I NAD probe pPK834 nucleolar associations in control and hexanediol-treated cells measured by 3D-DNA-FISH as in Figure 3. The changes observed did not achieve statistical significance (P = 0.11, Welch's t-test). (B) As in A, laminar associations were measured. The changes observed were not statistically significant (P = 0.24, Welch's t-test). (C) Representative 3D immuno-FISH localization of Type I NAD probe pPK834 (green), anti-nucleophosmin antibodies (NPM, red), and Anti-Lamin A antibodies (blue): (left) untreated control cells; (right) hexanediol-treated cells. (D) Quantitation of Type II NAD probe pPK828 nucleolar associations: (**) P = 0.0090, Welch's t-test. (E) Representative 3D immuno-FISH localization of Type II NAD probe pPK828 (green), anti-fibrillarin antibodies (Fib, red), and anti-NUP62 antibodies (blue): (left) untreated control cells; (right) hexanediol-treated cells. (F) Quantitation of Type II NAD probe pPK915 nucleolar associations: (**) P = 0.0051, Welch's t-test. (G) Representative 3D immuno-FISH localization of Type II probe pPK915 (green), anti-nucleophosmin antibodies (NPM, red), and Anti-Lamin A antibodies (blue): (left) untreated control cells; (right) hexanediol-treated cells. (H) Cells were either untreated or treated with 1% hexanediol for 2 h and stained with trypan blue to measure viability.

    Type II NADs house developmentally regulated genes

    To learn more about the potential roles for Type II NADs in regulated gene expression, we tested for enrichment of Gene Ontology (GO) terms among the genes in these regions (Fig. 7A–C; Supplemental Tables S2, S3). Within the “Biological Processes” classifications, we observed the lowest adjusted P-values are for terms associated with developmental processes, for example organ morphogenesis and sensory organ development. Individual genes within these categories included the developmental transcription factors Pax2 (Fig. 7D) and Sox1 (Supplemental Fig. S6A), and developmental signaling molecules such as Fgf1 (Supplemental Fig. S6B) and Bmp4. Most notable among the “Cell Components” GO terms were keratin intermediate filament proteins, which can occur in gene clusters within a NAD (Supplemental Fig. S6C). In addition, multiple types of ion channels were detected such as Scn sodium channel genes (Supplemental Fig. S6C). The most significant GO-derived “Molecular Functions” classes included transcription factors; other classes include solute transporters such as glutamate transporter Grik3 (Fig. 7E). In summary, these data indicate that a distinct subset of transcription factors related to differentiation, and hallmarks of differentiated cell lineages including ion channels, are enriched within the Type II NADs in MEF cells.

    Figure 7.

    Type II NADs contain developmentally regulated genes and can be differentially localized in MEF and ES cells. (A) GO enrichment analysis of Type II NADs (nonXL). In this panel, the seven most significant Biological Processes annotations and Q-values are shown. (B) As in A, showing the Cell Components annotations. (C) As in A, showing the Molecular Functions annotations. (D) Analysis of specific loci, using the same Integrative Genomics Viewer (IGV) (Thorvaldsdóttir et al. 2013) session as in Figure 2, except that tracks showing NADfinder-derived Z-scores are also included. Genomic region showing two Type II NAD peaks (the one on the right achieved a statistical significance of Q < 10 × 10−3 only in the XL but not nonXL data). The peak on the left contains the Pax2 transcription factor gene. The peak on the right contains the Tlx1 and Lbx1 homeobox transcription factor genes. (E) As in D, showing a Type II peak containing the Grik3 glutamate receptor gene. (F) Single optical sections of confocal microscopy images of 3D FISH analysis of mouse D3 ES cells. BAC DNA probes are labeled in green, fibrillarin in red, and the scale bar is 10 μm. In MEFs, pPK871 is a negative control probe, pPK828 is a Type II NAD, and pPK834 is a Type I NAD (see Fig. 3). (G) Quantification of D3 mES cell DNA-FISH. Graph (left), percentage of alleles from the indicated probes associated with nucleoli (mean ± SEM). Welch's t-tests for difference from pPK871: (**) pPK834 (P = 0.0018); pPK828 (P = 0.62). (H) A diagram that summarizes the distinct characteristics of Type I and II NADs.

    The observation that Type II NADs contain many developmentally regulated genes raised the question of whether there is developmental regulation of perinucleolar heterochromatic associations. As a first test of this idea, we performed DNA-FISH experiments using D3 mouse ES cells (Fig. 7F–H; Hiratani et al. 2010). As we observed in MEF cells, the negative control BAC pPK871 displayed low levels of nucleolar association, and Type I probe pPK834 was significantly more frequently associated. In contrast, a Type II probe (pPK828), which was associated with nucleoli but not lamina in MEF cells (Figs. 2, 3), was not significantly nucleolar associated in D3 cells. These data indicate that NADs can adopt different localizations in different cell types.

    Discussion

    NADfinder, a new Bioconductor package for large domain analysis

    In previous studies of NAD sequences (Németh et al. 2010; van Koningsbruggen et al. 2010; Pontvianne et al. 2016; Dillinger et al. 2017), bioinformatics tools for analyzing these data sets had not been developed. We found that tools developed for analyzing LADs and ChIP-seq data were suboptimal for the analysis of the mouse NADs. Specifically, multiple existing software packages failed to detect NADs in the regions distal to the centromere in long chromosomes and/or would label almost the entire short chromosomes as NADs in our MEF data set (Fig. 1C,D; Supplemental Fig. S3A). We therefore developed a Bioconductor package called NADfinder. We observed that the use of local baseline correction enhanced the recognition of peaks in the distal regions of large chromosomes, while improving the distinction of peaks on small chromosomes (Supplemental Fig. S4). Peaks called by NADfinder were successfully validated via DNA-FISH experiments (Fig. 3). We imagine that NADfinder will not only be useful for NAD studies, but will be globally useful for the analysis of large chromosome domains, particularly in cases in which local background signals will differ significantly across the genome.

    One advantage of deploying our software in the Bioconductor framework is that we can leverage the rich genomic annotation data and borrow functionalities already implemented in other R/Bioconductor packages. For example, we used the ATACseqQC package (Ou et al. 2018) to efficiently read in the alignments, the baseline package for baseline correction (https://cran.r-project.org/web/packages/baseline/index.html), the signal package for smoothing (http://r-forge.r-project.org/projects/signal/), the trackViewer package (Ou and Zhu 2019) for visualizing the signals, the limma package (Ritchie et al. 2015) for statistical analysis of the signal noise ratios across multiple replicates, and the ChIPpeakAnno package (Zhu et al. 2010; Zhu 2013) for annotating NADs with other genomic features of interest and for GO enrichment analyses. NADfinder is freely available (https://bioconductor.org/packages/release/bioc/html/NADfinder.html). Major functions are described in Methods. More detailed information on the commands and instructions for running the software are included in the associated user guide at https://bioconductor.org/packages/release/bioc/vignettes/NADfinder/inst/doc/NADfinder.html and https://bioconductor.org/packages/release/bioc/manuals/NADfinder/man/NADfinder.pdf.

    Diverse classes of perinucleolar heterochromatin

    Previous studies in human cells suggested that there should be extensive overlap among LAD and NAD regions of the genome. Specifically, several experiments revealed extensive reshuffling of the interphase localizations of heterochromatin during each mitosis. For example, in an interphase cell, photoactivation of caged fluorescent histones near the nucleolar periphery results in the labeled molecules being distributed both at the nucleolar periphery and at the lamina in the two daughter cells after mitosis (van Koningsbruggen et al. 2010). Likewise, imaging of a GFP-DamID reader domain protein showed that LADs redistribute not only to lamina but also to nucleoli in daughter cells (Kind et al. 2013). These experiments indicated that the interphase localization of heterochromatin is largely stochastic, with apparent random placement of heterochromatic loci either at the lamina or at the surrounding nucleoli. Therefore, the expectation was that LADs and NADs would be largely the same.

    Previous mapping of nucleolar heterochromatin has been done in human cells and in the plant Arabidopsis (Németh et al. 2010; van Koningsbruggen et al. 2010; Pontvianne et al. 2016; Dillinger et al. 2017). In human cells, the pioneering studies first describing perinucleolar heterochromatin (Németh et al. 2010; van Koningsbruggen et al. 2010) had low sequencing depth by today's standards, but generally supported the idea that NADs are heterochromatic, with reduced gene density and expression relative to the genomic average. More recently, analysis of human IMR-90 fibroblast heterochromatin emphasized the overlap and similarity between NADs and LADs, sharing enrichment for H3K9me3, gene, and repeat density (Dillinger et al. 2017). In summary, the cell biological and genomic experiments mentioned above indicated that the majority of mammalian heterochromatin is directed in a stochastic manner to either the lamina or the perinucleolar region.

    However, there have also been experiments suggesting that there are trans-acting factors that can specifically regulate particular heterochromatin destinations. For example, the lamin B receptor, an integral protein of the inner nuclear membrane, contributes to the laminar association of heterochromatin (Solovei et al. 2013). Conversely, the Chromatin Assembly Factor-1 p150 subunit governs the nucleolar localization of proliferation antigen MKI67 (also known as Ki-67), which in turn governs NAD localization to nucleoli (Smith et al. 2014; Matheson and Kaufman 2017). The specificities observed in these studies suggest that not all heterochromatic loci will localize stochastically, and that some recruitment to specific locations can occur.

    We also note that a recent study showed that an EZH2-containing ribonucleoprotein complex termed MICEE tethers specific loci to the perinucleolar region (Singh et al. 2018). However, EZH2 appears to have multiple roles, because inhibition of H3K27me3 inhibits both laminar (Harr et al. 2015) and nucleolar heterochromatin localizations (Fig. 5). Indeed, the latter data suggest that H3K27me3 could have a crucial and general role in heterochromatin localization, perhaps via contributing to phase separation (Fig. 5). However, there are likely to be multiple mechanisms, because the PRC2 complex that catalyzes H3K27me3 is not required for the nucleolar localization of the Kcnq1 locus (Fedoriw et al. 2012a).

    Facultative heterochromatin is a general term for regions that change their expression status during development. A major hallmark of facultative heterochromatin is H3K27 methylation catalyzed by the EZH subunits of Polycomb Repressive Complex 2 (PRC2). Conversely, constitutive heterochromatin is correlated with lamin association, late RT, and H3K9me3 modification (Trojer and Reinberg 2007; Beisel and Paro 2011; Becker et al. 2016). In these respects, Type II NADs display greater similarity with facultative than constitutive heterochromatin (Figs. 2, 4), and they are enriched in developmentally regulated genes (Fig. 7). We therefore hypothesize that the nucleolar localization of Type II NADs contributes to developmental regulation, and we are now testing these ideas.

    We also note that features of the MEF NADs we have described here may be found in other organisms. For example, to examine conservation from our mouse fibroblast data to humans, we lifted over our cross-linked (XL) NAD peaks from MEF cells to the hg38 human genome using the UCSC Genome Browser liftOver tool with default settings (Kent et al. 2002). We compared our MEF-XL-NAD-hg38 to published data sets for human IMR-90 NADs (Dillinger et al. 2017) which were also obtained from cross-linked samples. We observed that there are indeed many murine NADs that correspond to syntenic human NADs, as shown in a genome-scale Venn diagram in Supplemental Figure S7A. We also tested whether both the Type I and Type II subsets of the MEF NADs were represented among the NADs that were identified in human. Supplemental Figure S7B,C illustrate that both of these subclasses were detected; the distribution of these across two human chromosomes is illustrated; full lists of the lifted-over NADs are in Supplemental Table S5. These observations suggest that indeed there may exist two distinct classes of human NADs, which may have been overlooked, and need further investigation.

    The effect of hexanediol on NADs

    We observed that the nucleolar association by Type II NADs was disrupted by 2 h treatment with 1% hexanediol (Fig. 6). In contrast, the laminar association of a Type I NAD was stable during this treatment; nucleolar association by this probe was reduced without achieving statistical significance. Thus, not all NAD localizations are equally sensitive to this treatment. These data suggest that liquid droplet formation or other self-association phenomena shaping chromatin interactions may have distinct properties at different loci and different nuclear locations.

    Our interpretation of these observations is speculative at this point. Laminar associations could be intrinsically stronger because of factors directly tethering heterochromatin, such as the lamin B receptor (Solovei et al. 2013). However, additional factors likely contribute to the observed effects of hexanediol. For example, inhibition of EZH2 alters NAD association (Fig. 5), consistent with multiple other experiments indicating that altered histone modification can affect domain localization (Harr et al. 2015; Poleshko et al. 2017). We do not yet know the relative importance of altered histone modification levels versus effects of hexanediol on liquid demixing of heterochromatic proteins during our experiments. Of course, those two things could be interrelated. We also do not know the kinetic profile of events during the course of these experiments. Although beyond the scope of the current work, such experiments could be helpful in understanding the mechanistic basis for these observations.

    In summary, our data indicate that a significant subset of perinucleolar chromatin is not frequently distributed to the lamina. Specifically, through locus-unbiased biochemical purification of nucleoli based on their sedimentation properties and resistance to sonication, we found distinct classes of loci in the associated DNA. The majority of NADs identified (Fig. 2B) were the expected loci that overlap with LADs, which we term Type I NADs here. However, we also identify Type II NADs, which were defined as peaks in our data that overlap with regions previously classified as “ciLAD,” that is, regions that do not associate with lamina at multiple stages of cellular differentiation (Peric-Hupkes et al. 2010; Meuleman et al. 2013).

    Recently, other studies have shown that sedimentation behavior of short cross-linked chromatin fragments helps detect functionally distinct subpopulations that cannot be distinguished by histone modification patterns alone (Becker et al. 2017). Therefore, we believe that biochemical studies can generate important data for analyzing the structure and function of mammalian genomes.

    Methods

    MEF nucleoli isolation

    For each preparation, cells were grown in 10 15-cm plates until they reached 80%–90% confluence, for a total of 30–100 × 106 cells. Old media was removed and fresh media was added 1 h prior to nucleoli preparation. An additional plate grown in parallel was reserved for total genomic DNA extraction (Quick-DNA Universal Kit [Zymo Research]). Detailed methods for cell culture and treatments, measurements of cell viability, and 3D DNA FISH, immunocytochemistry, and microscopy are found in the Supplemental Methods.

    Non-cross-linked method of nucleoli isolation

    The method was adapted from the Lamond laboratory protocol (van Koningsbruggen et al. 2010). Each dish was rinsed twice with 5 mL prewarmed PBS, 3 mL 0.25% trypsin-EDTA was added per dish, followed by incubation at 37°C until cells were detached. Ten milliliters of complete media was added to each plate to quench the trypsin, and the cells were centrifuged in 50-mL conical tubes at 925 × rpm (200g) for 5 min (Sorvall Legend RT centrifuge) at 4°C. Cell pellets were washed with PBS, and the packed cell volumes (PCV) were estimated. Five pellet volumes of PBS were added, cells were centrifuged and resuspended in ice-cold PBS, and cell numbers were measured using a hemocytometer. At this point, the cells for the total genomic DNA sample were frozen. For the remainder of the cells, ∼3 × PCV of ice-cold Nuclear Extract Buffer A (10 mM HEPES-KOH, pH 7.9, 10 mM KCl, 1.5 mM MgCl2), plus 0.5 mM DTT was added freshly. Cells were incubated on ice for 10 min to promote swelling in this hypotonic buffer. Swelling was monitored using an EVOS microscope with a 20× objective, using the phase contrast filter. Swollen cells were transferred by pouring into a glass Wheaton Dounce homogenizer (15-mL size, prerinsed with Nuclear Extract Buffer A) with a Type B Pestle, prechilled on ice. Cells were homogenized with 30 slow strokes, and lysis was assessed by examination of a small aliquot on the EVOS microscope. The homogenization was stopped when >90% of the cells were burst, leaving intact nuclei and with various amounts of cytoplasmic material attached. A total extract sample of 0.1–0.5 mL was saved. Homogenized cells were centrifuged at 2000 rpm for 5 min at 4°C. A 0.5-mL aliquot of the supernatant was saved as “cytosolic extract.” The crude nuclear pellet was resuspended with 3 mL solution S1 (0.25 M sucrose, 10 mM MgCl2) and freshly added plus 0.5 mM DTT. Any pellet that could not be resuspended readily by gentle pipetting likely contained lysed nuclei and was discarded. Each resuspended pellet was layered over 4 mL of S2 (0.35 M sucrose, 0.5 mM MgCl2) solution in a 15-mL conical tube and centrifuged at 2000 rpm for 5 min at 4°C to enrich the purity of the nuclei. The purified nuclear pellets were then resuspended with 3 mL of S2 solution by gentle pipetting. A 0.1-mL sample of these resuspended nuclei was reserved as the “Nuclear” fraction. The nuclear suspension was sonicated with 10-sec bursts (with 10-sec intervals between each burst) using a Soniprep 150 (MSE) with a fine probe sonicator and set at power setting 5. The sonication was monitored using a phase contrast microscope. The sample was sonicated until almost no intact cells were detected and released nucleoli were observed as dense, refractile bodies. A 0.7-mL amount of the sonicated sample was layered over 0.7 mL of S3 solution (0.88 M sucrose, 0.5 mM MgCl2) in 1.5-mL microfuge tubes and centrifuged at 20 sec at 16,000g at 4°C. The pellet contains the nucleoli, and the supernatant (combining both the sucrose layer and the upper aqueous material) was retained as the “nucleoplasmic fraction.” The nucleoli were resuspended with 1.0 mL of S2 solution, combined in one microfuge tube, followed by centrifugation at 1430g for 5 min at 4°C. This wash was repeated twice. The pellet contains highly purified nucleoli. The nucleoli were resuspended in 0.2 mL of S2 solution, snap frozen in liquid nitrogen, and stored at –80°C.

    Cross-linked method of nucleoli isolation

    The preparation of formaldehyde cross-linked cell extracts was adapted from previous studies (Sullivan et al. 2001; Németh et al. 2010). Cells were grown as above and fixed with 1% formaldehyde added directly to the media and incubated at RT. After 10 min, the media was removed and the cross-linking was quenched by adding 5 mL 1 M glycine. The cells were washed with ice-cold PBS, collected by scraping in 40 mL PBS, and centrifuged at 200g for 5 min at 4°C. The cell pellet was resuspended in 1 mL high magnesium (HM) buffer (10 mM HEPES-KCl pH7.5, 0.88 M sucrose, 12 mM MgCl2, 1 mM DTT). Cells were then sonicated on ice (2–3 bursts for 10 sec at full power) using a Soniprep 150 (MSE) with a fine probe. The release of nucleoli was monitored microscopically. Nucleoli were pelleted by centrifugation in a microfuge at 15,000g for 20 sec and resuspended in 0.5 mL low magnesium (LM) buffer (10 mM HEPES-KCl pH 7.5, 0.88 M sucrose, 1 mM MgCl2, 1 mM DTT). The sample was subjected to a second round of sonication (1 burst, 10 sec full power) and centrifuged at 15,000g for 20 sec to pellet nucleoli. The pellet was resuspended in 20/2 TE buffer (20 mM Tris-Cl pH 7.5, 2 mM EDTA) for immediate use or in 20/2TE + 50% glycerol to snap freeze in liquid nitrogen and store at −80°C. Detailed methods for the analysis of preparations by quantitative PCR, immunoblotting, small RNA analysis, DNA isolation, and deep sequencing are found in the Supplemental Methods.

    NAD identification and annotation

    To facilitate the analysis of this type of data, we implemented a common workflow for NAD-seq data analysis as a Bioconductor package called NADfinder. We used version 1.6.1 for analyses in this paper.

    To make it easy to use, all the functionalities for analyzing data with multiple biological replicates have been wrapped into four functions: tileCount, log2se, smoothRatiosByChromosome, and callPeaks. An overview of the NADfinder workflow is depicted in Supplemental Figure S4A. Briefly, summarized counts for each moving window (default 50 kb) with a step size (default 1 kb) are computed for nucleolar and genomic samples using the tileCount function. Log2 ratios of counts or CPM are computed between nucleolar and corresponding genomic samples for each window using the log2se function. Chromosome-level baseline correction and smoothing of log2 ratios are performed using the function smoothRatiosByChromosome. Peaks/NADs across biological replicates are called, filtered, and merged using callPeaks. We used the default setting of NADfinder to analyze our data except the following parameter settings: lfc = log2(1.5), countFilter = 200, and cutoffAdjPvalue = 10−3 (10−2 for Chr X, where the signal to noise ratio is less than the other chromosomes).

    Nucleotide-level overlap analysis of NADs and other genomic features such as LADs, H3K9me3 (GSM1621023) (Delbarre et al. 2017), and H3K27me3 (GSM1621022) (Delbarre et al. 2017) were performed using GenomicRanges (Lawrence et al. 2013). Specifically, the functions intersect and setdiff were used to obtain the overlapping nucleotides between two peak ranges, and the nucleotides in only one of the peak ranges, respectively. In addition, GO enrichment analysis of different classes of NADs were performed using ChIPpeakAnno (Zhu et al. 2010; Zhu 2013). FPKM values based on RNA-seq data were from GSM1621026 (Delbarre et al. 2017). Detailed methods for NAD boundary, NAD-splitting region (NSR), and Jaccard analyses are found in the Supplemental Methods.

    Data access

    All sequencing data and quality control information for individual samples generated as part of this study have been submitted to the NIH 4DN Data Portal (https://data.4dnucleome.org) under accession numbers 4DNES15QV1OO, 4DNES8D85R8H, 4DNESR3K7VRN, 4DNES4U7AY3P, 4DNESNEBF9W4, and 4DNES7SVJXAI. See Supplemental Table S1 for details. The NADfinder software is available at https://bioconductor.org/packages/release/bioc/vignettes/NADfinder/inst/doc/NADfinder.html and as Supplemental Code.

    Acknowledgments

    We thank Thoru Pederson and Joan Ritland for review of this manuscript. Research reported in this publication was supported by the National Institutes of Health (National Institute on Drug Abuse, U01 DA040588 to P.D.K.) as part of the 4D Nucleome Consortium.

    Footnotes

    • Received December 2, 2018.
    • Accepted June 10, 2019.

    This article is distributed exclusively by Cold Spring Harbor Laboratory Press for the first six months after the full-issue publication date (see http://genome.cshlp.org/site/misc/terms.xhtml). After six months, it is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/.

    References

    | Table of Contents

    Preprint Server