Introduction

Over 100 years of cytology has recognized characteristic features of nuclear genome organization, including blocks of condensed chromatin enriched at the nuclear and nucleolar peripheries and varying degrees of condensation/decondensation throughout the nuclear interior. Foci of active transcription are concentrated at the edges of chromosome territories and sub-territory condensed chromatin masses and depleted at the nuclear periphery (Belmont, 2022; Bernhard and Granboulan, 1963; Cremer et al., 2015; Cremer and Cremer, 2006; Misteli, 2020; Politz et al., 2013, 2016; Takizawa et al., 2008; van Steensel and Belmont, 2017). Two current models for nuclear genome organization prevail (Belmont, 2022). The radial genome organization model describes a radial gradient of higher gene expression activity towards the nuclear center (Bickmore, 2013; Croft et al., 1999; Girelli et al., 2020; Kolbl et al., 2012; Kupper et al., 2007; Takizawa et al., 2008). The binary model of genome organization describes the differential intranuclear positioning of a peripheral, lamina-associated, late-replicating, transcriptionally inactive Hi-C B compartment versus an interior, early-replicating, transcriptionally active Hi-C A compartment (Guelen et al., 2008; Kind et al., 2013; Lieberman-Aiden et al., 2009; Peric-Hupkes et al., 2010; Ryba et al., 2010; White et al., 2004)

However, neither model acknowledges the considerable cell-type specific variability in nuclear shape and size, as well as the size, number, and relative positioning of nuclear bodies. This variability could dramatically affect chromosome trajectories and the distribution of active and inactive chromatin regions among nuclear locales. Indeed, high-resolution Hi-C has identified two active A subcompartments and three major repressive B subcompartments (Rao et al., 2014). Similarly, TSA-seq, NAD-seq, and SPRITE have distinguished between types of heterochromatin localizing preferentially at nucleoli versus the nuclear lamina and have revealed a type of active chromatin specifically localizing near nuclear speckles (Chen et al., 2018; Nemeth et al., 2010; Quinodoz et al., 2021, 2018; van Koningsbruggen et al., 2010; Vertii et al., 2019; Zhang et al., 2020). However, most previous studies have not attempted to combine multiple genomic readouts across multiple cell types.

Here we integrate light microscopy imaging of major nuclear locales-the nuclear lamina, nucleoli, and nuclear speckles-with DamID (Guelen et al., 2008; van Schaik et al., 2020), TSA-seq (Zhang et al., 2020), and high-resolution Repli-seq (Zhao et al., 2020) genome mapping to investigate how the genome is both spatially and functionally organized relative to nuclear locales across four human cell types hTERT-immortalized human foreskin fibroblasts (HFFc6, abbreviated hereafter as HFF), H1 human embryonic stem cells (hESCs), HCT116 (HCT) colon carcinoma epithelial cells, and K562 erythroleukemia cells.

Results

Cell types differ widely in morphology and arrangement of nuclear locales

Prior analyses have largely overlooked the possible contribution of nuclear and nuclear body morphology to nuclear genome positioning. To measure sizes and shapes of nuclei and nuclear bodies, cell lines were stained for nuclear speckles (hereafter referred to simply as speckles), nucleoli, and the nuclear periphery (NP) with anti-SON, anti-MKI67IP, and anti-glycosylated nucleoporin (anti-RL1) antibodies, respectively, and for DNA using DAPI. Deconvolved wide-field microscopy images (Fig. 1A) were then segmented to generate nuclear 3D models (Fig. 1B, SFig. 1A).

Figure 1:

K562 and H1 nuclei have larger volume (V), lower surface area (SA), and are rounder (lower SA/V) as compared to HCT and HFF nuclei (Fig. 1C). Whereas K562 and H1 speckles distribute radially in 3D, surrounding nucleoli and concentrating within the nuclear interior, HCT and HFF speckles preferentially distribute in the equatorial (center) z-planes while being depleted radially near the nuclear periphery within these equatorial planes (Fig. 1B). H1 and HCT116 have large, often single, nuclear-centrally located nucleoli, while HFF and K562 have smaller, multiple nucleoli distributed between speckles (Fig. 1A, B, D). Total nucleolar volumes per nucleus are largest in HCT116 and lowest in HFF (Fig. 1D).

We compared relative arrangements of these nuclear locales across cell lines by measuring the nearest distances between these locales (SFig. 1B). Rounder nuclei (K562, H1) have significantly larger speckle-to-NP and nucleolus-to-NP distances compared to flatter nuclei (HFF, HCT116) as expected (SFig. 1C, Fig. 1E). But speckle-to-nucleoli distances in K562 are lower than other cell lines, with HFF showing noticeably larger speckle-to-nucleoli distances (SFig. 1C, Fig. 1E). While total speckle and nucleolar volumes scale linearly with nuclear volume in all cell lines, curiously total speckle volume correlates more strongly with nuclear surface area versus volume (SFig. 1D).

Overall, the four cell lines divide approximately into round (K562, H1) nuclei with radial organization of nuclear bodies versus flat nuclei (HCT116, HFF) with nuclear bodies distributed relative to both x-y radial and z-axes. This division is further supported by Principal Component Analysis (PCA) using 14 morphological measures (Fig. 1F, SFig. 1F).

Global differences in nuclear genome organization in different cell types revealed by DamID and TSA-seq

We next surveyed genome organization relative to the NP, using nuclear lamina (lamina) as a target, speckles, and nucleoli by applying LMNB1 (lamina) DamID and LMNB1 (lamina), SON (speckle), and MKI67IP (nucleolus) TSA-seq (Fig. 2A). Speckle TSA-seq is highly conserved across all four cell types (Fig. 2A), as described previously (Zhang et al., 2020) and supported here by the high correlations for speckle TSA-seq comparisons between pairs of cell lines, despite a lower dynamic range in HCT116 and especially HFF (Fig. 2B). In contrast, genome-wide correlations in pair-wise cell line comparisons were substantially lower for lamina DamID and for both lamina and nucleolus TSA-seq as compared to speckle TSA-seq (Fig. 2A-B, SFig. 2A-C). For example, LADs near the ends of long chromosome arms switch from a strong lamina association in HFF and HCT116 to a more interior and nucleolar association in H1 and K562 (Fig. S2A). While in HCT116 cells local maxima in nucleolar TSA-seq align with LADs, in HFF, H1, and K562 cells small peak and “peak-within-valley” nucleolar TSA-seq local maxima instead align with speckle TSA-seq peaks (Fig. S2B).

Figure 2:

These results suggest that genome positioning relative to speckles is highly conserved across cell lines, despite variations in speckle intranuclear positioning (Fig. 1). In contrast, major variations exist across different cell types for genome positioning relative to the lamina and nucleoli.

Gene expression and DNA replication timing exhibit distinct correlations with nuclear locales

We previously identified genomic regions 100 kb or larger, comprising ∼10% of the genome, that showed statistically significant differences in speckle TSA-seq smoothed scores in pairwise comparisons between the same four cell lines (Zhang et al., 2020). Genes within these regions showed a strong positive correlation between differences in their gene expression and speckle positioning. Chromosome regions with different position relative to speckles typically show inverse differences in their position relative to the lamina (SFig. 2D, top panels), though individual genome regions have variable magnitude and even direction of these differences (SFig. 2E). However, the reverse is not true. Pairwise cell-line comparisons reveal ∼40-70% of the genome shows statistically significant differences in lamina TSA-seq over regions 100 kb or larger, with most of these regions showing little or no differences in speckle TSA-seq scores (SFig. 2D-E, bottom panels).

Genes within genomic regions displaying significant differences in speckle TSA-seq in cell line comparisons tend to show increased expression when they are closer to speckles and decreased expression when they are further from speckles (Fig. 2C), as described previously (Zhang et al., 2020). In contrast, genes within genomic regions which in pair-wise comparisons of cell lines show a statistically significant difference in lamina TSA-seq show no obvious trend in their expression differences (Fig. 2C).

Binning all genes as a function of differences in their normalized speckle, nucleolar, and lamina TSA-seq scores show a positive correlation of differences in gene expression with changes in speckle TSA-seq scores but no apparent correlation with differences in lamina or nucleolar TSA-seq scores (Fig. 2D). Linear modeling for correlation between changes in gene expression and changes in TSA-seq reveals a substantially stronger correlation for positioning relative to speckles, as compared to positioning relative to the lamina or nucleoli (Fig. 2E, SFig. 2F); this pattern holds genome-wide (Fig.2E) as well as for correlations over specific genome region types (LAD, iLAD, speckle TSA-seq maxima) (SFig. 2F). Instead, pair-wise cell type comparisons revealed significant correlations between changes in DNA replication timing and changes in TSA-seq with respect to all three locales-speckle, lamina, and nucleoli (Fig. 2F-G).

In summary, whereas gene expression primarily correlates with relative distance to nuclear speckles, DNA replication timing correlates with relative distance to all three nuclear locales, despite the variable positioning of nuclear speckles and nucleoli across cell types.

Both Type I and Type II speckle attachment regions are gene expression hot-zones and DNA replication IZs but show varying gene composition and DNA replication timing

We next sought to determine the relationship between gene expression “hot-zones” (Chen et al., 2018) embedded within larger iLAD/A compartment regions and speckle association. We first identified larger amplitude (Type I) SON TSA-seq local maxima embedded within A1 Hi-C subcompartments versus smaller amplitude (Type II) local maxima embedded within A2 Hi-C subcompartments (red versus blue tick marks, Fig. 3A).

Figure 3:

Type I peaks showed a near unimodal distribution of distance to the nearest speckle with a peak at ∼100 nm distance (Fig. 3B), based on IMR90 fibroblast multiplexed FISH data (Su et al., 2020), used as a proxy for HFF fibroblast measurements. In contrast, Type-II peaks showed a bimodal speckle distance distribution (Hartigan’s dip test (Hartigan and Hartigan, 1985), p-value < 0.05, (SFig. 3A)), with one narrow peak similarly located at ∼100 nm distance and a broader peak centered at ∼500 nm from speckles (Fig. 3B). Notably, Type II peaks show a lower speckle association frequency (<250 nm threshold, Fig. 3B, dashed line) than Type I peaks (Fig. 3C), suggesting that while Type II peaks still associate specifically with speckles their interaction with speckles is weaker and/or less frequent as compared to Type 1 peaks.

Both Type I (red) and Type II (blue) average peak widths are ∼400kb (Fig. 3D); however, Type I peaks have higher speckle TSA-seq and Hi-C compartment scores (Fig. 3D) and contain shorter genes, genes with smaller exon and intron length, and genes with lower fractional intron composition (SFig. 3B). Across the four cell types, Type I peaks also are more conserved than Type II peak (Jaccard indexes of 0.85 versus 0.67 for Type 1 versus Type 2) (SFig. 3C).

Although gene expression levels are similar over Type I and II peaks (Fig. 3D, 2^nd from right), the IMR90 multiplexed image data reveals a trend of gene loci within Type I peaks showing a higher fraction of “ON” states among the cell population as compared with gene loci within Type II peaks (SFig. 3D).

Both Type I and Type II peaks align with DNA replication initiation zones (IZs), identified as peaks in 16-fraction Repli-seq (Fig. 3A). Type I peaks replicate slightly earlier than Type II peaks (Fig. 3E-F, SFig. 3E), although this difference is noticeably larger in HCT116 versus H1 cells (Fig. 3A,E,F). Speckle TSA-seq scores progressively decrease with later replication timing of IZs (Fig.3F). However, DNA replication timing variability (Twidth) is notably larger for Type II peaks, particularly in HCT116 (Fig. 3E, SFig. 3E) suggesting a possible link between variability in speckle association (Fig. 3B) and variability in DNA replication timing (Fig. 3E).

We next used live-cell imaging to corroborate that chromosome regions close to nuclear speckles, primarily Type I peaks, would show the earliest DNA replication timing. Indeed, in HCT116 cells expressing EGFP-SON and mCherry-PCNA, the first PCNA foci appearing within the first 30 mins of S-phase were nearly 100% localized within 500 nm of speckles (SFig. 3F). From 30 to 135 mins after initiation of S-phase, more PCNA foci appeared at increasing distances from speckles (SFig. 3F).

These results show that while Type I and II speckle-associated regions are both DNA replication IZs and show similarly elevated levels of gene expression, they contain different types of genes, show different speckle association frequencies, and show different timing of DNA initiation.

LAD to iLAD conversions mostly correspond to transitions from a repressed/ late-replicating to an intermediate repressed/ late-to-middle DNA replicating chromatin state

The low correlation between differences in gene expression with differences in lamina and/or nucleolar TSA-seq and lamina DamID scores (Fig. 2D&E, SFig. 2D&E) suggested that transitions from facultative LAD (fLAD) to facultative iLAD (fiLAD) might not strictly represent a binary transition from repressed to active chromatin. Analysis across cell lines reveals that lamina DamID scores of a fLAD genomic region can continuously transition from positive peaks in lamina DamID (red rectangles), to peak-within-valley (p-w-v) local maxima (yellow rectangles), or to valleys (v) (blue rectangles) (Fig. 4A, SFig. 4A). A subset of p-w-v fiLADs also showed p-w-v local maxima in their lamina TSA-seq signals together with p-w-v local maxima or small positive peaks in their nucleolus TSA-seq signals (SFig. 4A).

Figure 4:

We implemented an algorithm to define these three domain classes genome-wide based on their local relative increase (”enrichment”) in DamID score (SFig. 4B) (See Methods), revealing that p-w-v fiLADs are more common than v-fiLADs in all cell lines except HFF, which has few fiLADs of either type (Fig. 4B-C). The transitions between LADs, p-w-v fiLADs and v-fiLADs are highly dynamic across the four cell lines (Fig. 4A, SFig. 4C), with some regions transitioning between all three domain types among the four cell lines (Fig. 4A, SFig. 4A & C). We observe an increasing number (Fig. 4B) and total genome coverage (Fig. 4C) of LADs versus p-w-v fiLADs and v-fiLADs comparing H1, K562, HCT116, and HFF cell lines.

p-w-v fiLADs show late to middle DNA replication timing and reduced gene expression levels, rather than the early DNA replication timing and higher gene expression typical of constitutive iLADs (ciLADs) and a subset of v-fiLADs (Fig. 4A-B, SFig. 4D). Moreover, transitions between cell lines from p-w-v-to v-fiLADs show an increase in gene expression of comparable magnitude to that observed in LAD to p-w-v fiLAD transitions (Fig. 4D, SFig. 4E). While p-w-v fiLADs are classified as A compartment, with positive Hi-C compartment scores, in pile-up plots they show lower compartment scores than iLAD or v-fiLAD regions (Fig. 4E). Both p-w-v-fiLADs and v-fiLADs have earlier replication timing than LADs, but p-w-v fiLADs replicate later than v-fiLADs (Fig. 4E-F). Additionally, p-w-v fiLADs typically have an intermediate chromatin state between LADs and v-fiLADs, as defined by their relatively increased levels of “repressive” and decreased levels of “active” histone marks (SFig. 4F).

In summary, p-w-v fiLADs, comprising the majority of fiLADs, retain some heterochromatic marks and repressive functional properties as compared to v-fiLADs, which instead show an earlier-replicating, more active chromatin state comparable to typical iLADs.

p-w-v fiLADs compete with LADs for nuclear lamina association; linking nuclear lamina association with later DNA replication timing but not lower gene expression

To probe the functional consequences of nuclear lamina genome association, we made LMNA, LBR, and double LMNA/LBR K562 knockout (KO) lines. Whereas the LMNA KO shows little change in lamina DamID from K562 wild-type (wt) or the parental clone, the LBR KO and the LMNA/LBR double knockout (DKO) lines show similar but partial reduction in the lamina DamID for most LADs and increases for some iLADs (SFig. 5A). Therefore, we focused our attention on comparisons between wt K562 and the DKO line.

DKO cells showed loss of the nuclear lobulation characteristic of wtK562 nuclei, the appearance of a DAPI-dense perinucleolar rim not apparent in the wtK562 nuclei, and an increased number of DAPI-dense foci in the nuclear interior (Fig. 5A). Based on DamID and TSA-seq, most LADs in DKO cells decrease their interaction and relative proximity with the lamina and increase their proximity to nucleoli but maintain their positioning to speckles (Fig. 5B, SFig. 5B). Unexpectedly, most p-w-v-fiLADs instead increase their lamina association while decreasing their proximity to nucleoli (Fig. 5B, SFig. 5B-C). Indeed, most of these p-w-v fiLADs become actual LADs in the DKO (SFig. 5D); some v-fiLADs also become LADs in the DKO (which may be due to a misclassification of v-versus p-w-v fiLADs) (SFig. 5C-D).

Figure 5:

Scatterplots comparing changes in lamina and nucleolar TSA-seq scores between wt and DKO showed an inverse relationship (Fig. 5C). This shift correlates with ChIP-seq data, showing that regions enriched in H3K27me3 generally shift away from nucleoli towards the lamina in the DKO, while the opposite trend is observed for regions enriched in H3K9me3 (Fig. 5D) (consistent with the enrichment of H3K27me3 in p-w-v fiLADs versus H3K9me3 in LADs (SFig. 4F)).

Despite the wide-spread genomic shifts relative to the lamina and nucleolus, only a small number of genes significantly change their gene expression in the KO and DKO lines (Fig. 5F). Moreover, there is no obvious bias in the direction of gene expression changes as a function of changes in the lamina DamID between the DKO and wt (Fig. 5G). Similarly, there is no obvious bias towards increased gene expression among those genes contained specifically within wt LADs that change their expression and also show reduced lamina interaction (decreased lamina DamID) in the DKO (SFig. 5E). (However, there is an intriguing bias of an increased fraction of differentially expressed genes in LADs versus iLADs in all KO cell lines (SFig 5F-G)).

In contrast, plotting 2-fraction Repli-seq versus changes in DamID shows a consistent trend of slightly later DNA replication timing for regions (primarily p-w-v fiLADs) moving closer to the lamina (Fig. 5B&G). However, LADs that shift further away from the lamina, do not show an obvious progressive shift towards earlier DNA replication timing (Fig. 5B&G, SFig. 5H).

Differential positioning relative to speckles and lamina identifies different types of LAD regions

We next used differences in lamina versus speckle TSA-seq scores to identify heterochromatin regions with differential biochemical and functional properties (Fig. 6, SFig. 6).

Figure 6.

First, whereas an overall inverse relationship between lamina and speckle TSA-seq observed in K562 cells is largely preserved in H1 cells, scatterplots reveal off-diagonal genomic bins (“H1 ROI” orange bins in Fig. 6A, “H1 ROI C1&C2” in SFig. 6A) with disproportionally low lamina TSA-seq scores (Fig. 6A&B, SFig A&B). These genomic bins are largely within several Mbp from centromeres (Fig. 6C, SFig. 6E), are H3K9me3-enriched (Fig. 6D), and are more associated with the nucleolus specifically in H1 cells in which they are less associated with the lamina (SFig. 6B-C). Many genomic regions contained within the H1 ROI-C1 display a similar position in the K562 DKO scatterplot (SFig. 6D), suggesting a LMNA and/or LBR dependence for nuclear lamina association for this subgroup of LADs.

Second, deviations in the inverse relationship between SON and lamin TSA-seq in HCT116 and HFF cells with flat nuclei reveals four types of LAD regions, defined by DamID (SFig. 6F), varying in their histone modifications (Fig. 6E-F): 1) H3K9me3-enriched LAD Cluster C1: low SON, high LMNB1; 2) H3K27me3-enriched LAD Cluster C2: moderate SON, high LMNB12; 3) H3K9me2- and H2A.Z-enriched LAD Cluster C3: low SON, low LMNB1; 4) low in multiple histone marks LAD Cluster C4: moderate SON, low LMNB1. Superimposing these histone mark enrichments over the HCT116 scatterplot more precisely revealed the differential spatial distribution of LAD regions with different types of histone mark enrichments (Fig. 6G).

These four types of LAD regions show different correlations with both gene expression levels and DNA replication timing (Fig. 6H). H3K9me3-enriched C1 LAD regions show the lowest levels of gene expression. H3K27me3-enriched C2, H3K9m2/H2A.Z-enriched C3, and C4 LAD regions show gene expression levels significantly lower than iLADs but slightly higher than C1 LAD regions. DNA replication timing is latest and most uniform for the C1 LAD regions and shows progressively earlier DNA replication timing for the C3, C2, and C4 LAD regions which all still replicate later than iLADs. We defined a “constitutive LAD (cLAD) score” as the number of cell lines out of 7 overall in which a chromosome region is a LAD. cLADs are enriched over LAD regions with the lowest SON TSA-seq scores (C1 and C3) while fLADs are enriched over LAD regions with moderate SON TSA-seq scores (C2 and C4) (Fig 6J).

In summary, LAD regions that segregate differentially relative to speckles and the lamina show different histone mark-enrichments and functional properties.

Polarity of Nuclear Genome Organization

The loss of the strong inverse relationship between speckle and lamina TSA-seq in the HCT116 and HFF cell lines was unexpected and puzzling. The noticeably weaker inverse relationship between speckle and lamina TSA-seq in cells with flat (HCT116, HFF) versus round nuclei (K562, H1) was due specifically to changes in lamina TSA-seq as lamina DamID scores showed less variation across cell lines.

We hypothesized that LADs lying within the nuclear equatorial plane in flat nuclei would have the highest lamina TSA-seq scores: after lamina TSA-staining they would be exposed to a higher concentration of tryamide free-radicals generated and diffusing from the nearby side, top, and bottom of the nuclear lamina as compared to LADs at either the top or bottom of the nuclear lamina. We simulated this effect by convolving the HFF anti-RL1 immunostaining image with a kernel corresponding to the exponential diffusion gradient of tyramide-biotin TSA labeling, yielding a notably elevated TSA signal in the equatorial plane of the nuclear periphery (Fig. 7A).

Figure 7:

Using the IMR90 multiplexed FISH data set [43], we compared the localization in HFF cells of 26 LADs mapping to regions of highest lamina TSA-seq (Cluster 1, Fig. 7B) versus 27 LADs mapping to regions of lowest speckle TSA-seq but with low to moderate lamina TSA-seq (Cluster 2, Fig. 7B). Cluster 1 LADs showed a biased localization in IMR90 fibroblasts towards the equatorial plane as compared to Cluster 2 LADs that distributed more uniformly over the x-y nuclear projection (Fig. 7C). More generally, Cluster 1 and 2 LADs show a similar differential localization in lamina versus SON TSA-seq scatterplots in HFF and HCT116 cells with flat nuclei which is notably different than in K562 and H1 cells with round nuclei (SFig. 7A). Cluster 1 LADs correspond to LADs located towards the ends of long chromosome arms (SFig. 7B). Overall, in fibroblasts there is a trend of distal chromosome arms localizing preferentially towards the nuclear equator versus centromeric regions localizing preferentially towards the nuclear center (Fig. 7D-E and SFig. 7C-D).

We next examined the nuclear polarity of genomic regions associated with different nuclear locales in fibroblasts. Single cell data (SFig. 7E) and average trends over all nuclei (SFig. 7F), show speckle-associated regions concentrated towards the equatorial plane, consistent with the direct localization of speckles within the equatorial plane (Fig. 1B). Nucleolus-associated regions were distributed towards the nuclear center but throughout the z-plane (SFig. 7E-F). Type I speckle attachment regions localize more centrally in the x-y plane and closer to the equatorial plane than Type II speckle attachment regions (Fig. 7F-G, SFig. 7G). The closer lamina proximity for Type II peaks is consistent with the shorter genomic distances of Type II versus I peaks to the nearest LADs ([1] and data not shown). fLADs localize to equatorial plane (Fig. 7G) consistent with the enrichment of H3K27me3 modification in these regions, while cLADs localize furthest away from this plane (Fig. 7G). Additionally, v-fiLADs cluster closer to the equatorial plane and further from the nuclear center as compared to p-w-v fiLADs (SFig. 7I).

To more comprehensively survey how the genome distributes relative to nuclear polarity, we turned to the use of SPIN states, which combine TSA-seq and DamID with Hi-C data to segment the genome into 9-10 states based on their intranuclear localization; SPIN states show differential epigenetic marks, gene expression, and DNA replication timing (Wang, 2021; Wang et al., 2021). SPIN states ordered according to highest to lowest gene expression, and earliest to latest DNA replication timing, showed a pronounced, nearly monotonic increasing distance from the equatorial plane (Fig. 7G) as well as differential localization relative to the x-y nuclear center (SFig. 7H).

Thus, in adherent cells with flat nuclei the distance from the equatorial plane emerges as an additional key axis in nuclear genome organization in addition to the x-y plane radial position.

Discussion

Here we integrated imaging with both spatial (DamID, TSA-seq) and functional (Repli-seq, RNA-seq) genomic readouts across four human cell lines. Our results provide new insights that significantly extend previous nuclear genome organization models, while also demonstrating a cell-type dependent complexity of nuclear genome organization (Fig. 8).

Figure 8.

First, across the four cell lines examined, differences in gene expression correlate primarily with differences in the genome positioning of these genes relative to nuclear speckles (Fig. 8A) rather than distance from the nuclear lamina and despite variations in speckle positioning relative to the nuclear lamina and nucleoli. Consistent with previous observations showing a more interior nuclear localization of genes after their activation (Bickmore, 2013; Takizawa et al., 2008), we demonstrate an inverse trend between changes in speckle and lamina distances for the small fraction of genomic regions that show significant changes in their relative speckle distances. For these regions, there is a strong bias towards increased (decreased) gene expression with decreased (increased) speckle distance. However, several-fold more genomic regions significantly change their relative lamina distances without significantly changing their speckle distances or gene expression. Similarly, LBR/LMNA DKO in K562 cells results in many genomic regions increasing or decreasing their lamina association without changes in their speckle distance or gene expression. In contrast, changes in DNA replication timing correlate with changes in positioning relative to all three nuclear locales-speckles, lamina, and nucleoli (Fig. 8A).

Second, Hi-C A compartment/iLADs are punctuated by strong (Type I) and weaker (Type II) speckle attachment regions which align with DNA replication initiation zones and show elevated gene expression levels above flanking iLAD regions (Fig. 8B). Type I regions contain shorter, more exon-rich genes as compared to Type II regions and are more highly conserved across cell lines than Type II regions. IZs within strong speckle attachment sites replicate earlier than IZs within weak speckle attachment zones with the very earliest DNA replication foci appearing adjacent to nuclear speckles (Fig. 8E).

Third, Hi-C A compartment/iLADs are also punctuated by partially repressed fiLADs which retain weak lamina interaction and, in some cases, elevated nucleolar interactions. Identified as local maxima in lamina TSA-seq and/or DamID, these “p-w-v-fiLADs” comprise ∼2/3 of all fiLADs across the four cell lines examined. P-w-v-fiLADs show gene expression levels intermediate between LADs and typical iLADs while retaining late or shifting to middle DNA replication timing. Only approximately 1/3 of fiLADs, identified as valleys in lamina DamID (v-fiLADs), convert to an active, early-replicating, interior chromatin state. The same fLADs across different cell lines can associate strongly with the lamina as a LAD, retain weak lamina association as a partially repressive p-w-v-fiLAD, or show a true transition to an interior, active, early DNA replicating state (Fig. 8C). Upon LMNA/LBR double KO, p-w-v-fiLADs shift closer to the lamina, suggesting they compete with LADs for lamina association in wt cells (Fig. 8F).

Fourth, we show how LAD regions showing different histone marks-either enriched in H3K9me3, H3K9me2 plus H2A.Z, H3K27me3, or none of these marks-can differentially segregate within nuclei. These results support the previous suggestion of different “flavors” of LAD regions, based on the autonomous targeting of BAC transgenes to the lamina (Bian et al., 2013). Differential nuclear localization also was recently inferred by the appearance of different Hi-C B-subcompartments, which similarly were differentially enriched in either H3K9m3, H3K27me3, or the combination of H3K9me2 and H2A.Z (Spracklin et al., 2023). Here, our use of TSA-seq directly measures and assigns the intranuclear localization of different LAD regions to different nuclear locales.

Fifth, beyond the previously described dependence of genome organization relative to the radial axis in the x-y plane, there exists an unexpectedly high degree of nuclear polarity with respect to an orthogonal z-axis in HCT116 and HFF cells with flat nuclei (Fig. 8D). Strong speckle attachment regions localize closer to the equatorial plane than weak speckle attachment regions. A subset of LADs, primarily at the ends of long chromosome arms, preferentially localize at the lamina within this equatorial plane. Using SPIN states, ordered by DNA replication timing and levels of gene expression, we show that this overall nuclear polarity applies to nearly the entire genome, suggesting a corresponding overall gradient in both later DNA replication timing and lower gene expression with increased distance from the equatorial plane.

In conclusion, by examining genome organization relative to multiple nuclear locales, we were able to demonstrate a higher correlation than previously apparent between nuclear location and functions such as DNA replication timing and gene expression. We anticipate future mapping of the genome relative to additional nuclear locales will reveal a still more deterministic relationship between nuclear positioning and genome function.

Materials and methods

Cell culture

H1, HCT116, and HFF cells were cultured according to standard operation protocols established by the NIH 4D Nucleome Consortium (https://www.4dnucleome.org/cell-lines.html). K562 cells were cultured according to the ENCODE Consortium protocol (http://genome.ucsc.edu/ENCODE/protocols/cell/human/K562_protocol.pdf).

Immunostaining

Adherent cells (HCT116, HFF, and H1) were seeded on 12mm round coverslips and cultured for two days prior to immunostaining. For K562 cells, 150µl of cell suspension was added to coverslips treated with poly-L-lysine as described previously (Chen et al., 2018). All washing steps were 3 ξ 10 mins in Ca, Mg-free PBS at room-temperature (RT), unless specified elsewise. Cells were fixed with freshly made 2% paraformaldehyde (PFA) in PBS for 10 mins at RT and washed. Cells were permeabilized for 5 mins on ice with 1X PBS containing 0.2% Triton X-100 (Sigma-Aldrich T8787) and 1% natural goat serum (NGS) and washed with 1XPBS/1% NGS (PBS/NGS) 3 x 10 mins at RT. Each coverslip was then incubated with 100µl of the primary antibody cocktail in a dark humid chamber at 4°C overnight. The cocktail contained moue-anti-RL1 (Santa Cruz cat# sc-58815, 1:100 dilution), rabbit-MKI67IP (Atlas antibodies cat# HPA035735, 4ng/ml dilution) and rabbit-anti-SON (Pacific Immunology Corp, custom-raised, 4ng/ml dilution) (Chen et al., 2018) diluted in PBS/NGS. The anti-MKI67IP and anti-SON antibodies were directly labelled with CF-594 and CF-640R fluorophores using Mix-n-Stain™ CF® Dye Antibody Labeling Kits (Biotium cat #92256 and #92258) according to manufacturer’s instructions. After immunostaining, coverslips were washed and RL1 antibody was immunostained with goat-anti-mouse-Alexa-Fluor-488 secondary antibody (Jackson ImmunoResearch Laboratories, cat#115-545-003) diluted 1:500 in PBS/NGS) for 1 hr in a dark humid chamber at RT and washed. The DNA was counterstained with 200 ng/ml DAPI in PBS for 10 mins at RT and washed. Finally, the coverslips were mounted on glass slides using an antifade mounting media (10% w/v Mowiol 4-88(EMD Millipore)/1% w/v DABCO (Sigma-Aldrich)/25% glycerol/0.1 M Tris, pH 8.5) and cured overnight before proceeding to imaging.

Microscopy and image segmentation

Locales were imaged using the OMX-V4 microscope (GE Healthcare) equipped with a U Plan S-Apo 100×/1.40-NA oil-immersion objective (Olympus), two Evolve EMCCD cameras (Photometrics). Z-sections were 125 nm apart. Images were segmented in 3D using a combination of in-house softwares and the Allen Cell and Structure Segmenter (Segmenter) toolkit in Python (Chen et al., 2020). Briefly, nuclei were first identified by segmenting DAPI using Otsu thresholding method. RL1 signal was first segmented using a modified nup153 workflow from Segmenter. In H1 cells, there are some cytoplasmic signals for RL1 which were removed using the DAPI nucleus mask. RL1 signals are puncta and the alpha-shapes algorithm (Edelsbrunner et al., 1983) was used to calculate a 3D mesh-surface enclosing all the puncta (SFig 1A). This surface was used to measure the nuclei’s volume (V), surface area (SA), and SA/V. We segmented SON using the SON workflow, and MKI67IP using “masked object thresholding” workflow from Segmenter. The binary segmented images were used to measure the total number of locales in each nucleus and average and total V, SA. Nuclei images and segmentations were manually curated for quality control and are available online. To measure locale distances in each cell line, the SON and MKI67IP binary images were converted to mesh surfaces using a marching cube algorithm. An asymmetric locale to locale distance was then calculated by measuring the closest distance between a pair of vertices from locales mesh-surfaces (SFig. 1B). The PCA analysis was performed in R using 15 numerical metrics for each nucleus (Fig. 1F and SFig. 1F).

HCT116 EGFP-SON KI

A Cas9/sgRNA RNP complex and linear dsDNA donor were transfected into HCT116 cells using the Amaxa Nucleofector II device (Lonza) set for program D-032 with the Cell Line Nucleofector Kit V. We obtained the guide and CRISPR RNA using the GeneArt Precision gRNA Synthesis Kit (Thermo Fisher Scientific) with our designed SON guide sequence 5’-gagagaacggagcggacgcca -3’. The Cas9 protein Alt-R S.p. Cas9 Nuclease V3 was purchased from Integrated DNA Technologies (IDT). Donor DNA was amplified by PCR from a previously modified version of the SON-containing BAC RP11-165J2 (Invitrogen) with EGFP sequence at the NH2 terminus of SON (Khanna et al., 2014). We used modified primers containing a 5’ amine group with a C6 linker (Yu et al., 2020) 5’-GTGCTCACTGATTGGTCCCTC - 3’ and 5’ - GAACGACTGCGTCTCCGAAG - 3’ (amC6, IDT) to obtain a 1044bp fragment that includes the EGFP sequence flanked by a 160bp upstream and a 165bp downstream homologous region of SON. After 5 days of recovery, EGFP-positive cells were sorted using the FACS ARIA II sorter at the UIUC Roy J. Carver Biotechnology Center. Individual HCT116 clones expressing SON-GFP were genotyped by PCR with primers 5’-CAAGAGAGACGGCTCCTGTAATG -3’ and 5’-TCGGAAAAGGCGAAGTTCCTCG -3’ which amplify the complete EGFP integration to identify the double allele knock-in clone “C4”. HCT116 EGFP-SON clone C4 was transduced with a lentivirus carrying the Lenti-mCherry-PCNA construct (Xiong et al., 2023) by supplementing media with 8 μg/ml polybrene (Santa Cruz Biotechnology).

Live-cell Imaging and analysis

HCT116 EGFP-SON clone C4 cells expressing mCherry-PCNA were seeded in 35-mm dishes with a #1.5-thickness glass coverslip bottom (MatTek, P35G-1.5-14-C) and grown to ∼90% confluency. Live cell imaging was performed using a OMX V4 (Applied Precision) microscope with a 100×/1.4 NA oil immersion objective (Olympus), two Evolve EMCCDs (Photometrics), and a live-cell incubation chamber. The chamber contains a temperature control (37°C) for the incubator and the objective lens, and a humidified CO2 supply. 3D images were acquired using a z-spacing of 300nm once every 5 mins for 12 hrs. Each z-slice was imaged at 2.0% transmittance and 10ms exposure for 488nm excitation (GFP, Speckles), and 5.0% T with 10ms exposure for 568nm (mCherry, PCNA). We used FIJI to smooth images and remove noise. A median filter of 0.24μm and “rolling ball” background subtraction of ∼0.5μm was used. All measurements were done from corresponding individual Z-sections per timepoint. To identify PCNA foci we used the FIJI maximum entropy threshold plugin followed by watershed segmentation of the grayscale image. FIJI’s particle analysis was used to identify the center of each outlined PCNA foci, and the distance to nuclear speckle was measured from the PCNA foci center to the visible edge of the nearest nuclear speckle.

K562 LMNA and LBR single and double knock-out generation

To create LMNA and LBR knockout (KO) lines and the LMNA/LBR double knockout (DKO) line, we started with a parental K562 cell line expressing an inducible form of Cas9 (Brinkman et al., 2018). The single KO and DKO were generated according to the procedure described previously (Schep et al., 2021).

TSA-seq

SON and MKI67IP TSA-Seq was performed using Condition E (labeling with 1:300 tyramide biotin, 50% sucrose and 0.0015% hydrogen peroxide) (Zhang et al., 2020) with the following minor modification: 150ul of Dynabeads M-270 streptavidin (Invitrogen, catalog no. 65306) was used to purify the biotinylated DNA. For LMNB1 TSA-Seq, Condition A2 (Zhang et al., 2020) (labeling with 1:10000 tyramide biotin, 50% sucrose, 0.0015% hydrogen peroxide and reaction time 20 min at RT) was used.

DamID-seq and data processing

The DamID-seq and subsequent processing was performed according to a previously published protocol (Leemans et al., 2019).

RNA-seq (for K562 DKO)

RNA was isolated using the Qiagen RNeasy column purification kit, after which sequencing libraries were prepared using the Illumina TruSeq polyA stranded RNA kit.

Repli-seq and data processing (for K562 DKO)

The 2-fraction Repli-seq used in analysis of the K562 DKO cell line was done according to a previously published protocol (Marchal et al., 2018). Briefly, BrdU was added to a final concentration of 100 μM and cells were incubated at 37°C for another 2 h. Cells were fixed in ice-cold ethanol and processed for Repli-seq as described before.

Genomic segmentations

For detecting SON TSA-seq Type I and Type II peaks, the SON TSA-seq score (25 kb bins) first was smoothed using locally weighted regression (LOESS). Next, a local maximum filter with window size (w) was applied, which recorded the maximum SON TSA-seq score within this window for each genomic locus. Then the data after the maximum filter were subtracted from the original smoothed data. Peaks were retained only if they showed a reduction in value, which indicates true local peaks. A parameter optimization procedure was conducted in order to maximize the peak agreement between two SON TSA-seq replicates. As a result, a smoothing factor (span=0.005) and a window size (w=50) were selected. To distinguish local peaks into Type I and Type II, we overlapped them with Hi-C subcompartments (Rao et. al, 2014). Genomic loci overlapping A1 subcompartments were classified as Type I peaks, while those overlapping A2/B1 were designated as Type II peaks.

ROIs in H1 cells corresponding to off-diagonal 100 kb genomic bins in SON versus LMNB1 TSA-seq scatterplots were defined as those bins below the line LMNB1 = -1.5ξSON – 1.3, in these scatterplots. H1 ROI genomic regions were further subdivided into H1 ROI-C1 and C2 sub-regions (Supp. Fig. 6A), in which the H1 ROI-C2 regions were defined as the genomic bins below the LMNB1= -1.5ξSON – 2.2 line, while ROI-C1 regions were all remaining ROI bins.

To define LADs, p-w-v-fiLADs, and v-fiLADs, first we identified a consensus set of all chromosome regions which are LADs in at least one of the four examined cell types. Segmentation into these three domain classes was implemented heuristically using the local relative increase (”enrichment”) of the domain DamID score relative to its flanking regions (SFig. 4B): LAD domains, defined by their positive DamID scores using a Hidden Markov model, showed the largest positive local enrichment. We defined v-fiLADs as those domains which changed from LADs in a different cell line to iLADs and which showed lower local enrichment than 95% of the LADs. P-w-v-fiLADs were then defined as all other fiLAD domains and showed intermediate local enrichment levels.

To define LAD clusters in HCT116 cells based on LMNB1 vesus SON TSA-seq scatterplots, all 100kb genomic bins corresponding to LADs based on the HCT116 DamID (4DN file) were divided into 4 quadrants (Fig. 6E) centered at SON = -3.0 and LMNB1= 0.25. The C1, C2, C3, and C4 clusters were defined as the top-left, top-right, bottom-left, and bottom right quadrants, respectively.

Equatorial (Cluster 1) and non-equatorial (Cluster 2) LADs were defined in HFF cells (Fig. 7B) as those 100 kb bins which were LADs as defined by the LMNB1 DamID which had LMNB1 TSA-seq (4DN file number) scores higher than 0.75 for the equatorial LADs or SON TSA-seq (4DN file number) scores less than 0.98 for the non-equatorial LADs.

Gene expression analysis

Regarding Fig. 6H, the genes overlapping with C1-C4 LAD clusters were identified and their expression was evaluated in HCT116 RNA-seq.

Correlations with other genomic features

Data sets used for this study (Supplementary Table 1) include TSA-seq, DamID-seq, Repli-seq, RNA-seq, ChIP-seq and HiC compartment scores. TSA-seq, DamID, Repli-seq, and ChIP-seq data were mapped to hg38 and binned in 25 kb genomic windows, except for analysis included in Fig. 6 and SFig. 6.

For the ChIP-seq analysis in Fig. 6, fold-change over input was calculated for 25kb bins and was converted to percentile values across the genome. The percentile values were averaged over all the bins corresponding to the C1-C4 LAD clusters in the Fig. 6F heatmap or were overlayed on the LMNB1 versus SON TSA-seq binned scatterplots in Fig. 6G.

A cLAD score was calculated using the DamID LAD calls in 7 available human cell lines (STable 1). First the union of all LAD calls among the 7 cell lines was computed. For each 100kb genomic bin a number between 1 to 7 was assigned depending on the number of cell lines in which the bin was called as a LAD. The average cLAD score was overlayed on the LMNB1 versus SON TSA-seq binned scatterplot in Fig. 6J.

Correlating omics data with the IMR90 multiplex FISH data

For the LMNB1 TSA simulation (Fig. 7A), a 3D image of the segmented HFF nuclear periphery (anti-RL1 staining) was used. The binary image was convolved with a symmetrical 3D kernel based on the exponential decay function, Be^-R*d, of the TSA-biotin signal in which R is the decay constant and d is the distance to the central element of the kernel (B= 5.86, R= 3).

To correlate HFF LAD clusters with the FISH data from Su et al. (Su et al., 2020) (Fig. 7C-D) the x-y coordinates of FISH probes were first rotated using Principal Component Analysis (PCA) such that the length (PC1) and width (PC2) of each nucleus aligns with its PC1 (X) and PC2 (Y) axes. The size of each nucleus was then normalized to fit in a 4×2×1 (X*Y*Z, arbitrary units) cuboid using a min-max normalization. Probes mapping within equatorial (EQ, Cluster 1) and non-equatorial (non-EQ, Cluster 2) LAD clusters were selected and overlayed to generate Fig. 7C. Similarly, to correlate high-throughput FISH data with distance to centromeres, probes mapping to chromosome 4 (Fig. 7D) were selected and their distances to centromeres were calculated using UCSC chromosome banding data (STable 1).

Data access and code availability

Bed files corresponding to the genomic domains and datasets that we are uploading or previously have uploaded to 4DN can be found on the manuscript webpage on the 4DN website (https://data.4dnucleome.org/belmont_lab_nuclear_locale).

The K562 LMNA, LBR, and LMNA/LBR DamID-seq, Repli-seq, and RNA-seq raw and processed data has been deposited to GEO (GSE263012). The locale imaging data and the corresponding segmentations can be found on Illinois Data Bank (https://doi.org/10.13012/B2IDB-9792611_V1). The code used to generate figures in this manuscript, along with processed input and output files can be found have been reposited to the Illinois Data Bank (https://doi.org/10.13012/B2IDB-4383352_V1).

Acknowledgements

This work was supported by the National Institutes of Health Common Fund 4D Nucleome Program grants U54DK107965 (A.S.B., B.v.S., D.M.G., and J.M.) and UM1HG011593 (J.M., A.S.B., and D.M.G.). Research at the Netherlands Cancer Institute is supported by an institutional grant of the Dutch Cancer Society and of the Dutch Ministry of Health, Welfare and Sports. The Oncode Institute is partially funded by the Dutch Cancer Society.

Author contributions

OG, TvS, YW, PK, LZ, YZ are co-first authors based on their contributions described below. This study was a highly collaborative group effort between the laboratories of ASB, BvS, DMG, and JM as a research center during the of NIH’s 4D-Nucleome’s consortium phase I. The analysis performed here were led by two PhD students, OG and TvS from the laboratories of ASB and BsV, respectively. They organized weekly online meetings among all the trainees with monthly guidance from the corresponding authors.

Data gathering: OG (HCT116 and K562) and PK (H1, HFF and K562-LMNA/LBR DKO) performed the immunofluorescence staining of locales and collected images. LZ (HCT116, H1 and HFF) and PK (K562, H1 and HFF) performed LMNB1 TSA-seq experiments. PK performed and processed MKI67IP TSA-seq experiments. PK performed MKI67IP, LMNB1 and SON TSA-seq in K562 WT and K562-LMNA/LBR DKO cell lines. TvS performed LMNB1 DamID experiments. AEV performed the early/late repli-seq in K562-LMNA/LBR DKO cell lines. YZ and YW curated and performed quality control on RNA- and ChIP-seq datasets from external sources. They also created the infrastructure needed for organization and hosting of data generated for this study and external data for trainees to use.

Hypothesis generation and analysis: All authors contributed to hypothesis generation with the following specifics. OG created software for the segmentation and the analysis of nuclear locales microscopy data (Fig. 1 and Supp. Fig. 1). TvS performed most of the analysis for Fig. 2 and supp. Fig. 2. LZ performed analysis for Fig. 2C and SFig. 2D-E. YW performed most of the analysis for Fig. 3 and SFig. 3. PAZ performed the initial analysis for Fig. 3F, which was then reproduced by YW. GAHG performed live-cell imaging and the analysis of related to SFig. 3F. Initial idea for p-w-v fiLADS (Fig. 4) was developed by PK and ASB which was then realized in Fig. 4 and supp. Fig. 4 by analysis done by TvS. Analysis in Fig. 5 and supp. Fig. 5 were done by TvS. Analysis in Fig. 6 and supp. Fig. 6 were done by OG. Analysis in Fig. 7 and supp. Fig. 7 were done by OG and TvS. Ideas behind Fig. 8 were developed by OG and ASB and was realized by OG. Writing: Original draft was written by OG, TvS, YW, and ASB. Later drafts and final manuscript were written by OG and ASB with inputs from others.

Supplementary Figure Legends

Supplementary Figure 1:

Supplementary Figure 2:

Supplementary Figure 3:

Supplementary Figure 4:

Supplementary Figure 5:

Supplementary Figure 6:

Supplementary Figure 7:

Note

This reviewed preprint has been updated to add an equal contribution statement.