Replication dynamics identifies the folding principles of the inactive X chromosome

Chromosome-wide late replication is an enigmatic hallmark of the inactive X chromosome (Xi). How it is established and what it represents remains obscure. By single-cell DNA replication sequencing, here we show that the entire Xi is reorganized to replicate rapidly and uniformly in late S-phase during X-chromosome inactivation (XCI), reflecting its relatively uniform structure revealed by 4C-seq. Despite this uniformity, only a subset of the Xi became earlier replicating in SmcHD1-mutant cells. In the mutant, these domains protruded out of the Xi core, contacted each other and became transcriptionally reactivated. 4C-seq suggested that they constituted the outermost layer of the Xi even before XCI and were rich in escape genes. We propose that this default positioning forms the basis for their inherent heterochromatin instability in cells lacking the Xi-binding protein SmcHD1 or exhibiting XCI escape. These observations underscore the importance of 3D genome organization for heterochromatin stability and gene regulation.

In mammals, one X chromosome is inactivated in females to equalize X-linked gene expression relative to males during early development 1,2 . Microscopic studies have revealed that this Xi is more compact than the active X (Xa) 3,4 , leading researchers to explore the relationship between three-dimensional (3D) genome organization and transcription 5 .
A genome-wide chromosome conformation capture technology, Hi-C 6 , has revealed that mammalian autosomes are composed of Mb-sized topologically associating domains (TADs) 7,8 . TADs can be in either active (A) or inactive (B) nuclear compartments 6 , which are further subdivided into several subcompartments 9,10 . It was reported that when the Xi forms during mouse embryonic stem cell (mESC) differentiation, neighboring small A or B compartment domains on the Xi initially fuse with each other to form larger S1 or S2 compartment domains, which further merge to form the compact Xi structure through the actions of SmcHD1 (ref. 11), a global Xi-binding protein known for its role in X-chromosome inactivation (XCI) maintenance [11][12][13][14][15] . Unlike the autosomes, Xi lacks TADs and compartments but instead has a unique bipartite 'megadomain' structure separated by a tandem macrosatellite repeat, Dxz4/DXZ4 (refs. 9,16-18). While many XCI regulators have been identified, how the Xi acquires this unique, compact 3D organization is still unclear 2,5 .
An evolutionarily conserved chromosome-wide early-to-late (EtoL) replication timing (RT) switch of the Xi has long been known [19][20][21][22] . While its biological significance remains obscure, a tight relationship between RT, A or B compartments [23][24][25] and subcompartments 10 suggests that Xi's RT dynamics might reflect its compartmentalization. However, Xi's RT regulation has been studied primarily by microscopy or bulk genome-wide replication assays 13,19,[25][26][27][28] . Here we used bulk and single-cell DNA replication sequencing [29][30][31] to address the Xi's dynamic RT changes during mESC differentiation. Haplotype-resolved 4C-seq revealed that the Xi's RT indeed reflected its compartment organization. Because SmcHD1 is involved in the Xi's RT and compartment regulation in humans 12 and mice 11,[13][14][15] , we also used SmcHD1-mutant cells. Our results are consistent with the idea that the default 3D architecture of the X chromosome forms the basis for regional differences in Xi heterochromatin stability. Article https://doi.org/10.1038/s41594-023-01052-1 the postimplantation epiblast 19 . The Xi's EtoL RT change probably also occurs around days 5-7 of JB4/EI7HZ2 mESC differentiation, although we did not pursue this further (as distinguishing gradual RT changes in all cells versus cells with abrupt RT changes that were gradually increasing was a challenge).

RT delays and advances generate a uniformly late RT Xi
To further dissect the Xi's replication kinetics, we performed single-cell Repli-seq (scRepli-seq) with cells throughout the S-phase to construct what we call the 'whole-S' RT profiles (Fig. 1e). In NSCs, the whole-S RT profile of the JF1-Xa looked similar to the autosomes, as expected ( Fig. 1e and Supplementary Fig. 2d). By contrast, the B6-Xi initiated replication in the second half of S-phase and lacked clearly distinguishable early or late RT patterns (Fig. 1e).
Unlike the Xa, the replication score (percentage of replication) progression of the Xi was uncoordinated with the autosomes (Fig. 1f,g and Supplementary Fig. 2e,f). The steep rise in the Xi's replication score of the fitted sigmoid curve (Fig. 1f, red line based on scRepli-seq average) suggested fast completion of the Xi replication with a T 10-90% (defined as the time required for a chromosome to go from 10 to 90% replication, assuming a 10 h S-phase) of 3.4 h. By contrast, the T 10-90% of the Xa and the autosomes were significantly longer, on the order of 8 h ( Fig. 1f and Supplementary Fig. 2e). The replication scores of the Xi were variable among cells after 60-70% S-phase (Fig. 1g). This suggests large variability in the timing of replication initiation among cells and, in turn, the timing of replication completion of the Xi. Thus, the Xi's T 10-90% of 3.4 h is probably an overestimate.
While widespread EtoL RT changes were expected, the Xi's latest-replicating regions advancing their RT toward mid-late S in NSCs was unexpected (Fig. 1h), as they were classified as LtoL domains based on Repli-seq ( Fig. 1d and Supplementary Fig. 2c). Thus, scRepli-seq revealed that the Xi's RT change is chromosome wide, with both advances and delays to achieve its uniform RT in the second half of S. Although the Xi's RT is uniform and synchronous for a given cell, there is variability among cells with regards to when the Xi initiates replication, which is out of phase with the rest of the genome as if the Xi is disengaged from the genomic RT program.

SmcHD1 is required for maintaining the uniformly late RT Xi
SmcHD1 is required for Xi's late RT in human TERT-RPE1 (hTERT-RPE1) cells, mouse embryonic fibroblasts (MEFs) and mouse embryos [12][13][14] . To see whether SmcHD1 also affects Xi's RT during mESC differentiation, we generated SmcHD1-mutant mESCs by CRISPR-Cas9 with a guide RNA targeting its ATPase domain. We confirmed 5-and 10-bp deletions, premature stop codons biallelically on exon 3 and SmcHD1 loss by western blotting (Supplementary Fig. 1g,h). Differentiation states of SmcHD1-mutant and wild-type (WT) cells were similar (Supplementary Fig. 1a-f), and so were the RT profiles of their Xs in mESCs and day 7 or 9 cells ( Fig. 2a and Supplementary Fig. 3a-d). Thus, SmcHD1 was largely dispensable for the initiation of Xi's EtoL RT change, as in mouse

An mESC-based system to study the RT dynamics of the Xi
To explore the RT dynamics of the Xi as it forms, we need to distinguish the two Xs by single nucleotide polymorphisms (SNPs). We used an F1-hybrid mESC line, JB4/EI7HZ2, from a cross between JF1 female and B6 male 32 (Fig. 1a). In JB4/EI7HZ2 mESCs, CAGzeo is on the JF1-X at the Hprt locus and IRESneo is inserted on the B6-X in the 3′ untranslated region of Eif2s3x, an XCI escapee. In G418 + /zeocin + medium, XO cells lacking either of the drug-resistance genes are eliminated and the XX cells show 100% skewed XCI with B6-X being the Xi, because differentiated cells with JF1-Xi cannot survive due to CAGzeo inactivation ( Fig. 1a and Supplementary Fig. 1a). Cells with B6-Xi can survive because the escapee locus allows IRESneo expression from the Xi. We used two neural differentiation protocols to compare mESCs ('before XCI'), day 7 or 9 early neurectoderm cells ('during XCI') 25 , and neural stem cells (NSCs) obtained after roughly 3 weeks of differentiation ('after XCI') 33 (Fig. 1a). Their differentiation states were confirmed by PCR with reverse transcription, immunofluorescence and RNA-sequencing (RNA-seq) ( Supplementary Fig. 1a-e). Day 7 or 9 neurectoderm cells and NSCs showed roughly 75 and 90% Xist RNA cloud formation, respectively, suggesting nearly uniform differentiation and XCI ( Supplementary Fig. 1f).

The Xi becomes late replicating during mESC differentiation
Our routine genome-wide RT assay (Repli-seq) is based on BrdU immunoprecipitation (BrdU-IP) from early and late S-phase cells fractionated by flow cytometry, followed by next-generation sequencing (NGS) (Fig. 1b). Then, relative enrichment of early-and late-replicating DNA is analyzed to generate a genome-wide RT map. Repli-seq data of mESCs, day 7 or 9 neurectoderm and NSCs were highly reproducible (Supplementary Fig. 2a). Because Repli-seq allows cell-type profiling 25,27,34 , we compared our Repli-seq data with various other datasets 30 and confirmed distinct and proper differentiation states (Fig. 1c).
Haplotype-resolved Repli-seq of mESCs, days 7 or 9 neurectoderm and NSCs revealed that the B6-X is converted to a chromosome-wide late-replicating Xi by day 7 and maintained thereafter, while the JF1-X remains active and maintains similar RT profiles ( Fig. 1d and Supplementary Fig. 2b). We classified the ESC-to-NSC RT regulation into four groups based on the mean RT values of each 400-kb bin ( Fig. 1d and Supplementary Fig. 2c): early-to-early (EtoE, mean RT of more than zero in both), EtoL (mean RT of more than zero and less than zero in mESCs and NSCs, respectively, with a greater than 0.5 difference), late-to-early (LtoE, mean RT of less than and more than zero in mESCs and NSCs, respectively, with more than 0.5 difference) and late-to-late (LtoL, mean RT of less than zero in both). In contrast to the JF1-X, almost all early RT bins on the B6-X in mESCs became late in NSCs (99 out of 102 EtoL bins of 39.6 Mb), while EtoE and LtoE domains were nearly nonexistent ( Fig. 1d and Supplementary Fig. 2c). Analysis of CBMS1 mESCs using the same differentiation protocol found an EtoL RT change of the Xi during days 5-7 when the cells acquired a late epiblast fate 25 , consistent with the emergence of a late-replicating Xi in  Supplementary Fig. 1i), which are subject to immunoprecipitation (IP) using an anti-BrdU antibody. The ratio of early and late S-phase DNA is calculated to generate a genome-wide RT profile. c, Comparison of genome-wide RT profiles (sliding windows of 200 kb at 80-kb intervals, excluding the X) during differentiation of JB4/EI7HZ2 mESCs, CBMS1 mESCs 30 , EpiSCs and MEFs by hierarchical clustering. d, BrdU-IP RT profiles of the JF1-X and the B6-X during JB4/EI7HZ2 mESC differentiation (average of 2-3 replicates, 400-kb bins). RT classes represent four types of ESC-to-NSC RT regulation, which were described on the right. The blue line shows Dxz4 position. e, Binarized whole-S scRepli-seq profiles of the JF1-Xa and the B6-Xi from 94 NSCs throughout the S. The scRepli-seq profiles are ordered by their percentage replication scores (of the whole genome excluding the X). BrdU-IP RT of the NSCs is shown for comparison. The blue lines show Dxz4 position. f, Percentage replication scores of the JF1-Xa and the B6-Xi were plotted against those of the whole genome for each NSC (open circles). Cells were divided into 18 groups by their percentage replication scores and the average scores of the JF1-Xa and the B6-Xi in each group are shown by filled circles. To minimize errors during data fitting, we added two pseudo-percentage replication scores, 0 and 100%, at the earliest end and the latest end of the S-phase, respectively. We fitted a linear and a sigmoid model to the JF1-Xa and the B6-Xi, respectively, to calculate the T 10-90% values (the time required to go from 10 to 90% replication, assuming a 10 h S). g, Cells were divided into ten groups by percentage replication scores and the boxplots show distributions of percentage replication scores of the JF1-Xa and the B6-Xi in each group. h, An exemplary latest RT region on the Xa, which advanced its RT on the Xi in NSCs. CBMS1 mESC data are shown for comparison 30 .     Single-cell RT represents an estimated time when a given genomic bin replicates in S (assuming a 10 h S). P values were obtained from one-sided paired Wilcoxon signed-rank test. N, the number of bins in each RT class. d, Whole-S scRepli-seq plots in Fig. 2b were sorted by the BrdU-IP RT values of SmcHD1-mutant NSCs (from early to late RT). The red box corresponds mostly to the SD and CE classes and is magnified on the right. e, Comparison of percentage replication scores of the Xs in WT (black) and SmcHD1-mutant NSCs (red) with those of the whole genome, as in Fig. 1f.
Article https://doi.org/10.1038/s41594-023-01052-1 Using whole-S scRepli-seq, cell-to-cell RT heterogeneity of the SD domains was assessed. The SD domains replicated earlier in SmcHD1-mutant cells (Fig. 2b,c). By sorting the scRepli-seq RT bins according to the average BrdU-IP RT values, we found that most of the mutant NSCs completed the B6-Xi SD-domain replication before mid-S, while WT NSCs initiated SD-domain replication after mid-S (Fig. 2d). Thus, most cells exhibited SD-domain RT reversal without SmcHD1.
Because the SD domains occupied only about 10% of the Xi, the overall Xi replication duration judged by T 10-90% was similar between the mutant and WT NSCs (Fig. 2e). The mutant Xi initiated replication slightly earlier than the WT Xi (Fig. 2c,e). However, the non-SD portion of the mutant Xi still replicated rapidly and uniformly as in WT cells (Fig. 2b,d). Unlike the Xi, the JF1-Xa and the autosomes exhibited similar whole-S scRepli-seq profiles ( Fig. 1e and Supplementary Fig. 2d) and coordinated replication progression (Fig. 1f,g and Supplementary Fig. 2e,f), which were maintained in SmcHD1-mutant NSCs as assayed by BrdU-IP RT ( Supplementary Fig. 3a-d), whole-S scRepli-seq ( Supplementary  Fig. 3e,f) and single-cell RT (scRT) values ( Supplementary Fig. 3g).

SD domains protrude out of the Xi core in SmcHD1-mutant NSCs
While RT closely reflects the A or B compartments 23,25 and subcompartments 10 , Xi's compartments are elusive 38 . To test whether RT reflected the Xi's compartment organization, we used haplotype-resolved 4C-seq (refs. 39,40) ( Supplementary Fig. 4a,b) and analyzed nine viewpoints ( Fig. 3a and Supplementary Fig. 4c). The 4C-seq profiles of both Xs in WT and SmcHD1-mutant mESCs looked identical ( Fig. 3b and Supplementary Figs. 4d, 5 and 6). Early-replicating (SD, SI and CE) viewpoints interacted with other early-replicating domains, skipping the late-replicating domains in between, while late-replicating CL viewpoints showed interactions primarily within their resident late-replicating domains ( Fig. 3b and Supplementary Figs. 5 and 6), consistent with spatial segregation of early-and late-replicating compartments. On day 9, the B6-Xi in WT and SmcHD1-mutant cells still looked similar but showed less contrast between the peaks and valleys, resulting in less significant far-cis interactions than in mESCs ( Fig. 3b and Supplementary Figs. 4d and 6; far-cis contact numbers in the top right corner of each plot). In WT NSCs, the contrast was even weaker, consistent with Xi's reorganized RT ( Fig. 3b and Supplementary Fig. 6). In SmcHD1-mutant NSCs, the 4C-seq profiles of the SD viewpoints on the B6-Xi were markedly different and had peaks at SD-domain positions, consistent with their RT reversal ( Fig. 3b and Supplementary Figs. 4d and 6). In addition, the SD viewpoints in mutant NSCs showed weaker cis interactions with their surrounding regions compared to SI, CL and CE viewpoints ( Supplementary Fig. 7). These observations are consistent with the idea that SD domains protrude out of the Xi core and contact each other in SmcHD1-mutant NSCs.
We find much fewer defects for non-SD viewpoints on the mutant NSC B6-Xi ( Fig. 3b and Supplementary Figs. 4d and 6; except for the 45-Mb viewpoint (VP45Mb); discussed later). The JF1-Xa in WT and mutant cells maintained similar 4C-seq profiles during differentiation (Supplementary Figs. 4d and 5). These results indicate that the B6-X is gradually transformed into a uniformly compacted structure that lacks segregation into early-and late-replicating compartments during differentiation, confirming earlier studies 11,41 . In SmcHD1-mutant NSCs, however, the SD domains protrude out and contact each other.

SD domains are prone to reactivation in SmcHD1-mutant cells
SmcHD1-mutant cells exhibit Xi reactivation 11,13,14,35,42 . To examine its relationship to Xi compartmentalization, we reanalyzed RNA-seq data of WT and SmcHD1-mutant neural progenitor cells (NPCs) 11 , as their cell identity was close to NSCs ( Supplementary Fig. 8a). We compared the gene expression of different RT classes because RT accurately reflected the Xi compartmentalization ( Fig. 3b and Supplementary Fig. 6). The SD-domain genes showed the highest reactivation in SmcHD1-mutant cells (Fig. 4a,b), which largely overlapped with the SmcHD1-sensitive class I genes 11 ( Fig. 4a and Supplementary Fig. 8b). The SI domain showed significant reactivation but to a much lesser extent, while the CL and CE domains did not (Fig. 4a,b). Thus, Xi reactivation is strongly correlated with compartmentalization defects.
We further analyzed the relationship of Xi reactivation with epigenetic status. The SD domains coincided with Xi regions in SmcHD1-mutant NPCs showing a significant decrease in Xist RNA binding (Fig. 4a,c), a decrease in histone H3 lysine 27 trimethylation (H3K27me3) (Fig. 4a,d) and an increase in H3K4me3 (Fig. 4a,e). The SI genes did not show less Xist binding but showed a significant increase in H3K4me3, although to a lesser extent than the SD genes (Fig. 4a,c,e). The SD domains also coincided with regions depleted of H3K27me3 in SmcHD1-mutant MEFs 13 ( Supplementary Fig. 8c,d), although parts of the SD domains retained H3K27me3. Several SI domains close to the telomere also showed H3K27me3 depletion ( Supplementary Fig. 8d).
We analyzed SmcHD1 enrichment on the Xi using published data 11,13 . SmcHD1 was particularly enriched on the SD domains in MEFs but not NPCs ( Supplementary Fig. 8e), suggesting that SmcHD1 binding is cell-type specific. Nonetheless, SD-domain defects in SmcHD1-mutant MEFs and NPCs were similar based on H3K27me3 data (Fig. 4a,d and Supplementary Fig. 8c,d). Thus, the SD domain's susceptibility to SmcHD1 mutation is not due to the preferential binding of SmcHD1 to these domains.

DNA-FISH showed SD-domain protrusion from the mutant Xi core
To validate the Xist RNA binding data, we performed sequential Xist RNA fluorescence in situ hybridization (RNA-FISH) and DNA-FISH using two probe sets targeting neighboring SD, SI and/or CL domains (Fig. 5a). We categorized the DNA-FISH signal localization relative to the Xist cloud into four groups (Fig. 5b). In WT NSCs, most of the SD, SI and CL signals were positioned similarly close to the Xi surface (Fig. 5b,c,e and Supplementary Fig. 8f) or from the Xist cloud centroid (Fig. 5d,f), confirming earlier studies 43,44 . In SmcHD1-mutant NSCs, however, the SD probes frequently protruded out of the Xist cloud (Fig. 5c,e and Supplementary Fig. 8f) and became distant (Fig. 5d,f). The SI (but not CL) probes exhibited the same trend but to a lesser extent ( Fig. 5c-f).
To validate the protrusion and interaction of SD domains, we performed pair-wise DNA-FISH using three SD probes (Fig. 5g). Simultaneous protrusion of two probes was much more frequent in SmcHD1 mutant than WT NSCs (Fig. 5h,i). Focusing on the 'two-SD protrusion' cells, the 71-SD to 98-SD distance was significantly closer in SmcHD1-mutant NSCs (Fig. 5j,k), consistent with 4C-seq ( Fig. 5g, pair-1). The 98-SD to 148-SD distance was similar in SmcHD1-mutant and WT NSCs (Fig. 5j,k), again consistent with 4C-seq (Fig. 5g, pair-2 interacts in both WT and KO, although slightly higher in KO). The 71-SD to 148-SD distance did not decrease in SmcHD1-mutant NSCs (Fig. 5j,k), consistent with their weak interaction by 4C-seq (Fig. 5g, pair-3). Thus, while there is some cell-to-cell heterogeneity, FISH results were consistent with 4C-seq (see Supplementary Text 1 and 2 for further discussion).

SD but not SI domains on both Xs contact other chromosomes
Is the SD domain positioning on the Xi surface a cause or a consequence of SD gene reactivation in SmcHD1-mutant cells? The answer was brought about serendipitously when we analyzed interchromosomal interactions by 4C-seq (refs. 39,45). In SmcHD1-mutant NSCs, the B6-Xi SD viewpoints frequently contacted other chromosomes, while the non-SD viewpoints rarely did so ( Fig. 6a and Supplementary Fig. 9a), consistent with the idea that the SD but not SI and CL domains are on the outermost surface of the Xi, allowing contact with other chromosomes. This was true for the B6-Xi in WT NSCs and, also for the JF1-Xa ( Fig. 6a and Supplementary Fig. 9a), suggesting that this SD domain positioning is independent of transcription or the SmcHD1 genotype.      FISH experiments suggested that the SD, SI and CL domains are equally close to the Xi surface (Fig. 5c-f), possibly due to the resolution limit. Based on the interchromosomal interaction frequency, we predict that the SI and CL domains are just underneath the SD-domain layer in NSCs, avoiding contact with other chromosomes. The B6-Xi SD domains in mutant NSCs make more interchromosomal contacts than in WT NSCs ( Fig. 6a and Supplementary Fig. 9a), possibly reflecting their protrusion out of the Xi (Fig. 5c-f,i).
A similar interchromosomal interaction trend was observed for the mESC (Fig. 6b and Supplementary Fig. 9a). To test whether this is true chromosome wide, we performed virtual 4C to analyze interchromosomal interaction using Hi-C data. We analyzed 'virtual' viewpoints throughout the Xs using mESC and NPC Hi-C data 11 , after confirming the remarkable similarity of actual and virtual 4C profiles (Supplementary Fig. 9b).
Virtual 4C supported the layered X-chromosome organization in NPCs, with the SD and CE domains showing the strongest interchromosomal interaction, followed by SI, then CL domains on both Xs (Fig. 6c,d; Xi's CE contains Xist, which is highly expressed; Supplementary Text 3). The SD domains in SmcHD1-mutant NPCs showed higher interchromosomal interaction frequency than in WT NPCs (Fig. 6c,e), consistent with their protrusion out of the Xi core in the mutant. Such an interchromosomal interaction trend was conserved in mESCs (Fig. 6c,d).
However, we found that the SI viewpoints, especially VP45Mb, and the CE viewpoint showed more interchromosomal interactions in mESCs than NSCs on both Xs (Fig. 6a-c and Supplementary Fig. 9a). As mentioned earlier, VP45Mb was an exceptional SI viewpoint that resembled the SD viewpoints' intrachromosomal interaction pattern, acquiring interactions with other SD domains in SmcHD1-mutant NSCs ( Supplementary Fig. 6). Therefore, VP45Mb could be a marginal SI viewpoint located close to the SD layer, showing frequent long-range intra-but not interchromosomal interactions in SmcHD1-mutant NSCs and NPCs. Likewise, the SD and CL layers could be further subdivided.
The viewpoint near VP100Mb maintained frequent interchromosomal interactions in mESCs and NPCs by virtual 4C (Fig. 6c) but not actual 4C (Fig. 6a,b), possibly due to the difference in the viewpoint size and resolution (virtual 4C, 0.4 Mb; actual 4C, 8.6 kb). In addition, some SD or SI virtual viewpoints showed interchromosomal interactions in NPCs but not mESCs ( Fig. 6c; for example, 130-160 Mb). Overall, however, virtual 4C data are consistent with a layered organization of both Xs regardless of cell types, XCI states or SmcHD1 genotypes. The SD domains showed earlier RT than the SI domains (Fig. 6f), which might reflect their distinct compartmentalization, extrapolating from the tight RT-subcompartment relationship 10 .

The outermost SD domains on the Xi are rich in XCI escapees
The XCI escapees are located near the Xi surface 16,39,43,44 . Given the outermost position of the SD layer, we asked how well the SD domains overlap with escapees. Using a comprehensive NPC escapee list 17 , the SD domains showed more than 54-and sixfold higher escapee density than the CL and SI domains, respectively ( Fig. 7a and Supplementary  Fig. 9c), a trend also observed in non-NPC cells 46,47 ( Supplementary  Fig. 9c,d). Consistently, escapees were significantly enriched in regions with frequent interchromosomal interactions (Fig. 7b). The escapees were positioned on the surface of both Xs based on interchromosomal interaction frequency (Fig. 6a-d), meaning that this is their default positioning regardless of expression states.
To test whether this is conserved in other species, we asked whether human Xi regions with frequent interchromosomal interactions are also escapee rich. We used a published hTERT-RPE1 Hi-C data 18 and performed virtual 4C. Although escapees are not necessarily conserved between mice and humans 1 , human escapees were indeed significantly enriched in Xi regions with frequent interchromosomal interactions ( Fig. 7c and Supplementary Fig. 9e), suggesting their positioning close to the Xi's outermost surface.
Thus, it is possible that escapees escape XCI because their default positioning on the Xi surface makes their repression state inherently unstable, predisposing them to be easily reactivated. Likewise, SD-domain gene reactivation in SmcHD1-mutant NPCs (Fig. 4b) can be interpreted as preferential derepression of Xi surface domain genes by the absence of SmcHD1. Consistent with this idea, even the nonescapees are reactivated more strongly in SD domains compared to SI and CL domains in SmcHD1-mutant NPCs (Fig. 7d). However, escapees were more strongly reactivated in all domain categories analyzed (Fig. 7d), suggesting that both domain-level and individual gene-level factors contribute to XCI escape.
Last, we asked whether strong gene promoters could drive the characteristic SD-domain positioning. As strong promoters exhibit high CpG density 48,49 , we used this as a proxy for strong promoters and analyzed their distribution on the Xi. We observed a similar distribution of high CpG-containing promoters in the SD and SI domains and their closer proximity in SmcHD1-mutant NSCs by FISH. a, SD, SI and CL domain BACs were used to examine their subnuclear localization relative to the Xist cloud. b, Representative images of RNA-and DNA-FISH using Xist RNA and SD, SI or CL BAC probes in JB4/EI7HZ2 NSCs. We categorized the DNA-FISH signal localization relative to the Xist cloud into four groups (Methods). Scale bars, 5 μm. c,e, Percentages of DNA-FISH signal localization relative to the Xist cloud in NSCs: set 1 (c) and set 2 (e). d,f, Distance between the DNA-FISH signal and the Xist cloud centroid in NSCs normalized with Xist Feret diameter in the same nucleus: set 1 (d) and set 2 (f). g, Three SD BACs and their interaction frequencies as observed by 4C-seq ( Fig. 3b and Supplementary Fig. 6). Red lines, viewpoints; blue lines, Dxz4; gray bars beneath each plot, significant far-cis contacts; colored bars, RT classes. h, Representative images of RNA-and DNA-FISH using Xist RNA and SD BAC probes in JB4/EI7HZ2 NSCs. The SD-SD localization patterns were categorized into three groups. Scale bars, 5 μm. i, Percentages of SD-SD localization patterns in WT and SmcHD1-mutant NSCs. j,k, Normalized distances between two SD probes were plotted as cumulative plots: comparing WT versus SmcHD1-mutant NSCs (j) and comparing the SD protrusion group (k). N, the total number of cells analyzed from two to three independent experiments. P values in c, e and i were obtained from a chi-square test for all groups and a Fisher's exact test for the protruded group. P values in d and f were obtained from one-sided Wilcoxon signed-rank test, with data that used a paired test indicated. P values in j and k were obtained from a two-sided Kolmogorov-Smirnov test.  11 . c, Similar plots to b were made after virtual 4C (388 viewpoints, 400-kb each) of the human X chromosome using a hTERT-RPE1 Hi-C data 18 . The inactive (genes subject to stable XCI), variable (genes variably escaping from XCI) and escapee genes were defined based on a systematic human transcriptome study 56 . Shown are their densities normalized by UCSC gene density. P values in b and c were obtained from Kendall's rank correlation test with the Bonferroni correction. d, Comparison of Xi probability of X-linked gene expression in SmcHD1-mutant (KO) and WT NPCs as in Fig. 4a,b. Genes were classified into nonescapees and escapees 17 and into SD, SI, CL and CE classes. P values were obtained from one-sided paired Wilcoxon signed-rank test. N, the number of X-linked genes in each class. e, Hi-C heatmaps of the mus-Xi in WT and SmcHD1-mutant (KO) NPCs from Wang et al. 11 (250-kb resolution). Black circles on Hi-C heatmaps indicate strong SD-SD interactions found in SmcHD1-mutant NPCs. The bottom panel shows enlarged views of blue boxes in the top panel. Magenta, SD domains; blue arrows, Dxz4. f, Aggregated plots of interactions between SD (left) and random bins (right) by Hi-C. g, A proposed model for the formation of the multi-layered 3D structural organization of the Xi. During XCI, an EtoL compartment switch of the Xi occurs first (SD and SI domains) and is followed by a further 3D reorganization of the Xi through the actions of factors such as SmcHD1. Without SmcHD1, the Xi fails to maintain its late replication and 3D structure later during differentiation, resulting in frequent protrusion of SD domains located close to the surface of the Xi (but not SI domains), which occasionally interact with each other. Whereas the figure shows the contact between protruded SD domains, it is also possible that SD-SD domain interactions could occur inside the Xi core (Supplementary Text 2). Article https://doi.org/10.1038/s41594-023-01052-1 ( Supplementary Fig. 9f), confirming a previous report 50 . Thus, strong promoters cannot account for the SD-domain positioning.

SD-SD domain interactions are captured by Hi-C
One confounding enigma was that the previous Hi-C did not identify the SD-SD interactions in SmcHD1-mutant NPCs 11 . However, Hi-C does capture what 4C captures, as the virtual 4C data derived from these Hi-C data resembled our actual 4C results and showed SD-SD interactions ( Supplementary Fig. 9b). This led us to perform virtual 4C for viewpoints throughout the Xs in NPCs and plot their z scores (Supplementary Fig. 10a-f), which immediately visualized strong SD-SD interactions on the mutant but not WT NPC Xi (Supplementary Fig. 10b,d,f,  circled regions). Thus, we revisited the original Hi-C heatmaps 11 and successfully identified the SD-SD interactions in mutant but not WT NPC Xi (Fig. 7e, circled regions). Aggregated plots also revealed strong SD-SD interactions specifically in SmcHD1-mutant NPCs (Fig. 7f).
We also found that the red and blue patterns of the Xa z-score plots resembled the Hi-C principal component 1 (PC1) profile to some extent ( Supplementary Fig. 10a,c). This resemblance was also observed on the mutant NPC Xi (Supplementary Fig. 10d), perhaps reflecting the S1 and S2 compartments 11,51 . However, the red and blue patterns on the mutant and WT NPC Xi were similar ( Supplementary Fig. 10b,d), suggesting the presence of S1 and S2 compartments on the WT NPC Xi. When we calculated PC1-4 on the Xs and focused on those with more than 8% contribution rates, PC3 of WT NPC Xi resembled the S1 and S2 compartments, which was also true in mouse Patski cells 52 ( Supplementary Fig. 10g-i).
Thus, although the Xi is more compact than the Xa, they still seem to share certain structural features. While the megadomain structure becomes prominent as the Xi becomes compact ( Supplementary  Fig. 10g,j), the S1 and S2 or A and B compartment features still remain on NPC Xi ( Supplementary Fig. 10g,h). When the Xi heterochromatin is disturbed by perturbations such as SmcHD1 mutation, the relative contribution of different structural features change, leading to, for instance, the emergence of S1 and S2 compartments as PC1 (ref. 11) or enhanced TAD boundaries in SmcHD1-mutant cells 11,13 . In addition, new structural features, namely protrusion of SD domains and their interactions, emerge on the loss of SmcHD1.

Discussion
In this study, we used scRepli-seq and analyzed the Xi's RT regulation during mESC differentiation, and explored its compartment organization. Our results demonstrate that (1) the entire Xi is replicated rapidly and uniformly in late S in a given cell but with cell-to-cell RT heterogeneity ( Fig. 1e-h); (2) the Xi has a 'layered' organization, which is present already in mESCs, suggesting that this is the default architecture of the X independent of XCI (Fig. 6); (3) SmcHD1 is required for maintaining the XCI state, late-S replication and proper 3D architecture of domains located on the Xi's outermost layer (Figs. 2a-c, 3b and 4-7); (4) the Xi's outermost layer is escapee rich in mice and humans, which becomes preferentially reactivated in SmcHD1-mutant mouse cells (Fig. 7a-d) and (5) the 3D organization defects of the Xi surface domains in SmcHD1-mutant cells can be captured by Hi-C (Fig. 7e,f). Taken together, while the Xi appears to be uniformly compacted in 3D, our results indicated that the positioning relative to the Xi surface can explain regional differences in Xi heterochromatin stability, which was manifested on SmcHD1 mutation or XCI escape (Fig. 7g).
The Xi's chromosome-wide late-S replication has been known since 1960 (ref. 20). However, because roughly 70% of the Xa is late replicating, the EtoL RT switch was assumed to affect only the remaining roughly 30% of the future Xi (Fig. 1d). This was not the case, however, and both RT advances and delays generated a uniformly late-replicating Xi that completed replication within a few hours (Fig. 1e-h), which was consistent with and explain the rapid replication of the mouse C2C12 myoblast Xi in mid-to late S by live imaging 26 . Because developmental RT changes reflect preceding A and B compartment changes 25 , we predict that the Xi's RT changes also reflect its compartment reorganization. Our results validate and update the 'fast and random' Xi replication model 28 and reveal the Xi's chromosome-wide RT reorganization process.
We found relatively large cell-to-cell Xi RT heterogeneity within the second half of S ( Fig. 1e-g). Because different Hi-C subcompartments show distinct RT 10 , the cell-to-cell RT heterogeneity of the Xi and its RT being out of phase with the rest of the genome likely reflect Xi's compartment states. For instance, the Xi may form its own compartment, which could be slightly variable among cells regarding its position relative to the late-replicating B compartment, resulting in cell-to-cell RT heterogeneity.
The Xi constantly changed its conformation during differentiation, even after becoming late replicating ( Fig. 3b and Supplementary Fig. 6). Thus, the Xi may undergo a two-step structural change in which the EtoL domains first alter their compartments in an SmcHD1-independent manner, followed by a process that requires SmcHD1 for further structural reorganization and/or maintenance (Fig. 7g). While the SmcHD1-mutant Xi exhibits pleiotropic defects in gene expression, DNA methylation and histone modifications 14,[35][36][37] , recent reports indicate a direct role of SmcHD1 in regulating the higher-order chromosome architecture 11,12 . In particular, SmcHD1 knockout (KO) in cells that have established the Xi clearly showed that the higher-order Xi architecture is altered without gross changes in gene expression 13,15 .
SmcHD1 shows binding affinity throughout the Xi 11,13,53 and yet the Xi's RT reversal in SmcHD1-mutant NSCs was specific to the SD domains. Our data indicate that being near the Xi surface makes genes in such regions relatively unstable for maintenance of silencing, possibly because it is easier to make them physically loop out of the Xi chromosome territory than genes underneath the surface layer. Looping out may be a cause or a consequence of gene reactivation. However, being on the X-chromosome surface before XCI cannot be a consequence of reactivation. Therefore, we favor the view that such default 3D architecture of the X forms the basis for regional differences in heterochromatin stability. It has been a longstanding debate whether a given 3D genome architecture is a cause or a consequence of transcription 54 . We believe that our observations represent one of the best examples corroborating the causal role of the 3D genome organization on genome function.
If the default architecture causally affects local differences in Xi heterochromatin stability, it should predefine the sensitivity to perturbations other than SmcHD1 KO. We found that XCI escapees were overrepresented in Xi's outermost surface in mice and human (Fig. 7a-c). By contrast, CL domains rarely contained escapees ( Fig. 7a and Supplementary Fig. 9c,d). These observations comprehensively validated the widespread notion that escapees are located near the Xi surface 16,43,44,55 .
Escapees were more strongly reactivated than nonescapees in all domain categories analyzed (Fig. 7d), indicating that some intrinsic factors render the escapees more easily reactivated 1,2,5 . However, SD-domain escapees were more strongly reactivated in SmcHD1-mutant NPCs than SI-domain escapees (Fig. 7d). In addition, some SD-domain nonescapees were strongly reactivated in mutant cells, which was much less prominent in SI domains (Fig. 7d). Thus, as the SD domains are closer to the Xi surface than the SI domains, this might make them inherently more unstable. In such ways, SmcHD1-mutant cells can serve as an excellent model to explore XCI escape. Moreover, our work demonstrates the potential functional significance of subcompartments, which is tightly coupled with distinct RT regulation 10 .
Interchromosomal interaction analyses by actual and virtual 4C played a critical role in our work. Moreover, we could identify SD-SD interactions in the SmcHD1-mutant NPC Hi-C data (Fig. 7e,f). Thus, routine 4C and Hi-C analyses can overlook important features, and 4C and Hi-C data contain more information to be discovered in the future.
Article https://doi.org/10.1038/s41594-023-01052-1 Given the discovery of the potential importance of the default structure of the X, the next challenge would be to decipher the molecular basis of such default 3D genome organization. However, while the interchromosomal interaction patterns were similar between the two Xs in WT mESC and NPCs or SmcHD1-mutant NPCs, there were small differences between mESCs and NPCs (Fig. 6a-c), indicating that the default structure can change during differentiation. In addition, SmcHD1-mutant cells exhibit slightly different RT in different cell types 13 , suggesting that SmcHD1's contribution to the Xi architecture is cell-type specific. This could be partly explained by cell-type specific Xi compartment organization. It will be important to distinguish the intrinsic ('default') and extrinsic ('epigenetic') aspects of compartment regulation and understand how they are coordinated.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41594-023-01052-1.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/. Article https://doi.org/10.1038/s41594-023-01052-1 with 2× SSC and dehydrated by sequential 5-min washes with 70, 90 and 100% ethanol at room temperature before being air-dried at 58 °C for 1 h. Slides were then denatured in 70% formamide in 2× SSC at 80 °C for 3 min, dehydrated by sequential washes with cold 70, 90 and 100% ethanol, and were again air-dried at room temperature until hybridization. After an overnight hybridization at 37 °C, slides were washed with 50% formamide in 2× SSC at 45 °C three times, washed again with 0.1% SSC at 60 °C three times and counterstained with 1 μg ml −1 DAPI in 2× SSC before mounting with Vectashield. We recorded DNA-FISH signals by DeltaVision at the same xy coordinates from Xist RNA-FISH. Images were analyzed using the Fiji software. Because xy coordinates had slightly shifted during the experiments, we corrected the positions of Xist RNA-and DNA-FISH images by TurboReg macro 59 using DAPI signals of each image as references.

Imaging analysis of Xist RNA-and DNA-FISH signals
The contour of the Xist RNA signals (Xist cloud) was automatically drawn in each nucleus using a custom Fiji macro. Briefly, an image channel corresponding to Xist RNA was subject to 'Enhance Contrast' by saturated at 0.3 with normalization, followed by smoothing through the 'Mean Shift' plugin with spatial value of 10 and color at 25. The image was then converted to binary, and the Xist cloud region of interest (ROI) was identified using the 'Analyze Particles' tool, which was defined as the outer edge. The inner edge of the Xist cloud was subsequently defined by drawing a scaled ×0.65 ROI. The xy positions of DNA-FISH signals were detected using the 'Find Maxima' tool with a prominence threshold of more than 120. Only nuclei with a single Xist cloud and two clusters of DNA-FISH signals were analyzed. The DNA-FISH signal localization relative to the Xist RNA cloud was manually scored based on four categories: fully overlapped (signal inside the inner edge of the Xist cloud), inner edge (signal inside the Xist cloud but in between the inner and outer edge), outer edge (signal just outside the Xist cloud at its outer edge) and protruded (signal outside the Xist cloud). The distance between the DNA-FISH signal relative to Xist centroid and between two DNA-FISH signals were calculated based on their xy coordinates. The distance was normalized by Xist Feret diameter in the same nucleus.

Sample preparation for RT profiling by BrdU-IP Repli-seq
We followed the BrdU-IP protocol as described 30,60 . Cells or EBs were incubated in a medium containing 10 mM BrdU for 2 h before cell collection. After trypsinization, single-cell suspension was fixed in 75% ethanol. For fluorescence-activated cell-sorting (FACS), we stained fixed cells with Propidium iodide (Nacalai, 29037) and used a Sony SH800 cell sorter (ultra-purity mode) to sort early-and late-S-phase cell population (at least 10,000 cells per fraction). We used a Bioruptor UCD-250 (Sonic Bio) for genomic DNA sonication (high-output mode), with ON/OFF pulse times of 30 s/30 s for 6 min in ice-cold water. BrdU-incorporated DNA was immunoprecipitated using anti-BrdU antibody (BD Biosciences Pharmingen, 555627). After BrdU-IP, immunoprecipitated DNA samples were subject to whole-genome amplification with a SeqPlex kit (Sigma, SEQXE). NGS libraries were constructed from early-and late-replicating DNA after whole-genome amplification with an NGS LTP Library Preparation Kit (KAPA, KK8232) according to the manufacturer's instructions and were subjected to NGS with the HiSeq X Ten system.

Sample preparation for scRepli-seq profiling
Single cells from the whole-S-phase were sorted with a Sony SH800 cell sorter using single-cell mode. Three gates corresponding to early-, mid-and late-S-phases were set before sorting to select the mode for binarization analysis (Computation associated with the RT profiling of single cells). Sample preparations were performed as described in ref. 31. In total, 120 and 96 single cells for WT and SmcHD1-mutant NSCs (112 S-phase cells and eight G1-phase cells for WT NSCs, and 88 S-phase cells and eight G1-phase cells for SmcHD1-mutant NSCs) were prepared and analyzed. After filtering out cells with X chromosome abnormalities and abnormal median absolute deviation scores, 95 S-phase cells and four G1-phase cells were subjected to scRepli-seq for WT NSCs, while 85 S-phase cells and four G1-phase cells for SmcHD1-mutant NSCs were subjected to scRepli-seq.

Allele-specific and nonspecific NGS mapping for RT profiling
Paired-end reads were used as single-end reads. The raw fastq files were trimmed to remove adapter sequences by using the trim_galore v.0.6.6 (-quality 20-phred33-length 35) and cutadapt program v.1.15 before mapping. We performed two-step adapter trimming, first removing the Illumina adapter on the basis of the index of each NGS library and then removing the SEQXE adapter 31 . For mapping to non-haplotype-resolved mouse mm9 reference genome, bwa (v.0.7.17-r1188) was used (command, bwa mem). The Picard tool was used to remove duplicated reads and defined with a MAPQ ≥ 10 as uniquely mapped reads. For haplotype-resolved analysis, we constructed the B6/JF1-specific diploid genome as described in ref. 58, which was used as a reference for the JB4/EI7HZ2 cell line. For haplotype-resolved mapping, bwa (v.0.7.17-r1188) was used (command, bwa aln ≥ bwa samse). Our in-house hTERT-RPE1 phased genome based on SNP information was used as a reference for the hTERT-RPE1 cell line. To obtain the SNP information, haplotype-phasing analysis of the hTERT-RPE1 genome was performed by using a combination of 10X Genomics linked reads, Hi-C data 18 and Strand-seq data 61 .We defined MAPQ ≥ 16 as allele-specific uniquely mapped reads. For reads uniquely mapped to the B6/JF1 diploid genome, we used the liftOver tool (UCSC Genome Browser) to convert the genome coordinates to the mm9. Among the unique reads, we filtered out duplicated reads with the chromosome start position and strand information identical to an existing read. We also filtered out reads that overlapped with the mm9 and hg19 blacklists 31 . Reads per sample are shown in Supplementary Table 2.

Computation associated with RT profiling of cell populations
We followed our established pipeline for BrdU-IP population RT analysis 30 . For non-haplotype-resolved analysis, after mapping, removing duplicate reads and filtering mm9 blacklist, we counted the reads of early-and late-S-phase BrdU-IP samples in sliding windows of 200at 80-kb intervals and performed reads per million normalization. Then, the ratio of early-S-phase to total read counts ((early-S reads)/ (early-S reads + late-S reads)) was calculated for each bin, and were further converted to make their distribution to fit within ±1. This value was defined as the BrdU-IP RT score of each bin. We filtered out bins whose total read counts were within the bottom 5% of all bins. For haplotype-resolved RT profiling, we followed the exact same procedures, using nonoverlapping 400-kb bins. All BrdU-IP RT profiles were quantile normalized before downstream analysis. We used published BrdU-IP data of CBMS1 mESCs (GSM2904968 and GSM2904969) and CBMS1 day 7 (GSM2905017 and GSM2905018) from Takahashi et al. 30 for comparison. Hierarchical clustering was done using Ward's method with Euclidean distance in R.

Computation associated with the RT profiling of single cells
We followed our established pipeline for scRepli-seq RT analysis 31 . First, we analyzed non-haplotype-resolved scRepli-seq data to obtain the percentage replication scores of the whole genome for each single cell in 100-kb nonoverlapping bins. Here, X chromosomes were excluded from the analyses to avoid the bias due to the late replicating Xi. For binarization, different options were applied to each cell depending on their FACS sorting gates (2-HMM option for early-S FACS gate: most. frequent.state = '1-somy'; 2-HMM option for mid-S and late-S FACS gates: most.frequent.state = '2-somy'). To analyze the X chromosomes, we performed haplotype-resolved scRepli-seq as described 30 in 400-kb nonoverlapping bins. Percentage replication scores of scRepli-seq results are shown in Supplementary Table 3. To obtain the single-cell

Nature Structural & Molecular Biology
Article https://doi.org/10.1038/s41594-023-01052-1 RT value of a given genomic bin, we first calculated the percentage replication value of each genomic bin (that is, the percentage of cells that have replicated the bin) using non-haplotype-resolved binarized whole-S scRepli-seq data. Then, the genomic bins were subdivided into one-percentile groups based on their percentage replication values, from the earliest (the top one-percentile group) to the latest-replicating group of bins (the bottom one-percentile group). Each of these one-percentile groups has a range of percentage replication values with upper and lower limits (note that the range is variable between different groups) and was assigned an average percentage replication value, which was converted to an S-phase time, that is, the single-cell RT value (0-10 h; with a resolution of 0.1 h), assuming a 10 h S-phase. Therefore, the single-cell RT value of a given one-percentile group represents the average time during the 10 h S-phase when the genomic bins within this group replicate. To obtain the single-cell RT value of each genomic bin on the X chromosomes, we first calculated the percentage replication value of each genomic bin using haplotype-resolved binarized whole-S scRepli-seq data in a way similar to the method described above. Then, a given genomic bin was assigned an single-cell RT value of the one-percentile group with a percentage replication value range that contains that of the given genomic bin.

Allele-specific 4C-seq experiments
The primary primer sequences for the nine viewpoints on the mouse X chromosome are shown in Supplementary Table 4. These viewpoints were selected based on RT domains containing SNPs that allow allele-specific 4C-seq to be performed, as previously described 24 . The 4C-seq inverse primers were designed to have the first read of the paired-end read (P5) to cover a portion of the viewpoint region (near HindIII) and read into the target ROI. The second read of the pair (P7) was designed to cover a portion of the viewpoint region (near DpnII) containing a SNP to distinguish the homologous chromosomes ( Supplementary Fig. 4a,b). The 4C-seq inverse primers were designed for the two-step PCR strategy as described in ref. 40. The primary primers are complementary to the ends of a viewpoint facing outward. These primers contained a portion of the Illumina adapter sequence required for the second PCR step. The second primers are universal primers carrying Illumina indexed adapter sequences (Supplementary Table 4). Thus, they can hybridize to the adapter sequences introduced earlier by the primary primers and amplify the first-round PCR products. This resulted in the complete Illumina adapter sequence for pair-end sequencing. 4C-seq was performed essentially as described 40,45 with modifications below. Briefly, 5-10 × 10 6 cells were cross-linked with 1% formaldehyde for 10 min at room temperature. Cross-linked cells were lysed and the nuclei were digested with a sino.-base cutter, 400 U HindIII (NEB, R0104T), overnight at 37 °C. Fragmented DNA ends were ligated by 50 U of T4 ligase (Roche/Sigma, no. 79900901, more than or equal to 5 U μl −1 ). Then, the purified DNA was cut again with a four-base cutter, 50 U DpnII (NEB, R0543S), overnight at 37 °C. DNA was ligated again by T4 ligase (Roche/Sigma, no. 79900901, more than or equal to 5 U μl −1 ) overnight at 16 °C to generate 4C-DNA. The 4C-DNA was purified by phenol and chloroform extraction and ethanol precipitation. 4C-seq libraries were first amplified from 800 ng of 4C-DNA per viewpoint with 16 cycles of inverse PCR using primary PCR primer sets. A typical 50 μl PCR reaction was performed with 200 ng 4C-DNA per reaction using Phusion High-Fidelity kit (F-553L). PCR products were purified using 0.8× Agencourt AMPure XP beads (Beckman Coulter, A63881), and eluted in 50 μl eluting buffer (QIAGEN, 19086). Then 5 μl of purified PCR products were used to amplify again with 20 cycles of PCR using secondary PCR primer sets. Final PCR products were cleaned up by the QIAGEN PCR purification kit. Diluted 4C-seq libraries were mixed with other libraries and subjected to paired-end sequencing using the HiSeq X Ten system. We read roughly 4 million reads per sample, but due to the nature of the 4C-seq library, which contained large size products, we usually obtained 1-2 million paired-end reads per sample (Supplementary Table 2).

Allele-specific 4C-seq data analysis
The 4C-seq reads were first separated specifically to B6 or JF1 based on SNPs in read 2 (from P7) using cutadapt (cutadapt -e 0-trim-n -g ^(Fw primer sequence) -G ^(Rv primer sequence + SNPs)-no-trimdiscard-untrimmed; Supplementary Table 5). Only reads with 0% mismatches to the expected sequence and SNPs were kept. Our 4C-seq libraries had an almost equal fraction of reads between two alleles (roughly 50% each; see Supplementary Tables 1 and 5), indicating an unbiased PCR amplification of libraries. P5 reads of each assigned allele were mapped to mouse mm9 genome as single-end alignment and analyzed using an R pipeline 40 with analysis mode: all,-wSize 201. To analyze significant far-cis interactions, we used established R pipelines 45 with parameters -w = 100, -W = 3000 and false discovery rate (FDR) of 0.01. To analyze significant interchromosomal interactions, we used established R pipelines 45 with parameters -w = 500 and FDR = 0.01. Hierarchical clustering was done using Ward's method with Euclidean distance in R.