Epigenetic-structural changes in X chromosomes promote Xic pairing during early differentiation of mouse embryonic stem cells

Determination of A-, B-, and M-Regions in Chromosomes from ES Cells and 2-Days Cells

A- and B-compartments for the first and second replicates of ES cells and 2-day cells were obtained as previously proposed, using Hi-C data (GSM3127755, GSM3127759, GSM3127756, GSM3127760) and GC-content of local chromatin regions of the mouse genome (mm9) [28].

Each genome region was defined as an A- or B-region if the region was labeled as an A- or B-compartment in the results from two biological replicates of each experiment. The genome region was defined as an M-region if the region was labeled as a different compartment in results from two replicates for each experiment. Regions with no Hi-C data were defined as NA regions and were expected to correspond to centromere or telomere regions.

Estimation of the Basic Conformations of Chromosomes Using Polymer from Hi-C Data and Radial Distributions of Each Region

A polymer exhibiting the average three-dimensional conformation of each chromosome was constructed using Hi-C data for ES cells and 2-day cells (https://doi.org/10.5281/zenodo.3371884) [28] based on a recently proposed consensus method [29–31] implemented in PASTIS as an MDS2 method [32]. Here, the MDS2 method is an extended multidimensional scaling (MDS)-based method that constructs the model of each chromosome as a strand of beads. In the first step of this method, the contact frequency between each pair of genome regions is converted to the pairwise distance between them by assuming a biophysical model of DNA, where the pairwise distance is assumed to be inversely proportional to the cube of the contact frequency. Next, a multidimensional scaling problem is solved to obtain the 3D structure that best matches the pairwise distance. This MDS-based method is known to well reproduce the global structure of chromosome territories and TAD-like structures.

For both ES and 2-day cells, merged Hi-C data obtained from two biological replicates were analyzed. The three-dimensional spatial position of each monomer in each polymer for each chromosome provided that of the corresponding genome region analyzed in the previous section (Supplementary Table S1). Based on our results, each monomer was attributed to the A-, B-, or M-region (Supplementary Figure S1a).

According to these position data, ND_A, ND_B, and ND_M were evaluated in ES and 2-day cells (Figure 2). The overall chromosome shapes and their epigenetic features were confirmed to be unchanged qualitatively, even if the three-dimensional structure of each chromosome was constructed using other methods implemented in PASTIS with the NMDS, PM1, and PM2 methods [32] (Supplementary Table S2–4). Therefore, the polymer inferred by the MDS2 method was defined as that of basic conformations of chromosomes for the following construction of coarse-grained blob models.

Determination of A-, B-, and M-Domains and Domain Groups in Each Chromosome

Each polymer model for each chromosome was divided into groups of successively connected monomers, called domains (Supplementary Figure S1a). Monomers acting as domain boundaries (boundary monomers) were obtained based on boundary scores of each genome region (monomer) that could be calculated from merged Hi-C data using a previously proposed method with the same parameter values (Supplementary Table S1) [28]. Here, the boundary score was obtained by smoothing the insulation score along the genome regions. The insulation score of each genome region was estimated as the interaction probability between the upper 200 kbp part and the lower 200 kbp part of each region calculated from the Hi-C contact matrix data, where the regions with lower scores exhibit stronger insulation (see [28] for more details). In this study, each domain was assumed to consist of one boundary monomer and the monomers between this boundary monomer and the upstream boundary monomers (Supplementary Figure S1a).

The epigenetic attribute of each domain was defined as the A-domain if more than 50% of the domain occupied the A-region, as the B-domain if more than 50% of the domain occupied the A-region, or as the M-domain if the domain was neither an A-domain nor a B-domain (Supplementary Figure S1a). Notably, the NA region was considered a B-region when determining the domain attributes because the NA region was expected to correspond to the centromere or telomere regions forming heterochromatin. Additionally, the domain group was defined as the group containing domains with the same epigenetic attributes successively connected along the polymer (Supplementary Figure S1a).

Description of Each Domain Group Using Blobs to Construct Coarse-Grained Models of Chromosomes

Each domain group was divided into blobs, forming a similar envelope shape as the monomer population in the blob (Supplementary Figure S1b), with elimination of small structural noise as follows (Supplementary Figure S1a).

First, k-means clustering of the three-dimensional coordinates of monomers in each domain group was performed. For this clustering, the number of clusters dividing each domain group in the n-th chromosome was assumed as the integer value with [number of domains in the domain group]/CND_n truncated. Here, CND_n was the characteristic number of domains in the n-th chromosome and was defined as the number of domains for which the physical distance between two loci increased successively on average with the increases in the domain number between the two loci. In other words, CNDn could be considered to have a value similar to that of the persistent length of n-th chromosome obtained when the chromosome was described as a chain consisting of domains as the units. By this definition, each domain group containing more domains than CND_n was divided into multiple clusters, enabling imitation of its curved shape by connection of these cluster centers. Here, each blob is defined as a particle that represents the chromatin region containing the domain of each cluster, and CNDn gives the scale of typical blob size as the number of domains. CND_n was evaluated by comparing the decay rate of the average PDND ([physical distance between two domains]/[number of domain connections between these domains]) as a function of the number of domain connections between the two domains. It was noted that the average PDND in each chromosome decreased monotonically (Supplementary Figure S2). Therefore, in this study, PDND in each chromosome was simply fit using a function ~ exp(-[number of domains between two loci]/L), and CND_n was assumed as L.

Next, each domain group was described by distinct blobs based on the results of clustering. If the monomer cluster consisted of only a set of monomers successively connected along the polymer, the set of monomers was described as one blob. If not, the monomer cluster was divided to subclusters in which monomers were successively connected along the polymer in each subcluster, and each set of monomers in each subcluster was described as one blob.

The coordinate and radius of each blob were defined as the centroid coordinate of blobs containing monomers and the standard deviation of the distance from the blob centroid to the contained monomers. The epigenetic attribute of each blob was given by that of the corresponding domain group and designated the A-, B-, and M-blobs. By connecting the neighboring blobs along the polymer, a coarse-grained model of chromosomes was constructed as blob chains involving structural, genomic, and epigenomic features (Supplementary Figure S1a, Supplementary Table S5). In this model, the distance between each pair of blobs in the same chain was defined as that in the basic structure of the blob chain model of chromosomes.

In recently proposed coarse-grained models of large-scale genome structures, each chromosome has been represented by a chain of particles with a constant radius, assuming that the length of the DNA sequence contained in each particle is constant [33–36]. Conversely, in the present model, the radius of each blob was given non-uniform because each blob was assumed to contain a DNA region of varying length depending on the characteristics of local chromatin structures. Using such blobs with heterogeneous radii, the shape of each chromosome could be more appropriately described with fewer particles.

Equation of Motion for Each Blob in the Coarse-Grained Model of Chromosomes

The dynamics of all intranuclear chromosomes (chromosomes 1–19 and the X chromosomes) of ES cells and 2-day cells were simulated using the motion of each blob (Supplementary Table S5), which was affected by the interaction potential and noise from the nucleoplasm. The motion of the i-th blob was assumed to obey the overdamped Langevin equation, as follows:

where x_i=(x_i,y_i,z_i) is the position of the i-th blob, V indicates the potential of the system, and γ_i and R_i(t) were the coefficient of drag force and the random force working on i-th blob from the nucleoplasm, respectively, where γ_i=6πηr_i with a nucleoplasm viscosity η and a radius of the i-th blob r_i. R_i(t) was given as Gaussian white noise and satisfied 〈R_i(t)〉=0, and 〈R_i(t)R_j(s)〉=2γ_ik_BTδ_ijδ(t–s), where k_B is the Boltzmann constant, T is the temperature, δ_ij indicates the Kronecker delta, and δ indicates the Dirac delta function.

The first term of the right-hand side of Eq. (1) indicates forces working on the i-th blob provided by the potential of the system, as follows:

where V_col is the potential of the force excluding the volume effects among blobs, V_chr is the interaction potential of the force sustaining the shape of each chromosome, and V_mem is the potential of the force to confine blobs within the nuclear envelope.

The potential V_col was denoted by

where r_i indicates the radius of the i-th blob, $d_{ij}^0$ indicates the basic distance between the i-th and j-th blobs, d_ij=|x_i–x_j|, and θ is the Heaviside step function defined by the following function:

$d_{ij}^0$ was given as the distance between the i-th and j-th blobs in the basic structure of the blob chain if the blobs belonged to the same blob chain and were given as infinity if not. $k_{ij}^c=k^c\sqrt{q_i{q}_j}$ provided the coefficient of the elastic repulsion working when the i-th and j-th blobs came in contact with each other, where q_i indicates the nondimensional parameter of rigidity of the i-th blob dependent on its epigenetic attribute.

The potential V_chr was denoted by

where $k_{chr}^{ad}$ and $k_{chr}^{near}$ are the elastic constants sustaining the shape of the basic structure of each chain. The Σ in Eq. (6) indicates the sum over the set of i-th and j-th blobs that belong to the same blob chain and adjacent to each other along the chain. The Σ in Eq. (7) indicates the sum over the set of i-th and j-th blobs belonging to the same blob chain and coming in contact with each other on the basic structure of the chain but not adjacent along the chain.

The potential V_mem was denoted by

where d_i=|x_i|, R is the radius of the nuclear envelope, and $k_i^m=k^m\sqrt{q_i}$ is the coefficient of the elastic repulsion working when the i-th blob came in contact with the nuclear envelope.

Simulation Method

To simulate the model, the time integral of the Langevin Eq. (1) was calculated numerically using the Eular-Maruyama method with a unit step of 10^–5 s. The parameters of the nucleoplasm features were given as η=0.64 kg m^–1 s^–1, k_BT=4.141947×10^–21 kg m² s^–2 (T=300 K).

Because A-, B-, and M-blobs described the open flexible chromatin regions, the closed dense chromatin regions, and the region with intermediate features, respectively, M-blobs were assumed to be more rigid than A-blobs and less rigid than B-blobs. In the current study, q_i was assumed to be 0.1, 10.0, and 1.0 for A-, B-, and M-blobs, respectively. Other parameters were assumed to be $k^c={10}^{-4} kg s^{-2}$, $k_{chr}^{ad}={10}^{-3 }kg s^{-2}$, $k_{chr}^{near}={10}^{-5} kg s^{-2}$, and $k^m={10}^{-4} kg s^{-2}$. Note that physiologically precise determinations of these parameter values were difficult. Therefore, these values were chosen to be sufficiently large to avoid the crossing of two local blob chains with each other. The qualitative features of the simulation results and conclusions of the current study were not changed, even if the detail factors of these parameters were changed.

The blob chain models of homologous chromosomes were assumed to have the same basic structures, and the blob chains were randomly placed in a spherical shell with a radius of 5 μm for the initial condition of each simulation.

Statistical Analysis of Simulation Data

The probability distributions of the distance between two blobs containing Xic pairs were plotted with a bin interval of 1 μm, consistent with that used in previous reports of experimental results [1,2]. Those of distances between pairs of homologous chromosomes were also plotted with the same bin interval. For both models of ES cells and 2-day cells, 30 simulations with different initial configurations of blob chains were performed. The probability distributions of these distances for each simulation were evaluated using data for blob coordinates from 18 to 36 s. For both models, the average values and 95% confidence intervals (error bars) of the probability of exhibiting each distance were calculated using 30 probability distribution profiles. The significance of the difference in the probability of exhibiting each distance between ES cell and 2-day cell models was examined by Welch’s t-test.

A/B-Compartment Distributions of X Chromosomes in ES Cells Were Unstable but Stabilized during Differentiation

To evaluate chromosome states, A/B-compartments of ES cells and cells 2 days after the onset of differentiation (2-day cells) were evaluated using two biological replicates of Hi-C data drawn from these cell populations (GSM3127755, GSM3127759, GSM3127756, GSM3127760) and GC-content of local chromatin regions of the mouse genome (mm9) [28]. Here, regions corresponding to open chromatin were labeled as A-compartment regions, whereas those corresponding to closed chromatin were labeled as B-compartment regions [37]. In ES cells, the profile of A/B-compartment locations in X chromosomes showed large deviations among the results from two biological replicates, although those of autosomes were robust (Figure 1a). Accordingly, the genomic region in which both replicates exhibited A-compartments or B-compartments was defined as the A-region or B-region, respectively, whereas the genomic region in which different compartments were obtained among the results from two replicates was defined as the M-region.

The M-region was distributed broadly and occupied the largest genomic region among X chromosomes from ES cells (Figure 1b). The epigenetic state of the M-region was expected to change temporally and to exhibit intermediate features of the A- and B-regions on average. By contrast, in 2-day cells, the M-region was dramatically reduced across the entire chromosome, whereas the B-region was increased (Figure 1b). Compared with the X chromosome, the intrachromosomal occupation of the A-, B-, and M-regions did not exhibit large variations in autosomes during this cell differentiation process.

Spatial Distributions of Open/Closed Chromatin Regions on the X Chromosome Changed Following Differentiation

Next, the spatial distributions of A-, B-, and M-regions in each chromosome in ES cells and 2-day cells were measured based on previous results obtained from Hi-C data [28] (https://doi.org/10.5281/zenodo.3371884) and the following procedure: First, the average three-dimensional conformation of each chromosome as the basic conformation was inferred as a polymer by analysis of the Hi-C contact matrices based on a recently proposed consensus method [29–31] implemented in PASTIS as an MDS2 method [32] (see Materials and Methods). Here, since the bin width of the focused Hi-C matrix data was 40 kbp, each monomer composing the polymer was assumed as each local chromatin region with 40 kbp DNA. The position of each monomer in the obtained polymer provided the intrachromosomal spatial position of each region (Supplementary Table S1). Second, the number distributions of monomers with A-, B-, and M-regions as a function of the distance from the center of each chromosome to these monomers (DCC) were evaluated and were, respectively, designated ND_A, ND_B, and ND_M (Figure 2a). ND_α was defined as the number of monomers with α-region between DCC–0.05 (μm) and DCC+0.05 (μm). Third, the radial distributions of A-, B-, and M-regions, designated RD_A, RD_B, and RD_M, respectively, and the radial probability distributions of monomers, RPD, for each chromosome were evaluated as a function of DCC, where RD_X and RPD were defined as RD_X=ND_X/4π(DCC)² (X=A, B, or M) and RPD=(ND_A+ND_B+ND_M)/4π(DCC)²N_C (N_C=the number of monomers in each chromosome), respectively.

The RPDs of chromosomes in both ES cells and 2-day cells exhibited dramatic decreases, with DCC values greater than ~1.2, where the average RPD over all chromosomes showed a significant decrease (DCC≥1.25; Figure 2b, Supplementary Figure S3). Accordingly, monomers showing values of DCC≥1.25 were expected to be localized on the surface of each chromosome.

For DCC≥1.25, RD_B+M–A=RD_B+RD_M–RD_A was measured as a function of DCC (Figure 2c), and the chromatin closedness of each chromosome surface (CL=Σ_DCC>1.2RD_B+M–A) was evaluated in both ES cells and 2-day cells (Figure 2d). From the results, we observed that only X chromosomes in 2-day cells exhibited negative CL values, indicating that the surface of the chromosome was predominantly occupied by open chromatin regions, although other chromosomes were predominantly surrounded by closed or potentially closed chromatin regions.

These results were qualitatively unchanged, even when the basic conformation of each chromosome was inferred using other methods implemented in PASTIS with NMDS, PM1, and PM2 methods [32] (Supplementary Table S2–S4). Therefore, in the following arguments, the polymer constructed by MDS2 method was defined as that of basic conformations of chromosomes.

Simulations of the Coarse-Grained Blob Chain Model Showed Mutual Approach of Xic Pairs in Differentiated Cells

Simulations of models of intranuclear chromosome dynamics in ES cells (ES cell model) and 2-day cells (2-day cell model) were performed using coarse-grained blob chain models. These models were constructed using the previously obtained polymers that gave the basic conformations of chromosomes and the distributions of A-, B-, and M-regions. The populations of monomers successively connected along the polymer were divided into blobs containing local chromatin parts with around 0.1–10 Mbp based on estimations of boundary scores (see Materials and Methods) [28] and appropriate chromosome-dependent characteristic distances along genome sequences. The coordinate and radius of each blob was defined as the centroid coordinate of the blob containing monomers and the standard deviation of the distance from the blob centroid to the monomers, where blob radii were approximately 0.15–1 μm.

Here, the scale of typical blob sizes was similar to that of the persistent length of chromosomes obtained when each chromosome was described as a chain consisting of domains as the units; each domain was defined as a local part of each chromosome divided by the domain boundaries determined by boundary scores (see Materials and Methods). Such chains of blobs provided the basic structure for the coarse-grained model of chromosomes, showing an envelope shape similar to that of the original polymer by elimination of small structural noises expected to be due to the noise of Hi-C data. Additionally, the epigenetic attribute, A, B, or M, of each blob was added based on the epigenetic attributes of monomers contained in each blob (see Materials and Methods) (Figure 3a, b).

Each simulation was performed by considering the force required to maintain the basic structure of blob chains and the repulsion between blobs due to their excluded volumes. The A-blob consisted predominantly of open chromatin regions, which were sparse and flexible, whereas the B-blob consisted predominantly of closed chromatin regions, which contained condensed nucleosomes and proteins. Thus, the rigidity of the blob was assumed to increase in the order A < M < B. The nuclear membrane was assumed to be a spherical shell with a radius of 5 μm (Figure 3c, d).

Thirty simulations were performed for both the ES cell model and the 2-day cell model from different initial configurations, and the time course and probability distributions of the distance between two blobs containing Xic regions were plotted (Supplementary Figure S4, Figure 4a). Significantly higher probability values at small distances and significantly smaller probability values at large distances were obtained in the 2-day cell model compared with those in the ES cell model. This result indicated that two blobs with Xic regions in the 2-day cell model tended to approach each other significantly more frequently than those in the ES cell model, consistent with experimental observations [1–4]. The probability distributions of the distance between the pair of X chromosomes also exhibited similar model-dependent features (Figure 4b). By contrast, for most homologous autosomes except chromosomes 4 and 19, the features of the probability distribution of the distances did not show significant changes between ES cell and 2-day cell models (Figure 4c, Supplementary Figure S5).