Conformational Change of Nucleosome Arrays prior to Phase Separation

Chromatin phase transition serves as a regulatory mechanism for eukaryotic transcription. Understanding this process requires the characterization of the nucleosome array structure in response to external stimuli prior to phase separation. However, the intrinsic flexibility and heterogeneity hinders the arrays’ structure determination. Here we exploit advances in cryogenic electron tomography (cryo-ET) to determine the three-dimensional (3D) structure of each individual particle of mono-, di-, tri-, and tetranucleosome arrays. Statistical analysis reveals the ionic strength changes the angle between the DNA linker and nucleosome core particle (NCP), which regulate the overall morphology of nucleosome arrays. The finding that one-third of the arrays in the presence of H1 contain an NCP invaded by foreign DNA suggests an alternative function of H1 in constructing nucleosomal networks. The new insights into the nucleosome conformational changes prior to the intermolecular interaction stage extends our understanding of chromatin phase separation regulation.


Introduction
Throughout the cell cycle, chromatin undergoes a structural transition between two major states, i.e. a loosely packed transcriptionally active state and a tightly packed transcriptionally silent state 1 . The switching between these two states is achieved by the reorganization of nucleosome arrays. The nucleosome arrays, known as "10-nm" chromatin fibers, consists of tandem repeats 2 of nucleosome core particles (NCPs) 3 separated by 10-60-bp DNA linkers, while each NCP is composed of 147-bp DNA wrapped ~1.7 turns around the histone octamer. Several models have been proposed to explain the conversion between transcriptionally active "10-nm" fiber 4 and its inactive condensed form 5 . For instance, the "10-nm" fiber has been observed to condense into various forms of "30-nm" fibers, including one-start solenoid 6,7 , two-start helical ribbons 8, 9 , untwisted two-start helix 10 , and a crossed-linker 11 . Further hierarchical folding and intertwining of the "30-nm" fibers have been historically considered key intermediates in the attainment of the higher-order chromatin structures 12-14 and chromatids 15, 16 . However, the failure to consistently detect the "30-nm" fiber in vivo [17][18][19][20] has challenged this hierarchical model. In contrast, a reversible thermodynamic process termed liquid-liquid phase separation (LLPS) was proposed to explain chromatin compartmentalization and its condensation behavior. This theory has been supported by evidence from both in vitro and in vivo experiments 21, 22 .
To understand the mechanism of the reversible LLPS process, tetranucleosome arrays (containing the 601-nucleosome positioning sequence, NPS) were frequently used as a standard sample for studying chromatin structures 3,23-26 and their phase separation 21,27, 28 in response to external stimuli, including ionic strength, temperature and protein factors, such as linker histone H1 19,21,27,29, 30 . Recently, we have used the cryogenic electron tomography (cryo-ET) technique to study the early stage of tetranucleosome phase transition 28 (Fig. 1a and Extended Data Fig. 1a-b). Through three-dimensional (3D) reconstruction of tetranucleosome condensates, we described a two-step phase transition mechanism 28 in nucleosomal LLPS. i) A spinodal decomposition process yields irregular condensates via intermolecular interactions among the nucleosome arrays. ii) Those condensates undergo unfavorable conversion into more compact, spherical nuclei, followed by growing into larger spherical aggregates through accretion of spinodal material or by fusing with other spherical condensates 28 . In this process, Histone H1 catalyzes more than 10-fold the spinodal-to-spherical conversion 28 . This study implies that chromatin will naturally maintain or develop into a highly compacted state in physiological conditions unless hydrophobic histone modification or external stimulation are involved 21 .
Cells in their cycle spend 95% of its time in interphase 31 , where large portion of chromatins are loosely packed and accessible in the euchromatin region 32 . However, due to its intrinsic flexibility, a consensus structure model of the opening array is still missing. Thus it is necessary and of great interest to characterize the array structure at a loosely packed transcriptionally active state, and to visualize how the array conformation responds to external stimuli 29 prior to the intramolecular interaction and phase transition. Here, we extended the capability of cryo-ET to snap-shot the 3D structure of each individual particle of mono-, di-tri-and tetranucleosomes at a loosely packed transcriptionally active state, achieved under a low-ionic strength condition (5 mM Na + ) 33 . Through quantitively measuring the conformations, we investigated the conformational change of nucleosome array dynamics affected by two independent factors, i.e. an increased salt concentration (up to 50 mM Na + ) and the presence of linker histone H1. The changes observed could help us understand the reversible process of the LLPS mechanism for gene expression and regulation under external stimuli 29 , especially for those prior to the reported intermolecular interactions and phase transitions.

Flexible mononucleosome array particles
To obtain 3D snap-shot structures of individual nucleosome array particles requires the ability to directly image the NCP and the linker DNA without averaging. Due to its small physical diameter and flexibility, 3D reconstruction of individual DNA fragments is challenging [34][35][36] . Here, we employed both cryo-ET ( Fig. 1b-c) and negative-stain ET (NS-ET) (Extended Data Fig. 1c) to test the 2D-imaging capability to determine the nucleosome positioning along a flexible DNA string. Histone octamers were assembled on a 447-bp-long DNA containing a 601-nucleosome positioning sequence (NPS). The 601-NPS facilitates the formation of highly positioned nucleosome arrays and thus has been used as a standard in a great number of nucleosomal structure studies 3,23-26 and phase separation researches 21,27, 28 . To determine the orientation of the nucleosome, we designed the entry DNA arms to be longer than the corresponding exit DNA arms.
Surveyed cryo-EM and NS-EM micrographs and their representative particles of mononucleosome in low ionic strength buffer (5 mM Na + ) showed a "bead-on-string" structure with characteristic features, including i) a discoidal-shaped NCP with dimensions of ~10 × 10 × 4 nm (Fig. 1b-c) and DNA arms with a width of ~2 nm displaying a ~2 nm helical pitch corresponding to the major groove (Extended Data Fig. 1d-f; ii) nearly 2 turns of DNA wrapped around the histone disc (Fig. 1c, orange arrow); and iii) a low-density pore (~1 nm) present near the center of the nucleosomal disk ( Fig. 1c yellow arrow and Extended Data Fig. 1g-i). These features are consistent with the NCP crystal structure (PDB entry 1AOI) 3 ; likewise, the length of the DNA arms matched their expected values. This test demonstrates our nucleosome 2D-imaging capability, which became the basis for the following 3D reconstructions.

3D structure of an individual mononucleosome
As a proof-of-concept to perform 3D-reconstruction of individual particles, we used mononucleosomes assembled on 447-bp DNA fragments (200-147-100). Reference-free 3D reconstruction by both cryo-ET (Fig. 1d-e) and NS-ET (Extended Data Fig. 2a-c) were performed using a method that combined individual-particle electron tomography (IPET) 37 , image denoising 38 , and missing-wedge correction by low-tilt tomographic reconstruction (LoTToR) 39 . The estimated resolution of up to 43 Å of the resulting 3D reconstruction was estimated based on different criteria and validated by additional analyses (Extended Data Fig. 2d-f). The reconstruction confirmed the "bead-on-string" conformation with structural details, such as the discoidal shape of the NCP (Fig. 1f), the two approximate turns of DNA wrapped around the histone core with defined entry and exit positions (Fig. 1g), and a measurable DNA arm length (Fig. 1f). These features allow us to determine the structure of each individual particle through rigid-body docking of the NCP crystal structure and the flexible fitting of the DNA arms, followed by connecting them together under energy minimization of molecular dynamics (MD) simulation 40 (Supplementary Video 1). Since the structures obtained with NS-ET (Extended Data Fig. 2a-c) were relatively flattened due to the constraints imposed by the thin-layer stain, only the cryo-EM structures were used for further characterization.
The quantitative characterization of cryo-ET structures was conducted by measuring the spatial orientation of the DNA relative to the NCP according to the following parameters: i) wrapping length, i.e., the length of DNA wrapped around the histone core (Fig. 1g, reddashed line); ii) arm/linker angle θ, i.e. the angle formed between two vectors each aligned with 20-bp DNA segments of the entry and exit DNA arms extending from the NCP surface (Fig. 1g, orange arrows); iii) the two orthogonal projections of those vectors through angles, ‖ and ⊥ , with respect to the discoidal plane of the NCP (Fig. 1h); iv) the wrapping angle α of the DNA arm vectors' front-disc projections relative to the corresponding tangents of the discoidal plane (Fig. 1h, left); and v) the bending angle β of the vector's side-disc projections relative to the disc side plane (Fig. 1h right).These parameters provided a basis for the statistical analysis of the nucleosome dynamics.

3D structural dynamics of mononucleosomes
To delineate the dynamics of the mononucleosome 3D structure, a total of 47 density maps were reconstructed by cryo-ET (Fig. 2a, Supplementary Fig. 1 and Video 1), followed by modeling and characterization protocols described above. The ability to fit the DNA turns within the discoidal portion of the maps (Supplementary Fig. 2) enables the alignment of these models, which reveals highly dynamic DNA arms (Fig. 2b). Statistical analysis showed that the DNA wrapped around the histone octamer varied in length (Fig.  2c), consistent with cryo-EM single-particle averaging analysis 26 and reflecting the "breathing" property inherent to NCP 41 . Compared to the crystal structure (1AOI 3 ), the entry and exit DNA arms were unwrapped on average 5 ± 5 bp (mean ± standard deviation) and 11 ± 13 bp, respectively (Fig. 2d), which agree well with MD simulations showing DNA unwrapping in increments of 5 or 10 bp 42 . The NCP exit-side DNA arm displayed a more dynamic behavior than the entry-side, consistent with single-molecule fluorescence resonance energy transfer (FRET) experiments 43 , a property of the 601 DNA sequence with asymmetric affinities, in which the exit-side DNA is less flexible than the entry-side, leading to the former's decreased affinity to the histone octamer.
The angle θ of NCPs are distributed over a wide range from ~0°-140° (with a peak population at ~39°, and a mean ± std of 46° ± 27°), indicating that the two DNA arms were relatively free to orient themselves with respect to one another (Fig. 2e). Its two components ‖ and ⊥ showed their peak populations at ~9° (with a mean ± std of 10° ± 44°) and ~32° (32° ± 45°) respectively (Fig. 2f), indicating that the two DNA arms were relatively parallel to the NCP discoidal plane, and only modestly bent away toward the plane's normal direction. The weak correlation between the entry and exit α angles (Fig.  2g, green points), and between the entry-exit β angles (Fig. 2g, orange points) suggests that the two DNA arms move independently of one another. Thus, we merge entry-and exit-side measurements into a single distribution. Moreover, a very weak correlation between the α and β angles from the same DNA arm (Fig. 2h) implied no preferential DNA arm "swing trajectory" on the NCP surface (Fig. 2i). The α angles displayed a left-skewed distribution in a range of ~± 50° with major and minor peaks at ~39° (39° ± 5°) and ~15° (14° ± 16°), respectively (Fig. 2j), suggesting two low-energy states of the DNA arms along the NCP plane. The β analysis showed a near normal distribution in a range of ~30° with a peak population at ~ -9° (-9° ± 11°) (Fig. 2j), suggesting somewhat restricted dynamics of the DNA arms along the direction normal to the NCP plane.
To display the dynamics of the mononucleosome in solution, the 47 structures were pairwise aligned based on minimization of their root-mean-square deviations (RMSD), followed by ordering them according to a hierarchical classification 44 ( Supplementary  Fig. 3). Through interpolation of the intermediates between each two consecutive ordered structures by target molecular dynamics (TMD) simulations 40 , mononucleosomes are morphed from the most populated state to the rarest state, displaying the thermallyinduced dynamics of the particles in solution (Supplementary Video 1). These results demonstrate the capability of the use of individual particle cryo-ET methods to study nucleosome array dynamics in the following studies.

Effect of NCP number on nucleosome array 3D structural dynamics
To investigate whether increasing the number of NCPs on a DNA segment changes the intrinsic structural dynamics of the resulting array, we repeated the above cryo-ET experiment using di-, tri-and Tetranucleosomes. These arrays were designed in a similar configuration as the mononucleosomes, i.e. 200-bp entry and 100-bp exit DNA arms. Additionally, a 40-bp DNA linker was inserted between two consecutive 147-bp 601 regions. We reconstructed a total of 33, 45 and 31 density maps of the di-, tri-, and Tetranucleosomes, respectively, followed by modelling and quantitative analysis as described above (Extended Data Fig. 3, and Supplementary Fig. 4-8 and Video 2). These nucleosome array structures displayed an overall asymmetric zig-zag architecture ( Fig. 3a-b) consistent with previously described nucleosome fiber conformations found within oligo-nucleosome samples 45 , cell nucleus sections 46,47 , and decondensed mitotic chromosomes 36 , but inconsistent with the symmetric structures previously proposed, such as the twisted double-helix revealed by cryo-EM single-particle averaging methods 24 and crystallography 10,48 .
Unexpectedly, we find that for nucleosome arrays, the number of NCPs is negatively correlated to the θ angle peak position, i.e. ~55°, ~47° and ~41° for di-, and tri-and tetranucleosomes respectively, and that its θ∥ component is more sensitive to the NCP number than the ⊥ component (Extended Data Fig. 5a, row 2-4). The dispersion of θ∥ spans a range of ~33° (from -40° to -21° and -7° for di-, tri-, and tetranucleosomes, respectively), significantly larger than that of ⊥ , i.e. ~3° angle (from -38° to -37° and -35° for di-, tri-, and tetranucleosomes, respectively) (Extended Data Fig. 5a, row 2-4). Among these arrays, the dinucleosomes displayed the most deviated θ and θ∥ peaks compared to the mononucleosome, a feature that may result from the fact that it is the only array structure whose two negatively charged long distal (entry and exit) arms lie in close proximity and repel each other (Extended Data Fig. 5b, panel I, red arrow). In this case, the absence of a third and/or a fourth intermediate NCP separating the two repulsive entry and exit DNA arms induce a unique conformation of dinucleosome, displaying a noticeably restricted dynamic range and a twisting between two NCP discoidal planes (Extended Data Fig. 5b, panel I). By defining θ∥ < 0 as 'closed-arm' conformations (Extended Data Fig. 5c-e), ~83% of DNA linkers were found to adopt the closed forms (Extended Data Fig. 5a). These unique dynamics of dinucleosomes may relate to the ability of remodelers to recognize them. For example, it has been reported that Isw1a prefers to bind dinucleosomes rather than mononucleosomes and its binding affinity to the dinucleosome can be elevated by increasing the distal DNA arm length of the array regardless of histone modifications 49 . Also, SWI/SNF shifts the dinucleosome positions faster than those of trinucleosomes, as shown by electrophoresis assays 50 .
Although the β angle displayed no observable change with NCP number (all peaks were at −9°-−8°), the α angles of nucleosome arrays showed a negatively skewed distribution with two peaks: a major peak at ~32°, similar to the value for mononucleosomes (39°), and a minor peak between −8° and 0° not observed for mononucleosomes (15°) ( Fig. 3c and Extended Data Fig. 6, row 1-4). A large positive α angle corresponds to the DNA conformation that bends away from the NCP discoidal plane, while a small negative α angle suggest that DNA arms could also slightly curve inward after lifting off from the NCP discoidal plane. The shift of the minor peaks from positive to negative values with increasing NCP number likely indicates that the DNA bending at entry and exit positions of the NCP is slightly reduced by the presence of nearby NCPs 51,52 .
Further analysis of the distance between consecutive NCPs (i vs. i+1) among arrays at different lengths showed that NCP spacing is independent of the NCP number and mainly controlled by the distance between the NPS. ( Fig. 3d and Extended Data Fig. 7b, row 1-3). The analysis of the dihedral angle between consecutive NCPs showed no preferred orientation between their discoidal planes except for dinucleosomes, in which the two NCPs preferred a perpendicular arrangement (Extended Data Fig. 7d-e, row 1-3). The lack of a preferred orientation of the former, despite the high torsional rigidity of the DNA (torsional persistence length of ~90 bp), most likely reflect the dispersion of the α and β angles observed in our analysis. The defined orientation observed for dinucleosomes, on the other hand, may reflect the narrower dispersion in the α and β angles for this species, probably as a result of the interaction of the two distal DNA arms. This perpendicular conformation is consistent with the crystal structure of pseudo dinucleosome bound to Isw1a 53 , which suggest, in turn, that the factor recognizes or stabilizes it upon binding.
To display the dynamics of nucleosome arrays in solution, movies of di-and trinucleosomes were also created by the same protocol used for the mononucleosomes (Supplementary Fig. 5 and 7 and Video 2). Above experiments suggest the increased number of NCP on a DNA segment did not dramatically change the intrinsic structural dynamics of the resulting array, except the NCP at two distal ends.

Effect of ionic strength on the 3D structural dynamics of tetranucleosomes
Tetranucleosome has been used as a phenotype with its minimal structural unit 54 and its capability of forming higher-order chromatin structures including helical fibers 24,48 and phase separated condensates 28 . In the previous study 29 we showed, under physiological salt concentration (150 mM Na + and 5 mM Mg 2+ ), tetranucleosomes conduct a phase transition through intermolecular interaction. To investigate how tetranucleosomes initialize the intermolecular interaction prior to the early stage of phase transition 29 , we repeated the cryo-ET experiment described above by slightly increasing the salt concentration to 50 mM Na + (Extended Data Fig. 8a), and achieved a total of 34 density maps with their fitted models ( Fig. 4, Supplementary Fig. 9 and Video 3). These structures again confirmed the asymmetric zig-zag arrangement of the arrays but displayed a more compacted shape in contrast with the structures observed at 5 mM Na + .
The extent of DNA unwrapping at the entry-side remains invariant with ionic strength, i.e. -4 ± 9-bp vs. -4 ± 8-bp for 50 mM and 5 mM, respectively, but the interaction at the exitside was slightly stabilized by the increase of Na + (-9 ± 8-bp vs. -14 ± 12-bp, for 50 mM and 5 mM respectively) (Extended Data Fig. 4b, row 4-5). Although θ angle is not sensitive to the increase in ionic strength (changing its peak position from 41° to 47°), its two perpendicular components θ∥ and ⊥ display a significant peak shift with the increase of Na + (from -7° to -42° and 35° to 19°, respectively) yielding a high percentage of closedarm conformations (Extended Data Fig. 5a, row 4-5). Although the major peak of α remained unchanged, its minor peak became more apparent and shifted from -8° to -14° ( Fig. 4b and Extended Data Fig. 6c, row 4-5). This transition may arise from the stronger electrostatic screening effect 55 , which reduces the repulsive interactions between the two DNA linkers, allowing them to become closer together. Consistently, the DNA linkers bend toward the NCP discoidal plane (due to the change of β from -9° to -4°, Fig. 4b and Extended Data Fig. 6d, row 4-5). Although further analysis showed a slight reduction in the distance between consecutive NCPs (i vs. i+1), (from 237 Å to 218 Å) with increased ionic strength, a significant distance reduction was found between every-other NCPs (i vs. i+2, from 265 Å to 114 Å), and between the first and fourth NCPs (i vs. i+3, from 388 Å to 269 Å) ( Fig. 4c and Extended Data Fig. 7b, row 3-4), indicating a compaction process of the overall array conformations (Fig. 4a). Meanwhile, the dihedral angle between NCPs showed that, although the consecutive NCPs adopted a more perpendicular orientation with respect to each other, every-other NCPs turned to more parallel packing in contrast to the observations at low salt concentrations (Extended Data Fig. 7e, row 3-4). This experiment suggests the ions increased the electrostatic effect by inducing two DNA linkers to get closer to each other, resulting in the NCP packing towards to each other, resulting in shrinking the overall diameter.

Effect of H1 on the 3D structural dynamics of tetranucleosomes
Considering the Histone H1 catalyzes the spinodal-to-spherical conversion 28 for more than 10-time 29 faster than that absents H1, it is interesting to investigate how H1 alone contributes to speeding-up the conversion during the regulation of the structural dynamics of the arrays. We incubated tetranucleosomes with H1 in a molar ratio of 1:4 under the same low-salt conditions (5 mM Na + ). A total of 33 cryo-ET density maps with their models were obtained (Extended Data Fig. 8b, Supplementary Fig. 10 and Video 3). Once again, these arrays displayed a general asymmetric zig-zag conformation (Fig. 4d). However, these structures seem partially intertwined compared to arrays observed under low-salt without H1 (Fig. 3a), but less compacted than those imaged under high-salt without H1 (Fig. 4a). Interestingly, ~40% of the particles displayed a conformation in which the 200-bp DNA arm invades one of the NCPs by replacing part of its originally wrapped DNA at the exit side (Fig. 4e). This invasion can be observed within an array as well as between arrays, leading to the formation of larger structures (Extended Data Fig.  9a), reminiscent of the previously reported alternative function of H1 in strengthening the inter-array network 56 .
The analysis showed that the presence of H1 did not change the amount of DNA wrapping at the entry side (-6 ± 8-bp vs. -6 ± 8-bp) but slightly increased the unwrapping level and dynamics on the exit side (-17 ± 19-bp vs. -14 ± 12-bp) and generated a second peak at ~-45-bp (Extended Data Fig. 4a-b, row 6). This change in length implies that the DNA entirely unwrapped from the H2A-H2B heterodimer and reflects a weakened DNA-histone binding 57 leading to a spontaneous site exposure 58,59 .
The distribution of θ angles showed a slight shift of the peak position from 41° to 50° in the presence of H1 (Extended Data Fig. 5a, row 4 and 6, left). However, its two components display more significant changes. First, the major peak of θ∥ shifts from -7° to -48° with H1, which was similar to the shift observed in 50 mM Na + (from -7° to -42°) (Extended Data Fig. 5a, row 4-6, middle). Second, a minor peak of θ∥ emerged at 136° in 28% of the NCPs, adopting an alternative "wide-open arm" conformation (Extended Data Fig. 5a, red arrow). Third, the major peak of θ ꓕ showed the DNA linkers adopting a more parallel arrangement when viewed from the side of the NCP in the presence of H1 (from 35° to 6°), similar to the change caused by the increase of Na + (from 35° to 19°) (Extended Data Fig. 5a, row 4-6, right). Finally, a second peak of θ ꓕ at 152° suggests the co-existence of alternative conformations (Extended Data Fig. 5a, green arrow). Classification analysis (Extended Data Fig. 9b, orange arrows) confirmed that these minor peaks mainly resulted from DNA unwrapping around the second NCP (Fig. 4e, II).
The α and β angles show the major population of DNA linkers turned toward each other in the presence of H1, likely due to the "triangular-shape" formed between the NCP and its two linkers (Fig. 4e, red triangles) (i.e. the major and minor peak of α decrease from 32° to 27° and -8° to -16°, respectively, while the β peak decreases from -9° to -1°, Fig.  4f and Extended Data Fig. 6c-d, row 4 and 6). These changes were again, similar to those induced by 50 mM Na + (Extended Data Fig. 6c-d, row 5). Interestingly, these two different conditions modulate the DNA linkers dynamics on the NCP in a similar manner even though one involves electrostatics screening and the other protein binding 10, 25 .
Characterization of the nearest-neighbor distance of NCPs shows that the major peak of consecutive NCPs was slightly decreased by H1 (from 237 Å to 205 Å), an effect also seen by the increase of Na + (from 237 Å to 218 Å) (Extended Data Fig. 7b, row 3-5, left). The second peak appearing at 313 Å (Fig. 4g, left) likely reflects co-existing distal arm invading conformations (Fig. 4e), in which the linker between i and i+1 NCPs was extended due to the DNA unwrapping. Meanwhile, the major peak depicting the distance between i and i+2 NCPs only decreases from 236 Å to 213 Å compared to the reduction caused by 50 mM Na + (from 236 Å to 114 Å) (Extended Data Fig. 7b, row 4-5, middle), and the second peak moved even further to a distance of 453 Å, which indicates a volume expansion of the array instead of a compaction. This result is also confirmed by the distance between the i and i+3 NCPs (Extended Data Fig. 7b, row 4-5, right).
The analysis of NCP dihedral angles showed that H1 did not induce a preferred orientation among NCP discs compared to those acquired under 50 mM Na + (Extended Data Fig. 7d-e, row 5). Combined with the angle/distance analysis, these observations suggest that H1 can only partially regulate the arrays conformation, such as by changing the direction of DNA linkers on NCPs, but not the interaction among NCPs. The reason may be as follows: the presence of H1 bridges locks the two DNA linkers/arms of an NCP, yielding a similar outcome as high salt concentrations by drawing them closer together. However, H1 is not able to reduce the electrostatic repulsion among the charged NCPs due to the high ionic strength via charge screening. As a result, H1 alone at low-salt conditions cannot induce the full compaction of the array. The appearance of conformations in which the upstream DNA arm of the first nucleosome invades and wraps around the second nucleosome (Fig.4e, panel II) may reflect the fact that, at low ionic strength, NCPs are kept separated at "preferred" distances leading in some cases to partial unwrapping of the DNA, which in turns permits the wrapping of "foreign DNA". It also leads to additional DNA junctions that provide other possible binding sites for H1 (Fig. 4e, II and III, black arrow), which could further stabilize the system. It is worth noting that the invasion of the second NCP is the dominant conformation observed, compared to those of the third and fourth NCPs, which yield two dinucleosome regions separated by a longer linker (Fig. 4e, II, red-dashed line). This experiment suggests H1 reduced the freedom range of two DNA linkers as reported 60 , while increasing the connectivity among the NCPs but not significantly reducing the overall diameter.

Modelling the effect of linker DNA orientation and dynamics on chromatin morphology
To validate whether the angle change between two DNA linkers can really modify the large-scale morphology and density of chromatin fibers, hecta-nucleosome arrays were computationally generated based on the obtained statistics. The analysis of simulated hecta-nucleosomes at 5 mM Na + showed curvy irregular fibers ( Fig. 5a and Extended Data Fig. 10a) with an average length and width of 561 ± 44 nm and 77 ± 19 nm, respectively (Extended Data Fig. 11a, row 1). By using the parameters obtained at 50 mM Na + , both the chromatin length and width were reduced to 395 ± 35 nm and 38 ± 16 nm, respectively (Extended Data Fig. 11a, row 2), producing irregular/disorganized and condensed "30 nm fibers" (Fig. 5b and Extended Data Fig. 10b) as observed experimentally 61,62 .
To simulate possible chromatin condensates, the generated fibers were inserted into a 200 nm cube with two different approaches: "seamless assembling" vs "random docking" (for details see methods). Results showed that the estimated density of condensates were between 134 -266 μM at low-salt conditions. These values were nearly doubled (263 -341 μM) at 50 mM Na + conditions (Fig. 5c-d and Extended Data Fig. 11b, row 1-2). These chromatin densities are within the experimentally measured range of 80-520 μM in vivo 30 and ~340 μM in vitro 21,28 .
The 30-nm fibers generated in this simulation were not often observed in cryo-EM imaging 17,63 or cryo-ET cell sectioning slices 32,46 . One possible explanation for this discrepancy is the low probability of capturing fibers lying exactly on the plane of the thin slice cell section, especially when the latter approaches a thickness of 30 nm 12, 46 . Slicing through the simulated condensates across the z-axis with a 25 nm-thick slab greatly reduces the fiber pattern and replaces them with more irregular structures (Fig. 5c-d

left panel).
Because NCP "breathing" dynamics plays a critical role in biologically relevant DNA processes, such as enabling restriction enzyme access to the inner wrapped DNA target site 58 and nucleosome translocation without the participation of chromatin remodelers 64 , we repeated the above simulations by explicitly preventing DNA unwrapping. The resulting simulations increased the length of the chromatin fibers by 11-13% (from 561 nm to 646 nm at low-salt, and from 395 nm to 441 nm at high-salt) (Extended Data Fig.  10c-d and 11a-b, row 3-4), but did not induce an obvious change of the condensate's density in either salt conditions (Extended Data Fig. 10e-f), suggesting that breathing does not increase enzyme accessibility to the interior of the condensates. Instead, the DNA linker angles of the arrays are found to have a major effect on the density of condensates and the morphology of chromatin.
Unlike conventional averaging methods, in which a population of particles possessing homogeneous structural identity is selected from a large and diverse pool of particles in thermodynamic equilibrium, the individual-molecule approach of cryo-ET provides continuous particle conformational distributions by randomly sampling the system. In this case, the probability ( ) of finding a molecule in an energy level ( ) is proportional to − / based on the Boltzmann distribution, where is the Boltzmann constant and T is the temperature. Thus, the energy difference (∆ = − ) between two energy states can be calculated from ∆ = − = − ( / ) × . In this study, we observe two peaks of α angles in all the molecular populations studied (Extended Data  Fig. 11c), which correspond to two energy states. The calculated energy difference between the relatively opened and closed DNA arm conformation were 0.17 (given the corresponding populations of 54% vs. 46%), 0.27 (43% vs. 57%), 0.56 (36% vs. 64%) and 0.43 (39% vs. 61%) for mono-, di-, tri-, and Tetranucleosome arrays, respectively. However, the energy differences changed to 0.48 at 50 mM Na + (38% vs. 62%) and 0. 69 in the presence of H1 (33% vs. 67%). Since the bending energy of wrapped nucleosomal DNA is in a range of 0.34 to 0. 46 /bp 65 , the transition between the two conformations requires only the bending of one or two DNA bp per nucleosome, which results in the observed large-scale array conformational changes. These simulations confirmed the angle between two DNA linkers plays a key role in regulating the overall morphology and density of chromatin fibers.

Conclusions
The structural dynamics and conformational changes of loosely packed nucleosome arrays are critical to our understanding of transcriptionally active chromatin architecture prior to the formation of compacted nucleosomal condensates. To gain insights into the nucleosome array structure in the active state and their conformational change in response to chemical and biological factors, we have determined the conformation for each individual nucleosome array particle (without averaging) by cryo-ET. We found that an increased ionic strength reduced the angle between the DNA linker and the NCP, which resulted in a decreased size and increased density of the asymmetric tetranucleosome array unit. In comparison, H1 likewise changes the dynamics of the array but builds alternative nucleosomal networks. The surface physical properties of NCP, such as polarization and hydrophobicity, may be changed by altering the geometry of the nucleosome array, which may then strengthen the intermolecular interaction for the formation of a spinodal structures. The regulation of gene expression followed the development of the spinodal structure into compacted spherical condensates, which may serve as the underlying mechanism of how chromatin mediates the change from interphase to metaphase structures. The new insight into these conformational changes prior to the nucleosome arrays intermolecular interactions extends our knowledge to understand the regulation mechanism of nucleosome arrays via phase separation.

Expression and purification of histones
Vectors containing the genes of Xenopus laevis histones H2A, H2B, H3, and H4 under the control of a T7 promoter were expressed in E. coli BL21(DE3) and purified as previously described 66 . For each histone, a single E. coli BL21(DE3) colony was grown at 37˚C in LB media supplemented with 100 µg/uL ampicillin and 25 µg/uL chloramphenicol to an OD600 of ~0.6, and histone expressions were induced by adding 1mM Isopropyl b-D-1-thiogalactopyranoside (IPTG) to the cell culture. After 3 h, the cell culture was centrifuged at 7,000 × g for 30 min, and the cell pellet was suspended in icecold wash buffer (50 mM Tris-HCl pH 7.5; 100 mM NaCl; 1 mM EDTA; 5 mM 2mercaptoethanol (BME); 1% Triton X-100 [w/v]; and protease inhibitors, Roche). This procedure was repeated twice to wash away LB medium. The cell pellets were suspended in 4 volumes of wash buffer, flash freezing with liquid nitrogen, and stored at -80˚C for later use. To purify the inclusion bodies, suspended cell pellets were thawed and sonicated on ice seven times with 20-s bursts at 7.0 W with a Misonix 2000 sonicator.
The lysate was centrifuged at 30,000 x g for 1 h at 4˚C, and the pellets were rinsed by suspension with wash buffer. This procedure was repeated three times. To remove Triton X-100, pellets were rinsed by suspension using washing buffer without the detergent. Suspended pellet was centrifuged at 30,000 x g for 1 h at 4˚C and this procedure was repeated three times. Inclusion bodies were then solubilized in 20 mM Tris-HCl pH 7.5; 8 M Urea; 1 mM EDTA; 10 mM DTT, and purified by anion and cation exchange. Purification was checked by 15% SDS-PAGE and fractions containing pure histones were pooled and dialyzed against 1 L of 10 mM Tris pH 8.0. Histones were centrifuged to remove aggregates, concentrated (∼10 mg/mL) using a 10K Amicon Ultra-15, lyophilized, and stored at -80°C for later use.

Histone octamer purification
The reconstitution of the X. laevis histone octamer was performed as previously described (Wittmeyer). Briefly, lyophilized histones were solubilized with unfolding buffer (20 mM Tris-HCl pH 7.5; 7 M guanidine hydrochloride; 10 mM DTT), and combined at a H2A:H2B:H3:H4 ratio of 1.2:1.2:1:1. The volume was adjusted to a total histone concentration of 1 mg/mL and it was incubated with mixing for 3 h at room temperature. Solubilized histones were dialyzed four times for 12 hours each against 1 L of refolding buffer (10 mM Tris-HCl pH 8.0; 2 M NaCl; 1 mM EDTA; 5 mM DTT) using a 3.5 kDa dialysis membrane at 4°C. Refolded histone solution was centrifuged at 100,000 x g for 30 min to remove aggregates, concentrated by centrifugation to ∼0.5 mL using 10 kDa Ultra-15 (Milipore), and loaded onto Superdex 200 Increase 10/300 GL (Cytiva) previously equilibrated with refolding buffer. Fractions were checked by 15% SDS-PAGE and the gel was stained with AcquaStain protein staining (Bulldog Bio). Fractions containing three bands (corresponding to H3, H2A/H2B, which, because of their similar size they are not resolved in the gel, and H4) and in equimolar quantities (based on the intensity of the three bands) were pooled, concentrated to ∼8-10 mg/mL using 30 kDa Ultra-15 (Milipore), aliquoted, and flash-frozen with liquid nitrogen, and stored at -80°C for subsequent use. The synthesis of the H2A-H2B dimer followed the same procedure utilized for octamer reconstitution, with H2A and H2B combined at a ratio of 1:1.

Synthesis of DNA templates
The mono-, di, tri-, and tetranucleosome DNA templates consists of one, two, three, and four repeats of the 601-nucleosome positioning sequence (NPS), respectively. Each template is flanked to the left by 200 bp DNA and to the right by 100 bp DNA, and the 601 NPS in di-, tri-, and tetranucleosome templates are separated by a 40 bp linker length. Mononucleosome DNA templates were amplified by PCR from the pGEM-3Z/601 plasmid (addgene) and cloned back into the PGEM 601 vector using primers containing the restriction recognition site for the BsaI enzyme (NEB). Di-, tri-, and tetranucleosome DNA templates were synthetized by ligation of two, three, or four PCR products containing one 601 NPS each. The PCR products were amplified using the pGEM-3Z/601 plasmid (Addgene) and primers containing the restriction recognition sites for enzymes BsaI and BbsI (NEB). In our design, the di-, tri-, and tetranucleosome DNA are flanked by BsaI restriction sites, and the BbsI sites were used to ligate the two, three, or four PCR products, respectively. To generate the full templates, the PCR products were digested first with BbsI and ligated at an equimolar ratio using E. coli DNA ligase (NEB). The longest product of each ligation, corresponding to the di-, tri, or tetranucleosome DNA, were purified using 0.8% agarose gels, digested using the BsaI enzyme, and cloned into BsaI-restricted pGEM-3Z/601 plasmid using T4 DNA ligase (NEB). Ligation product was transformed into E. coli DH5α for plasmid extraction by miniprep. Each ligation product was checked by DNA sequencing. Plasmid containing the mono-, di-, tri-, and tetranucleosome templates were grown in dam -/dcm -E. coli (NEB), purified by maxiprep, and excised by restriction with BsaI. Mononucleosome DNA template was purified from the vector backbone by 5% preparative acrylamide electrophoresis using a Model 491 Prep Cell (Bio-Rad) electroelution system. Di-, tri-, and tetranucleosome templates were purified by 4% preparative acrylamide electrophoresis.

Nucleosome assembly and purification
X. laevis histone octamer and mono-, di-, tri, and tetranucleosome DNA templates were combined at a ratio of 1:1.2, respectively, in high-salt buffers (10 mM Tris-Cl pH 8.0; 2 M NaCl; 1 mM EDTA; 0.5 mM DTT; and 1 mM PMSF) at a final DNA concentration of 100 ng/uL. In the case of the arrays, H2A-H2B heterodimer dimer was also incorporated at a molar ratio of 0.2 compared to the octamer. Assembly solutions were dialyzed using a 3.5 kDa dialysis membrane at 4°C against 500 mL of high-salt buffer for 1 h at 4°C, followed by a 36-hour lineal gradient dialysis against 2 L of low-salt buffer (10 mM Tris-Cl pH 8.0; 1 mM EDTA; 0.5 mM DTT; and 1 mM PMSF) with continuous stirring. A final dialysis of 3 h was performed against 500 mL of low salt buffer. Nucleosome reconstitutions were checked by 4% acrylamide and 0.2X TBE buffer (Tris-borate-EDTA) native electrophoresis. Mononucleosomes were purified from hexasomes and bare DNA by 4% preparative acrylamide (59:1 acrylamide:bisacrylamide) electrophoresis using the 491 Prep Cell (Bio-Rad). Di-, tri-, and tetranucleosomes were purified by 10-40% lineal sucrose gradient (20 mM HEPES-NaOH pH 7.5; 1 mM EDTA; 1 mM DTT) and centrifuged for 16 h at 38,000 rpm at 4°C, using an ultracentrifuge Beckman Optima MAX-XP with the rotor MLS-50 (Beckmann), as previously described (Mol cell paper reference).

TEM specimen preparation
The cryo-EM specimens were prepared following the procedure described before 37 .
Briefly, an aliquot (~ 3 μl) of nucleosome array sample ~400 nM was placed onto the 200 mesh Quantifoil copper grid (Q210CR-06, Electron Microscopy Sciences) that had been glow-discharged for ~15 s by (PELCO easiGlow™ Glow Discharge Cleaning System) for 15 seconds. After incubating for ~ 10 sec, the grid was flash-frozen in liquid ethane at ~ 90% humidity and 4 °C with a Leica EM GP rapid-plunging device (Leica, Buffalo Grove, IL, USA) after being blotted with filter paper with a controlled blotting time (2 s). The flashfrozen grids were transferred into liquid nitrogen for storage.
The NS-EM specimens of nucleosome array sample were prepared using the optimized negative-staining protocol (OpNS) as described 67 . In brief, the samples were diluted to ~20 nM with sample buffer. An aliquot (~4 μL) of diluted sample was placed on an ultrathin carbon-coated 200-mesh copper grid (CF200-Cu-UL, Electron Microscopy Sciences, Hatfield, PA, USA) that had been glow-discharged for ~15 s. After ~1 min incubation, the excess solution on the grid was blotted with filter paper. The grid was then washed with water and stained with 1% (w/v) uranyl formate (UF) before air-drying with nitrogen.

TEM data acquisition
Cryo-EM specimens were screened by a Titan Krios (FEI) transmission electron microscope operated at 300 kV high tension with a Gatan energy filter. The untitled cryo-EM micrographs were acquired with a Volta phase plate under defocus at ~2.5 μm using a Gatan K2 Summit direct electron detection camera under a magnification of ~81 kx (each pixel of the micrographs corresponds to ~0.9 Å in specimens). The particles with a box size of 800 × 800 pixel will be used to display the morphology. Cryo-EM tilt image series of the samples were collected from -60° to +60° at 3° increments on a Titan Krios G2/G3 TEM equipped with a Gatan energy filter and a K2 Summit direct electron detection camera operated under 300kV high tension. During data acquisition, the SerialEM 68 software was used to automatically track the specimen and maintain defocus at ~2.5 -3.0 μm. The acquired tilt image series at magnification of ~81/50 kx (each pixel corresponds to 1.46/1.47 Å for K2/K3 camera) represents a total dose of ~102 -183 e-/Å 2 . For each tilt angle, a total number of ~8 frames were collected under the exposure of 0.25 s per frame.
NS-EM specimens were screened by using a Zeiss Libra 120 Plus TEM (Carl Zeiss NTS) operated at 120 kV high tension with a 10-20 eV energy filter. The OpNS micrographs were acquired under defocus at ~0.6 μm using a Gatan UltraScan 4K × 4K CCD under a magnification of 125 kx and 160 kx (each pixel of the micrographs corresponds to ~0.94 Å and ~0.74 Å in specimens respectively). Tilt image series of mononucleosome samples were collected from -60° to +60° in 3° increments using a Zeiss Libra 120 Plus TEM (Carl Zeiss NTS) equipped with an in-column energy filter and a Gatan UltraScan 4K X 4K CCD. During data acquisition, the Gatan tomography module (Gatan Inc., Pleasanton, CA, USA) operated in Advanced Tomography mode was used to track the specimen and maintain defocus at ~2.0 μm. The acquired tilt image series at magnification of 50 kx (each pixel corresponds to 0.24 nm) represents a total dose of ~60 e -/Å 2 .

Image preprocessing
The motion of the cryo-ET frames was corrected by MotionCor2 69 . To reduce the cryo-ET image noise, we followed a machine learning method (NOISE2NOISE, T2T method) as described 38 . The contrast transfer function (CTF) of both cryo-ET and NS-ET tilt series was measured by using ctffind3 software 70 and the phase and amplitude were corrected by SPIDER 71 , GCTF 72 and TOMOCTF 73 after the X-ray speckles were removed. The tilt series were initially aligned by using IMOD 74 . By using boxer from EMAN software 75 , a box of 256 × 256 pixels was used to select the particles of nucleosome from the ~113 micrographs imaged under 125 k× magnification, while a box of 200 × 200 pixels was used to select the segments of the upstream (or entry) and downstream (or exit) DNAs from the ~150 micrographs imaged under 160 k× magnification. All NS-EM particles were masked using a round mask generated from SPIDER software after a Gaussian highpass filtering. The reference-free class averages of particles were obtained using refine2d (EMAN software) based on 30,540 particles of DNA segment and 13,029 particles of NCPs.

Individual particle electron tomography (IPET) 3D reconstruction
The tilt series of each of the targeted particles were semi-automatically tracked and then windowed in square windows of ~ 1,000 × 1,000-pixel size using IPET software 37 , before being binned four to five times to reduce computation time in subsequent reconstructions. Following the pipeline of IPET reconstruction 37 , a local tilt series images containing a single nucleosome particle was extracted from the IMOD-aligned full-size tilt series to perform focused ET reconstruction (FETR), which can eliminate the effects from the large image distortion 37 . Briefly, an ab initio 3D density map was directly back-projected in Fourier space and served as an initial model. The refinement was then iteratively invoked to translationally align each tilted particle image to the computed projection. During the refinement, automatically generated Gaussian low-pass filters, soft-boundary circular masks and loose particle-shaped masks were sequentially applied to the tilt images and references to increase the alignment accuracy 37 . An improved model was then reconstructed based upon the refined alignment at the end of each refinement iteration. The 3D map was then reconstructed by back-projection of the filtered and masked particle tilt series. The final back projection was performed in Fourier space without weighting. To reduce the artifact caused by the limited tilt angle range, the final 3D map was submitted for a published missing-wedge correction method by LoTToR 39 , in which the low resolution mask used was generated by the Model-Based Iterative Reconstruction (MBIR) method 76 . All final density maps were low-pass filtered to 4.5 nm using EMAN 75 and displayed using UCSF Chimera 77 .

Estimation of the reconstruction resolution
The resolution for the reconstructions were estimated by two methods. i) Data-to-Data based analysis: the Fourier Shell Correlation (FSC) was calculated between two independently reconstructed 3D maps, in which each map was based on one-half of the tilt-series (split by even and odd tilt index) after particle alignment refinement of the IPET 37 . The frequencies at which the FSC curve first falls to values of 0.5 and 0.147 were used to represent the reconstruction resolution. Notably, the resolution estimated by this method could be severely under-estimated since the reconstruction from one half of the tilt-series significantly reduced the quality of the map compared to the final reconstruction. ii) Data-to-Model based analysis: the FSC curve between the final IPET reconstruction and the density map converted from the corresponding fitting model was calculated. The frequencies at which the FSC curve fell below 0.5 was used to estimate the resolution. The density map of the fitting model was generated by pdb2mrc in EMAN software 75 .

Modeling the structure of nucleosome arrays
To build a structural model for the reconstructed map of each individual nucleosome array particle, the pathways for the two flanking DNA arms were initially traced by sampling a group of 3D points located at the high-density loci of the map followed by sorting their order into a points list. Then, by fitting the crystal structure of Xenopus laevis NCP into the discoidal-shaped high-density region, the DNA pathways with the NCP were also defined and converted to a list of 3D points. These points series from different models were merged together after removing adjacent clashes and then fitted with a smooth quadratic Bezier curve followed by conversion into DNA model by using GraphiteLifeExplorer 78 . The total length of the DNA used to thread the model matched with our designed DNA construct, i.e. 456, 632, 818, 1008 for mono-, di-, tri-and tetranucleosome arrays, respectively. Notice that for some of the maps, a small portion of DNA densities (~10%) was missing near the middle portion of DNA arms (which may be due the orientation of DNA segments aligned near-perfectly perpendicular to the beam direction, resulting the lowest image contrast of the DNA). Fortunately, those missing portions were small and did not prevent the DNA model fitting, which can be circumvented by interpolating surrounding density at those loci. To further refine the model, the Molecular Dynamics Flexible Fitting (MDFF) 40 was applied to energy minimize the model under a force gradient created by the electron density map.

Evaluating the fitting models
To evaluate the self-consistency of the above fitted models, the correlation between the measured DNA linker length between two consecutive NCPs (named as ( )) and the estimated DNA unwrapping level between the same NCPs were calculated. The ( ) was acquired from the experimental data, in which distance between the DNA exit position of the nth NCP and the entry position of the (n+1) th NCP from the density map were measured. The DNA unwrapping level was estimated from the fitted model, in which the unwrapping angles from the exit side of nth NCP and entry side of the (n+1) th NCP were added together ( ( ) = ( ) + ( + 1)). The Pearson correlation coefficient, r=-0.9 was calculated from the linear regression fitting of the unwrapping length ( ) against the angle ( ) using the statsmodels package in Python. The off-template NCPs (binding to the 100 or 200-bp distal arms region) are removed from the following statistical analysis.
Defining the entry and exit linker DNA origins on the NCP By using the fitted nucleosome array models, the entry and exit DNA arm origins on the NCP can be estimated with the following procedure. The DNA portion of the fitted nucleosome array model was first converted into a list of 3D points (series m) by averaging the coordinates of C1 atoms of each base-pair. After aligning the nucleosome crystal structure to the fitted array model at its nth NCP region, a points list (series n) along the wrapped DNA from the nucleosome crystal structure was also generated by the same method. By comparing the two points lists, the overlapped region on the fitted DNA model can be identified if any points from the series n were found within its 8 Å radius of series m. The 8 Å criterion was chosen based on that, when two DNA helical centers separated away from each other for over one third of the DNA diameter (~24 Å), the separation of two aligned DNAs can be distinguished. The base pair indexes at the two distal ends of this overlapped region were used to define the entry and exit DNA linker arm origin along the array. The total number of the base-pairs between the entry and exit position were used to calculate the length of wrapped DNA on the histone surface.

Measuring the NCP wrapping dynamics on histone surface
By measuring the position and length of the wrapped DNA on each histone surface, the DNA unwrapping footprint along the array can be quantified as follows. For a fitted array model, a binary score was assigned to each DNA base pair along the DNA sequence depending on its contact with the histone octamer ("1" for contact if the atoms of the base pair fall in a 4 Å distance with any histone octamer atoms and "0" for non-contact). By averaging the scores of the i th base pair within same type of fitted array models, the mean score distributions along the DNA template (from upstream to downstream) for a mono-, di-, tri-, and tetranucleosome array were calculated separately. The score distributions around the regions contains the designed 146-bp Widom 601 position sequence 79 reflects the probability of finding DNA unwrapping events. The difference of the score distribution for the entry and exit side of each NCP loading position reflects the asymmetry property of the DNA sequence in binding to histone octamers and the dynamic relationship between DNA unwrapping and equilibrium assembling, such as transition states among the tetrasome, hexsome and full octasome 80 .

Defining the entry and exit linker DNA vectors on the NCP
Due to the fact that the 20-bp DNA segment (~6.8 nm) is relatively stiff compared to the persistence length of DNA ~50 nm, the base pairs residing at the entry side DNA arm origin and its 20-bp upstream were used to define the start and end-point of a vector, respectively (based on their helical center coordinates). This vector was used to represent the entry DNA arm pointing direction. Similarity, the base pairs residing at the exit side DNA arm origin and its 20-bp downstream were used to define the exit DNA arm vector.

Measuring the angle θ, wrapping angle α and bending angle β of the DNA arm vectors
The angles of the DNA linker arms that extended from the discoidal-shaped NCP surface were measured by the following steps. i) Defining the X-, Y-, and Z-axes of each NCP model. The Z-axis was defined along the helical axis of the wrapping DNA measured from its rotational symmetry. The Y-axis was defined by the dyad axis of the NCP. The crosspoint of these two axes was used as the center of NCP and defined the X-axis, where it simultaneously passed the NCP center and was perpendicular to both Z-and Y-axis. ii) Defining the relative DNA arm angle θ. θ measured the angle between the previously defined entry and exit DNA linker arm vectors on the same NCP. Because DNA arm conformational states ("open and close") must be defined relative to the NCP, the projections of the θ angle on the NCP discoidal plane (X-Y plane) and its perpendicular plane (Y-Z plane) were measured as θ∥ and θꓕ, respectively. iii) Defining the in-plane wrapping angle α of an NCP arm vector. α calculated the angle formed by two vectors within the X-Y plane. One vector is the projection of the NCP linker arm vector on the X-Y plane, and the other vector is defined by the tangential direction of the discoidal-shaped projection of NCP on X-Y plane, which crosses the origin of the corresponding NCP linker arm. iv) Defining the out-of-plane bending angle β of an NCP arm vector. β calculated the angle formed by the NCP arm vector and the X-Y plane. The measured angle distribution was fitted with either one gaussian or two gaussians with the sklearn.mixture.Gaussian Mixture package.
Measuring the intra-array NCP core-core distances and plane-plane angles To quantitatively define the spatial relationship among the NCPs from different types of nucleosome arrays, the core-core distances and plane-plane angles between each pair of NCPs were measured. The core-core distance, ( , + ), measured the distance from the center of the nth to the (n+m) th NCPs. The core-core angle, ( , + ), measured the dihedral angle between the discoidal planes (X-Y planes) of the nth and (n+m) th NCPs. The histograms were fitted by a Kernel Density Estimation (KDE) function 81 . The dihedral angle distributions were compared with a sine function, which represents a distribution of orientations of two random planes (the angle θ between their normal directions), which in turn are equivalent to the probability of finding a point on the unit hemisphere contained in a differential ring-shape area:

Insilco-construction of the chromatin-like higher order structure
The distribution of the wrapping and bending angles relative to the NCP were used to build longer in silico nucleosome array fibers. Each longer array fiber containing 100 NCPs was generated by randomly connecting nucleosome model units. To be specific, due to the large number of atoms within the NCP (>12k) and the fact that the length of 40-bp DNA (~14nm) is much smaller than the persistence length of DNA (~50 nm), NCP units were coarse grained where the histone octamer and linker DNA arms were treated as spheres and straight lines, respectively. Based on the distribution of the measured mean and standard deviation of the bending and wrapping angles (i.e. α and β) for the tetranucleosome array at low salt and high salt condition, two vectors in a 20-bp length followed the same angle distributions relative to the NCP were randomly generated and used to represent the entry and exit DNA linker arms by using VMD software 82 . By repeating this process, a pool of NCP units with various linker arm conformations but following the same spatial distribution to the experimental data were prepared. To assemble the NCPs units into a fiber, randomly selected NCPs from the pool were sequentially connected together based on the following procedures: i) The exit side 20bp linker for the ith NCP was linearly connected to the entry side 20-bp linker for the i+1th NCP, ii) The 40-bp linker causing ~70° left-handed DNA rotation has been considered 83 , in which the (i+1) th NCP unit was rotated along the entry side linker vector after connecting to its previous NCP. iii) if the (i+1) th NCP was identified to clash with any of the previous NCPs, a new NCP will be randomly selected from the pool and resembled onto the ith NCP by repeating step i through iii. In order to simulate more realistic dynamics of the nucleosome fibers, the unwrapping events were also considered into the system. The origin of the entry and exit DNA linker vectors on the helical NCP track were also varied based on the experimentally measured distribution when constructing the pool of the NCP units.

Insilco-construction of the nucleosome condensates
To reconstitute the chromatin-like higher order structures, a pool of 100 of the above constructed fibers was prepared by using low-and high-salt angle and unwrapping parameters. The assembling of these fibers followed two approaches: "seamless assembling" vs "random docking". In the first approach, the condensates were assembled with no space between array fibers. In this case, the density of the condensate is simply equal to the average density of the above simulated chromatin fibers, which yielded the upper bound to estimate the condensate density. In the second case, all fibers were treated as rigid-bodies and then sequentially fitted into a spherical volume density in a diameter of 200 nm given an inter-array distance. This distance was estimated by the experimentally measured i and i+2 NCPs distance, which mostly reflects the equilibrium NCP interaction distance within an array. To do so, array fibers were converted into lowpass filtered density maps with specific threshold cutoff values that followed the distribution of the measured i and i+2 NCPs distance. The produced fiber density maps were randomly fitted into the spherical volume by using the Chimera sequential fitting function without clashing. The central cubic volume in 20 nm length was cropped out from the generated large condensate for the final representation. The measured density from the second approach was used as the lower-bound to estimate the condensate density.

Classification of the nucleosome array conformations
To identify some of the low energy conformational state of the nucleosome array, which should be observed more often than others based on Boltzmann distribution, nucleosome arrays were classified based on their structural similarity. For each type of the array, 30 to 40 conformations were obtained and subjected to a pair-wise alignment through minimization of the root-mean-square deviation (RMSD) between each pair of models. Based on the values of the constructed distance matrix evaluated by RMSD, the array conformations were sorted and classified with an optimization algorithm that minimized the tree spanning (method option = 'single') using the hierarchical clustering in Scikit-Learn of Python package 84 . The final result was displayed in a dendrogram produced by the scipy.cluster.hierarchy package.

Visualization of the structure dynamics of the nucleosome array
The pseudo dynamics of nucleosome array was simulated by morphing through array conformations in an order followed by the hierarchical clustering. The morphing begins with the structures at the lowest branch of the dendrogram and stops at the highest branch. By using Targeted Molecular Dynamics simulation (TMD) of the NAMD2, the coarse-grained structure was steered from one conformation toward to another under the amber SIRAH force field. The moving forces used in the TMD were pre-calculated based on the distances of the corresponding CG model atoms. After 1,000,000 to 4,000,000 steps (corresponding to a total time of 20 ns to 80 ns, with each step of 20 fs) at 298 K temperature within the implicit solvent, the simulation was terminated when the real-time RMSD fell below 3 Å. Among these morphing structure pairs, ~10% of them were eliminated due to the DNA intertwining with each other. The morphing movie was also displayed with the same order as above from the most populated structure states to the rarest states.

Analysis of in-silico nucleosome array fiber stiffness, hydrodynamic radius, and density
The persistence length was obtained from the orientation correlation of segments along the polymer, expressed as a scalar product of two segment vectors separated by a contour length of L and predicted from the random coil statistics as an exponential decay governed by the persistence 85 . To estimate the persistence length of long simulated nucleosome array fibers containing = 100 NCPs, the NCP pathway of the array was first converted into smoothly connected fixed length line segments ( = 15 ). The persistence length was then estimated from the calculation of the autocorrelation decay along the trajectory of these line segments using the polymer.PersistenceLength.run script from MDAnalysis 86 . The segment length (r) was defined by the distance between two distal ends of the line segment ( = 15 ) and the resulted exponential decay fitting were averaged form 100 simulated array models. The hydrodynamic radius of the fiber was calculated from its persistence using the following equation based on worm-like chain model: where is the radius of gyration of the polymer, P is the persistence length of the simulated fiber, and L is the contour length of the fiber which is proportional to the number of line segments ( = ), which in turn is the number of threaded nucleosomes. The hydrodynamic radius ℎ is related to the radius of gyration by the equation: ℎ = 0.662 . The fiber density was calculated as: /( 4 3 3 ).