domain: sdv nlin phys

text:
Fabrizio Cleri
email: fabrizio.cleri@univ-lille.fr

Stefano Giordano

Ralf Blossey

Nucleosome array deformation in chromatin is sustained by bending, twisting and kinking of linker DNA

Chromatin in the nucleus undergoes mechanical stresses from different sources during the various stages of cell life. Here a trinucleosome array is used as the minimal model to study the mechanical response to applied stress at the molecular level. By using large-scale, all-atom steered-molecular dynamics simulations, we show that the largest part of mechanical stress in compression is accommodated by the DNA linkers joining pairs of nucleosomes, which store the elastic energy accumulated by the applied force. Different mechanical instabilities (Euler bending, Brazier kinking, twist-bending) can deform the DNA canonical structure, as a function of the increasing force load. An important role of the histone tails in assisting the DNA deformation is highlighted. The overall response of the smallest chromatin fragment to compressive stress leaves the nucleosome assembly with a substantial plastic deformation and localised defects, which can have a potential impact on DNA transcription, downstream signaling pathways, the regulation of gene expression, and DNA repair.

Introduction

Cells in tissues are affected by mechanical forces at different length and time scales. Both intracellular and extracellular mechanical signals affect the evolution of cell life, ranging from the internal mechanics of actin-myosin motors at the molecular scale, to the forces driving patterning and organogenesis in embryonic development, to cell-cell and cell-matrix contact interactions in fully developed tissues [1]. Mechanical tension can activate downstream signaling pathways that regulate cell division and growth, orient the transcriptional machinery and drive cell differentiation.

Cells generally respond to mechanical forces by altering their intracellular tension, through a coordinated cytoskeletal rearrangement and actomyosin contraction [START_REF] Gensbittel | Mechanical adaptability of tumor cells in metastasis[END_REF]. Cancer cells in particular are found to be especially sensitive to mechanical stimuli in their microenvironment [START_REF] Gensbittel | Mechanical adaptability of tumor cells in metastasis[END_REF][START_REF] Northcott | Feeling stress: the mechanics of cancer progression and aggression[END_REF], by responding with an increase in the overall stiffness of the tumoral tissue, and the occurrence of a "solid stress" from the growing tumor mass in a confined environment, which is not observed in healthy tissues [START_REF] Nia | Solid stress and elastic energy as measures of tumour mechanopathology[END_REF]. For this reason, tumours typically appear more rigid than most normal tissues, despite individual tumoral cells may display stiffness variations in either direction [START_REF] Guimar Ães | The stiffness of living tissues and its implications for tissue engineering[END_REF]. As a result, cancer cells at the same time produce, and are exposed to important physical forces, composed of complex and largely varying combinations of tension, compression and shearing [START_REF] Belaadi | Under pressure: mechanical stress management in the nucleus[END_REF].

Mechanical tension can affect not only components of the cell surface, but it can also regulate molecular processes within the nucleus, such as gene expression, or even induce DNA damage [START_REF] Shah | Nuclear deformation causes DNA damage by increasing replication stress[END_REF]. The Linker of Nucleoskeleton and Cytoskeleton (LINC) protein complex was found to directly connect the nucleus to the cytoskeleton and hence to the plasma membrane [START_REF] Crisp | Coupling of the nucleus and cytoplasm: role of the LINC complex[END_REF], the nucleus being part of a continuous physical network spanning the extracellular matrix, the cytoskeleton and the nuclear envelope [START_REF] Verstraeten | Experimental techniques for study of chromatin mechanics in intact nuclei and living cells[END_REF]. Hence, environmentally mediated forces can be transmitted to the nucleus and induce deformations of chromatin. While the nuclear structures display visco-elastic mechanical properties, experiments showed that chromatin is rather viscous and has liquid-like behavior [START_REF] Pajerowski | Physical plasticity of the nucleus in stem cell differentiation[END_REF]. In such experiments single cells are aspired in a glass pipette, realising a compression and shear stress that is very similar to that encountered, e.g., by metastatic cancer cells during extravasation and tissue invasion. Notably, fluorescent imaging of histones showed that chromatin is progressively compacted during aspiration; in some nuclei, chromatin bundles aligned and visibly extended within the pipette, indicative of a large physical remodeling that is mainly responsible for nuclear plasticity. Moreover, experiments in which cells are forced to pass through micrometric constrictions [START_REF] Shah | Nuclear deformation causes DNA damage by increasing replication stress[END_REF], showed the formation of strand breaks in the DNA directly correlated to the high mechanical stress, which could also lead to rupture of the nuclear envelope.

With the present study, we shed a first light on the molecular processes of stress transfer and relaxation at the scale of the individual chromatin units, the nucleosomes. By using molecular dynamics simulations of force-induced deformation of an array of three nucleosomes under ideally controlled conditions, we show that external forces acting on such a nucleosome cluster transmit a mechanical stress, which is mainly translated as elastic energy stored in the elastic and plastic response of DNA. (Here, elastic and plastic deformations refer to reversible and irreversible configurational transformations of the molecular structure.) The ability of the double-stranded DNA helix to absorb and release mechanical stress, most notably in the form of bending, twisting and kinking deformations controlled by localised molecular damage, may constitute a platform to elicit or repress the interaction with remodelers, by controlling the access to histone domains. More generally, the concerted action of mechanical deformations and remodeling enzymes may open the way for a new framework, to understand the microscopic control of chromatin organization by mechanical forces, and the attending gene expression as well as the activity of transcription factors and repair proteins.

The mechanical response of DNA in the elastic regime is well understood in terms of the worm-like chain (WLC) model [START_REF] Bustamante | Ten years of tension: single-molecule DNA mechanics[END_REF], which holds until the typical curvature is larger than the persistence length of DNA (∼40-50 nm). However, understanding what happens (well) beyond the elastic regime is interesting both theoretically and in practice, since DNA in the chromatin experiences curvature around and below ∼10 nm (the diameter of the nucleosome). Compression of long polymers [START_REF] Pilyugina | Buckling a semiflexible polymer chain under compression[END_REF] and dsDNA [START_REF] Bleha | Force-displacement relations at compression of dsDNA macromolecules[END_REF] has been recently studied, both by extensions of the WLC and computer simulation. However, such studies mainly focused on the theoretical response of ideal, long monomer chains, which however are scarcely representative of DNA in the chromatin context, since typical linker length is 20-80 base-pairs, and nucleosome-free regions are rarely longer than 150-200 base-pairs [START_REF] Valouev | Determinants of nucleosome organization in primary human cells[END_REF][START_REF] Ranjan | Nucleosome-free region dominates histone acetylation in targeting SWR1 to promoters for H2A.Z replacement[END_REF].

In recent years, nucleosomic DNA deformation has been most often studied in the context of single, isolated nucleosomes, both experimentally and theoretically [START_REF] Zhou | Nucleosome structure and dynamics are coming of age[END_REF]. Attention was focused on the endopening movements or "breathing" [START_REF] Polach | Mechanism of protein access to specific DNA sequences in chromatin: a dynamic equilibrium model for gene regulation[END_REF][START_REF] Polach | Effects of core histone tail domains on the equilibrium constants for dynamic DNA site accessibility in nucleosomes[END_REF], sliding [START_REF] Armeev | Histone dynamics mediate DNA unwrapping and sliding in nucleosomes[END_REF][START_REF] Brandani | DNA sliding in nucleosomes via twist defect propagation revealed by molecular simulations[END_REF], or loop formation [START_REF] Pasi | Extra views on structure and dynamics of DNA loops on nucleosomes studied with molecular simulations[END_REF], in which short DNA segments temporarily detach from the bulk of the nucleosome core particle. In the simpler case of breathing movements, the force providing such a fully reversible DNA deformation could in fact arise from thermal fluctuations. However, other more complex deformation modes require a combination of sliding and twisting along the DNA main axis, intrinsically irreversible and energy consuming. While such deformations might appear spontaneously in some computer simulations, however in the context of chromatin they require external forces, e.g. provided by ATP-powered chromatin remodelers [START_REF] Becker | Nucleosome sliding: facts and fiction[END_REF][START_REF] Bowman | Mechanisms of ATP-dependent nucleosome sliding[END_REF][START_REF] Morgan | Remodeler catalyzed nucleosome repositioning: Influence of structure and stability[END_REF].

2/23

3/23

Given the above considerations, we selected an array of three nucleosomes, the trinucleosome, as the minimal model system to observe the effects of stress and strain relaxation (Figure 1a). Notably, trinucleosomes have been studied experimentally with several techniques, e.g., FRET and SFM [START_REF] Bussiek | Trinucleosome compaction studied by fluorescence energy transfer and scanning force microscopy[END_REF][START_REF] Poirier | Dynamics and function of compact nucleosome arrays[END_REF][START_REF] De | Robustness of trinucleosome compaction to a-tract mediated linker histone orientation[END_REF], cryo-EM [START_REF] Cai | The in situ structures of mono-, di-, and trinucleosomes in human heterochromatin[END_REF] and SAXS [START_REF] Mauney | Solution structure(s) of trinucleosomes from contrast variation saxs[END_REF]; theoretically, only a few coarse-grained molecular dynamics study were reported to date [START_REF] Chang | Histone acetylation dependent energy landscapes in tri-nucleosome revealed by residue-resolved molecular simulations[END_REF][START_REF] Kenzaki | Linker DNA length is a key to tri-nucleosome folding[END_REF]. As starting points, we use one complete experimental trinucleosome structure from cryo-electron microscopy (in the following called T169, see Methods section below), and one structure assembled by repeating three identical units of an experimental mononucleosome (called T183, see Methods). In all cases, we simulate the application of a constant force to the center of mass of the histone octamer in the central nucleosome (C, in the figure), attached by its two linkers to two other nucleosomes (A1 and A2, identical for simplicity), which provide representative mechanical reactions of the background chromatin structure, to displacements and internal deformations induced by external forces. In the following it will appear that the largest part of the deformation under the various initial configurations is taken up by the DNA linkers, with a very active participation of the histone tails. However, such a complex configuration including the whole structure of three consecutivxe nucleosomes is necessary to reproduce the correct loading conditions. A simpler configuration in which, e.g., just a fragment of straight dsDNA is compressed [START_REF] Pilyugina | Buckling a semiflexible polymer chain under compression[END_REF][START_REF] Bleha | Force-displacement relations at compression of dsDNA macromolecules[END_REF], would not be enough to capture the whole range of possible loading variants, such as spontaneous torsion and bending of the whole nucleosome about its superhelical axis, which variously transmit and redistribute to the entire structure the loading force. Furthermore, a relevant part of the modelling will be dedicated to the subsequent relaxation phase, in which the force and all restraints are removed, and the trinucleosome is left free: in this stage it will be important to study if, and to what extent, the deformed structures can somewhat recover a stress-free configuration.

RESULTS

DNA linkers respond to compressive stress by bending, twisting and kinking

Compressive stress applied to the central nucleosome builds up a large amount of elastic energy, which appears to be mostly stored in the two DNA linkers. While C moves in the negative-z direction under the effect of the applied force (Fig. 1c), the two linkers become increasingly deformed, with respect to the initial C-B distance L 0 . Fig. 1b displays a summary of several steered-MD compression simulations at constant applied force, on the two trinucleosome systems T169 (thick lines, full vs. dashed) and T183 (thin lines). Force constants between about 40 and 400 pN are shown. The thick full or dashed lines refer to the T169 with two different initial arrangements of the H3 histone tails (see below, Section 2.3). The force threshold at which deformation actually starts can be estimated by looking at the force value for which the ratio L/L 0 starts to decrease from the (nearly) constant value ≃ 1. Such a threshold appears to be comprised between 40-85 pN for the T183, and between 85-125 pN for the T169 (note that we are always compressing two DNA linkers that equally share the load, at least in the initial stages of the deformation). In general, all curves appear to converge at long times to a deformation of about L ≃ L 0 /2, half of the initial C-B distance (the low-force simulations require much longer times). However, as it will be shown in the next Sections, the type of deformations observed can be rather different, among the different conditions. Notably, there seems to be no unique, final equilibrium configuration, despite all simulations starting from the same initial atomic structure. For the T169, Suppl.Mat. correspond, respectively, to the two dashed/full black lines at f =200 pN in Fig. 1b (obtained with different initial conditions, see below), the (g) notably showing the C nucleosome to roll about the A nucleosomes. Panel (i) corresponds to the full purple line at f =400 pN in Fig. 1b, showing that at large force the very rapid displacement just pushes the C nucleosome straight along z; in this case, both DNA linkers appear to undergo supercoiling, a rather exceptional occurrence given their short length. A summary of the set of MD runs discussed in this work is provided in Suppl.Mat. Table 1.

Notably, the overall deformation is almost exclusively supported by the DNA linkers, while the body of the moving nucleosome experiences but a slightly adjustement of its configuration (and even less the two anchoring nucleosomes) during the compression. This can be readily appreciated by looking at the RMS displacements and fluctuations during the whole compression time. For the representative case of T183 at 125 pN compression, Suppl.Mat. Fig. S3 shows in panel (a) the RMS displacement with respect to the initial configuration, as a function of time: it is seen that the two DNA linkers experience a very large excursion, whereas the DNA wrapped around the histones, and the histones themselves, quickly settle to a constant, low value, meaning very little deformation with respect to the initial state. Panels (b-f) display the time-averaged RMS fluctuation for each residue: it can be appreciated how all values fluctuate between about 0.1-0.2 nm, with the notable exception of the histone-tail regions that can wildly fluctuate with amplitudes up to 1.-1.5 nm (see also the striking difference between the dashed and full red lines in panel (a), for the histones with or without tails); this will be the subject of an accurate analysis in the following.

From a continuum-mechanics perspective, the uniaxial compression of a long cylindrical rod may lead to a sequence of mechanical instabilities, represented schematically in Figure 2a. The first bending instability (Euler) is encountered at small initial load; for a homogeneous and isotropic rod, the limiting load P c = π 2 EI/L 2 to initiate buckling is (i) proportional to the product EI of the Young's modulus E of the material, and the second moment of the cross section area, I = πR 4 /4 for a circular radius R, and (ii) inversely proportional to the square of the length L 2 of the rod. If we take E=300 MPa and R=1 nm, the critical load is P c ≃25 pN for a dsDNA of 10 nm length. This may be taken as an ideal lower bound for the force necessary to induce bending. At the same time, DNA clearly cannot be viewed as a dense homogeneous and isotropic cylinder. A more appropriate comparison could be with a hollow thin-walled 5/23 cylinder, filled with a softer material [START_REF] Luongo | The Brazier effect for elastic pipe beams with foam cores[END_REF]; also, the intrinsic twist and water embedding should modify the limiting load P c . This behavior indeed corresponds to the elastic regime of the DNA response, described by the WLC model.

As the applied load increases the approximately circular DNA cross section may turn elliptic (or squeezed) at some point. This induces a second instability (Brazier), in which the rod develops a sharp kink. The original analysis of Brazier for an infinite 2D homogeneous strip was reformulated for a 3D rod of finite length [START_REF] Takano | Interaction between Euler buckling and Brazier instability[END_REF], also expressing the result in terms of the critical Euler bending moment. However, for a heterogeneous polymer such as DNA, whose structure includes an intrinsic helical twist, predicting the Brazier instability is very difficult. Notably, the event of kinking may have a number of different possible origins, when molecular conformations of the stacked base-pairs are considered, as we will show in the following. This regime is beyond what can be described by the WLC model and is characterized by plastic (i.e., permanent, irreversible) deformation.

Moreover, a bifurcation toward a twist-bend deformation mode can alternatively develop (see Fig. 2a). For twist to occur as a result of simple bending, the structure must have different bending stiffness along (at least) two directions not parallel to its main axis. If, moreover, the longitudinal fiber also carries a helical symmetry axis, upon bending the spiral line will be elongated or shortened, according to whether the added twist is in the same direction of the unstressed helix, or opposite to it [START_REF] Hartog | Advanced strength of materials[END_REF]. This extra tension is mostly longitudinal, but has a component in the transverse plane and is directed along the tangent, thus adding up to a torque. Twist-bend coupling has been invoked in DNA under different conditions, more commonly under stretching and twisting (for analogy to single-molecule experiments), or spontaneous looping [START_REF] Marko | Stretching must twist DNA[END_REF][START_REF] Strick | Twisting and stretching single DNA molecules[END_REF][START_REF] Kulic | Twist-bend instability for toroidal DNA condensates[END_REF][START_REF] Nomidis | Twist-bend coupling and the torsional response of double-stranded DNA[END_REF]. In the present case, it will peculiarly result as a response to the excess elastic energy from the compression force.

By applying load to the DNA linkers, it is possible that the above instabilities develop to a different degree and at different points. Figure 2b shows the result of mechanical instabilities in the two (left/right) 6/23 [START_REF] Strick | Twisting and stretching single DNA molecules[END_REF]-bp linkers from a same MD simulation of the T183 trinucleosome: both linkers start with a Euler-like buckling at low load; however, upon increasing deformation one (green) develops a kink at a well defined site, near the middle of the double helix; while the other (blue) bends and over-twists at three different sites, one central and two close to the ends (which overlap with the entry site of the two nucleosomes A and C). The sketched configurations span the first 10 ns, from a 20-ns long simulation.

The development of Brazier-kinking can be readily appreciated by looking at the helical DNA parameters, "opening" and "buckle" [START_REF] Olson | A standard reference frame for the description of nucleic acid base-pair geometry[END_REF], defining the deformation between two bases in each base-pair (bp), as computed by the CURVES+ module [START_REF] Strick | Conformational analysis of nucleic acids revisited: Curves+[END_REF] (see Methods section). The time diagrams in Figure 3a,b describe the evolution of the deformed state in the same DNA linker shown in cyan color in Fig. 2b (or "left linker"). The vertical ordinates span the discrete pathlength s i along the DNA (see Suppl.Mat. Appendix I). It can be appreciated how the narrowing of the cross section of the DNA pseudo-cylindrical tube is associated with the in-plane distortion of the base-pairs 22-23, where the two bases loose their hydrogenbond coordination (upper panels in the Figure ; in the following Section we will analyse in more detail such distortion). The in-plane deformation is also correlated with the "tilt" (θ ) and "roll" (ρ) parameters (lower panels in the Figure ), which measure the orientational distortion between two consecutive base-pairs, 22 and 23 in this case.

The twist-bending instability, on the other hand, arises when a sharp increase in local curvature is accommodated by opening the double helix, rather than squeezing it into a kink. The local curvature about each bp (assuming the DNA as a tube with different bending stiffness parallel to b and n in the cross section, see Suppl.Mat. Appendix I) can be obtained from local "tilt" and "roll" measures, as κ ∝ (θ 2 + ρ 2 ) 1/2 . The time plots of Figure 4, relative to the other DNA linker (the blue one in Fig. 2b, or "right"), show the evident correlation between the continuous increase in the local curvature (upper left panel) and the under-twisting (upper right panel) of the double helix, centered about bp 7-8, and bp 22-23. Note that the normal straight DNA has an average twist of about 35 degrees between consecutive base pairs: the two purple bands show that the local curvature correlates with the under-twisting, which decreases to below 10 degrees in correspondence of bp 21; the two lower panels in the Figure also show that such a deformation tends to slightly open up the helix at the same site, by adding a component of stretching and staggering between the two facing bases (that is, the two bases tend to move away and slide along z from each other). Note that the peaks in stretch/stagger also roughly coincide with the blackish peaks of large under-twisting. An additional, broad anti-correlation can be deduced by looking at the opposite behavior of twist vs. roll parameters, shown in the Suppl.Mat. Fig. S4 (top row): the over-twisting (yellow) is correlated with a decrease in roll (purple), while under-twisting correlates with increase in roll; at the same time, the tilt parameter fluctuates between much smaller positive/negative values. As a result, the curvature is mostly dictated by the "softer" roll component, while the tilt component is "stiffer" and less subject to large variations. Both effects appear in broad agreement with the notion of twist-waves suggested by Carlon et al. [START_REF] Skoruppa | Bend-induced twist waves and the structure of nucleosomal DNA[END_REF], a typical signature of twist-bend coupling in DNA.

Base-pair flipping at the origin of Brazier instability

The mechanical instabilities observed in the previous Section have a clear origin in the molecular rearrangements following the distribution of the applied load. We show in this Section that the initiating event is often a relative sliding of the bases in a given bp, or between adjacent bp, leading to a stretched and staggered configuration. In this case, the two bases find themselves in a position which, via a subtle balance of interatomic forces and thermal fluctuations, may run into a kinematic bifurcation leading to different outcomes.

A comparative example is shown in Figure 5, from a compression MD simulation of the T183 system. One of the linkers, as shown in Fig. 5a from two neighboring, T-A (blue-green) and G-C (orange-red), base-pairs. By contrast, the other linker displays a shift between the two bases within a single base-pair. While the first results in sharp kinking, the second softens via twist-bending.

The left Fig. 5a describes the first case. In the leftmost MD frame, corresponding to about the first ns of simulation, the G-C pair (orange-red) splits off with a large staggering movement, while the T-A pair (blue-green) opens up in the bonding plane and the T (blue) starts to turn away from the helical axis. About 4 ns later (right MD frame) the guanine from the G-C pair has moved on top of the adenine from the other bp, while the thymine from the T-A pair is seen to completely flip out in extrahelical position. Such configuration is stabilized by the formation of a number of "exotic" hydrogen bonds, as show in the time-plot at the bottom of the figure: the two thin dashed lines reporting the H-bond lengths of the T-A pair soon jump to exceedingly large values (8-12 Å), thus signaling the opening of the bp with the thymine flipping out; at about the same time, two H-bonds form between the adenine and the cytosine of the two adjacent bp, one notably involving the O3' oxygen from the adenine ribose, and the other an adenine nitrogen; and a third H-bond forms within the G-C pair, but strangely involving a phosphate O2P oxygen. The result is a sharp kink at the site of the two adjacent bp, which is 'macroscopically' interpreted as a kind of Brazier instability.

The molecular rearrangement following the base flipping involves also a change in the solvation structure of the DNA double helix. We analysed the distribution of water molecules and ions around the G-C/T-A dinucleotide, in comparison with a normally stacked pair. All the water molecules whose oxygen atom lied within 3.25 Åfrom any other atom in the dinucleotide, and all ions within 6 Å, were used to build the distribution; the Suppl.Mat. Fig. S5 compares the H-bonded water network for the normal (a) and flipped-out (b) configuration of the dinucleotide. The overall number of water molecules comprised within the cutoff is nearly the same, as well as the number of H 2 O directly bound to the nucleobase atoms. However, it is observed that the dense network of H-bonds between waters and the oxygens in the backbone phosphates, is substantially reduced in the flipped-out configuration, compared to the normal 8/23 one (21 vs. 28 molecules); because of the extrahelical rotation of the T95, part of these backbone-bound water molecules tend to "invade" and further stabilize the kinked region of the DNA double helix. Also, H 2 O chains of 2, 3 and in some case 4 molecules are observed to link the stacked bases, across both the minor and major groove in the normal dinucleotide, whereas only a few shorter (single-molecule) links are surviving in the flipped-out configuration. Furthermore, in the normal configuration, Na + cations are observed to nest in the minor groove, while in the flipped-out configuration such ions are dispersed around the T-A pair, coordinating with both bound and unbound water molecules.

The right panel, Fig. 5b, describes the staggering and opening of the T-A bp 22 (see also Fig. 4) in the opposite linker. The T and A bases (blue and red), seen in their normal bonding configuration in the left MD frame, switch to a staggered position in which one lies above, but distant from the other. In this case the relative distance of the two bases, as monitored by the time variation of the 2+1 H-bond lengths in the bottom plot, experiences a temporary excursion coinciding with the time of twisting of the double helix about that site. However, the H-bond lengths never extend beyond ∼6 Å, and can turn back to nearly regular values; notably, it can be seen that while the two strong H-bonds are overstretched, the third, usually very weak bond between adenine O2 oxygen and the thymine C2 hydrogen, becomes shortened and helps maintaining a (later) recoverable conformation. This overall movement under-twists the DNA helix, which yields to the bending by adding a twist component to the deformation; while the near normal H-bonding configuration is eventually recovered, the twisting partly remains as embedded plastic deformation of the double helix.

While our statistics is of course limited, we always observed such events in each compression simulation, for both the T169 and the T183. Typically, 2-4 twist-bending events per simulation were observed, and often 1 kinking. Notably, once the extrahelical flip-out is realised in the kink, it is extremely resilient, and rarely observed to heal back during MD relaxation simulations in the range of up to 200 ns.

Such observations can contribute to the discussion about DNA softening at large deformation. Al-

9/23
ready F. Crick, back in 1975 [START_REF] Crick | Kinky helix[END_REF], had suggested kinks as a possible source of deformation in DNA (imagining DNA wrapped around histones in small straight segments joined by kinks). Later, by using a simple discrete-energy model for short dsDNA fragments, Ramstein and Lavery had brought to attention that bp opening is facilitated by bending and, conversely, once a bp is disrupted DNA can bend more easily [START_REF] Ramstein | Energetic coupling between DNA bending and base pair opening[END_REF]; later on, it was also suggested that denaturation "bubbles" of 3-5 bp, by forming short ssDNA strands, could rather be the origin of sharp kinks ("teardrop" shape) [START_REF] Yan | Localized single-stranded bubble mechanism for cyclization of short double helix DNA[END_REF]. Experimental data however were not conclusive, and appear to fit both interpretations [START_REF] Cloutier | DNA twisting flexibility and the formation of sharply looped pro-teinDNA complexes[END_REF]. Lankaš et al. [START_REF] Lankaš | Kinking occurs during molecular dynamics simulations of small DNA minicircles[END_REF] observed kinking in DNA minicircles, as a kind of extreme unstacking between consecutive bp. Our results represent a first, completely atomistic test of DNA bending at extremely large deformation under realistic embedding in chromatin, and show that very localized extrahelical base flipping can induce sharp kinking, accompanied as in Ref. [START_REF] Lankaš | Kinking occurs during molecular dynamics simulations of small DNA minicircles[END_REF] by a large negative roll. At the other end, local bp disruption in the form of intra-bp staggering may allow large under-twisting coupled to positive roll, which also leads to short regions of tight curvature with a radius of just a few nm. It seems therefore that melting bubbles should not be needed to explain the extreme DNA flexibility under large deformation, which instead seems to be mediated by very localized base-pair disruption events of various kind.

DNA linker length, symmetry, and the role of histone tails

By comparing the response of the T169 and T183 under the same loading conditions, we can observe possible differences of the mechanical response to the applied load, according to the geometrical characteristics of the trinucleosome system. In the two systems, each pair of DNA linkers have exactly the same bp sequence, for a length of respectively 22 and 36 bp, that is a fraction of about 15% and 30% of the DNA persistence length λ P . Clearly, as the linker length approaches λ P DNA becomes a flexible string, with small mechanical resistance to compressive stress. Genomic analyses reveal typical spacings of 14-28 bp in non-transcribed regions, of which the T169 may be representative, while the T183 could be at the lower limit of active-region spacings [START_REF] Valouev | Determinants of nucleosome organization in primary human cells[END_REF] (in which, however, also much larger inter-nucleosome distances can occur, up to 100-200 bp, even longer than λ P ).

In practice, only a minor length dependence is observed between the qualitative response of the T169 and T183. The force thresholds from Fig. 1b are slightly different, the longer T183 starting to yield at a somewhat lower force, in the interval 40-80 pN, vs. 80-120 pN for the T169 (despite the estimation being subject to a considerable fluctuation). However, in both systems the DNA linkers experience similar mechanical instabilities upon perpendicular loading. Within the limited statistics, no clear sequence dependence of the base-pair disruption was recorded. What may be interesting, conversely, is the fact that even a DNA length as short as 22 bp can display the same flexibility of longer fragments, if under severe mechanical constraints.

Notably, the linkers in the two systems have by construction a difference in symmetry. We must remember that in bending a "tube" of finite cross-section, the half tube distant from the centre of curvature is in tension, while the half tube closer to the centre is in compression (the center line of the tube being the neutral axis). Now, the 36-bp linkers in the T183 are symmetric with respect to the loading (z-direction), whereas the 22-bp linkers in the T169 are almost perfectly mirror-reflected about the z-axis. This means that in the T169 one of the linkers has the major grooves in tension and the minor grooves in compression, but the grooves in the other linker are inverted. In the T183, instead, both linkers are always in the same tension/compression condition. The two movies in the Supplementary Material (see discussion below) clearly display the consequence of such symmetries.

Another important difference, which however affects in the same way linkers of either length, is observed with the interaction between DNA and the histone tails. The role of the histone tails in nucleosome dynamics has been investigated in detail recently in 15 µs simulations ( [START_REF] Armeev | Histone dynamics mediate DNA unwrapping and sliding in nucleosomes[END_REF]; for a review of earlier work, 10/23 see Ref. [START_REF] Huertas | Breaths, twists, and turns of atomistic nucleosomes[END_REF]). In particular, the H3 histone tails protrude from the two sides of each nucleosome and extend toward the linkers. However, histone tails are extremely flexible and disordered regions, with a length of 37 amino acids for the H3, 28 a.a. for H4, 20 a.a. for H2B, and 15 a.a. for the H2A N-terminals, this latter comprising also a mobile tail of 10 a.a. at the C-terminal; that is, theoretical contour lengths comparable or even longer than each DNA linker. The large mobility of the tails is clearly visible in panels (b-f) of Suppl.Mat. Fig. S3, where the tail regions are shaded in light blue and light red; also, panel (a) shows that all the RMSD of the histones is contributed only by the tails, see dashed vs. full red lines.

The key role of histone tails in initiating (and/or constraining) the mechanical destabilization of the DNA linker, is supported by a series of simulations in which we removed the tails of all histone copies in the three nucleosomes. (Tail-less nucleosomes are often used in experiments, although this is a purely artificial in-vitro construct.) During 20-ns MD compression runs of the T169 at forces of 50 and 100 pN, we could not observe any typical signs of the mechanical instabilities inducing the defective DNA structures observed above. In all cases, the top (C) nucleosome moves against the two lower ones (A1,A2) by tilting and turning about its center of mass, such that the DNA linkers are flattened and broadly curved, with a rather sharp bending at each of the entry sites, and accommodate moderate elastic deformations along their respective helical axis. The Suppl.Mat. Fig. S6 shows the final configuration at t=20 ns of the green/blue linkers, for the 50 pN (a) and the 100 pN (b) MD simulation; in the lower panel (c), a few of the key DNA helical parameters for the green linker are shown (the ones for the blue being qualitatively similar), it can be appreciated that no precursor signs of instabilities appear (compare the scale values with those of Fig. 4 and5). Upon free relaxation, all compressed structures are observed to retrieve defect-free DNA conformations.

By comparison, the choice of the relative position of the H3 tail and the DNA linker at the beginning of the compression, affects dramatically the subsequent dynamics of the stress transfer. Figure 6 shows in the left panel the details of the H3-tail-linker interaction in the T169, for two different initial conditions (Fig. 6a, blue vs. red). The time evolution is shown in Fig. 6b-d, at times t=1,5,10 ns: the blue tail conformation, initially adherent to the DNA linker, evolves smoothly and "helps" the linker to bend in a smoothly curved loop; the red one, initially away from the DNA, moves abruptly toward the middle of the linker, and "pushes" it, like an intruding finger, into a drastic curvature of about κ=0.25 nm -1 . The key event leading to this extreme bending is detailed in the panels of Fig. 6e,f, where the G-C pair contacted by the H3 tail tip undergoes a "forced flip-out", resembling the initiation of the Brazier-like instability: the pair is disrupted, and the guanine-80 is forced into extrahelical configuration, as also demonstrated by the sudden breaking of two of its G-C hydrogen bonds.

The short length of the linker in the T169 also shows another interesting cooperative effect, when the H3 tails from opposite nucleosomes interact with the same DNA linker. Figure 7 depicts such a situation, by showing a DNA linker sandwiched between a red and a blue H3 tail, respectively from the upper and lower nucleosome. For the snapshots at subsequent times, we calculated the contact surface area (in nm 2 ) and the free energy of adhesion ∆G (in kcal/mol), with the PDBePISA utility [START_REF] Krissinel | Inference of macromolecular assemblies from crystalline state[END_REF]. Both tails increase steadily their contact with the linker, the blue one changing more evidently by increasing the contact surface by a factor of 3, and more than doubling the ∆G of adhesion. The effect is that of forming a sort of "cradle", that holds the DNA fragment and accommodates its bending. Due to the close proximity of the two tails, there is little room for a smooth curvature in this case, and the final effect is that of a kinking, similar to the Brazier instability but without the help of the extrahelical flipping of a nucleobase.

Due to the linker asymmetry, the relaxation to the initial state after the compression phase is very different for the T183 and the T169. As shown in Suppl.Movie 1, in the T183 the two symmetrically arranged linkers bend "on the same side", and are finally able to relax back to a configuration close to the initial one. In particular, the two (cyan vs. blue) linkers end up in a fully extended geometry, like the initial one, albeit the overall trinucleosome conformation is not identical, the upper nucleosome ending up with its superhelical axis tilted by about 60 degrees with respect to the initial configuration. In the movie it can also be observed the behavior of the thymine-95 (shown in red), which flips out of the DNA main axis some time after the start of the compression, thus originating the kink for the Brazier-like instability. This structural defect does not disappear completely, since the nucleobase flips back into the correct position at the start of the relaxation phase, but subsequently continues to move in and out to partly-flipped positions.

The T169 instead shows a dramatic effect of the linker asymmetry upon relaxation, in combination with the role of the H3 tails, which contribute to the asymmetric response. The Suppl.Movie 2 shows the parallel simulations of fast compression (10 ns) of the T169, followed by slow relaxation (200 ns), with two, very different initial configurations of the H3 tails. In both simulations, the two mirror-arranged linkers bend toward opposing directions, because of the mirror symmetry of their arrangement. In the left panel, at the beginning of the simulation the H3 tails are symmetrically contacting the upper end of the cyan linker (one red tail from the C nucleosome), and the lower end of the blue linker (yellow tail from the A nucleosome); the other two tails are starting just close to the other two DNA ends. At the end of the relaxation stage, the left simulation has the cyan linker getting back to the fully extended geometry, while the blue linker regains only in part its straight conformation. This is due to a minor extent to the interaction of the H3 tails with the DNA, but more importantly to the lock formed between the C and A nucleosomes by the other H3 tail of A, and the H4 tail of C, which maintain the three nucleosomes in a closely compact cluster. In the right panel, instead, all the four H3 tails (two from the upper C nucleosome, and one each from the bottom A nucleosomes) start from an open, non-contact configuration. After about half compression, the blue linker is smoothly curved (κ=0.12 nm -1 ), while the cyan linker has been intruded by the tip of the (red) H3 histone tail. Such a deformation is in fact irreversible, and perturbs the subsequent relaxation. Whereas the blue linker slowly regains the straight configuration, with the opposite (red and yellow) H3 tails interacting with its ends, the cyan DNA linker is instead tightly supercoiled around the (red) H3 tail (as also shown in Fig. 6d) and shows a Brazier-like kink at the site of maximum supercoiling.

DNA deformation locally redistributes the elastic energy

A key question in molecular simulations is how to correlate the observed deformations with the molecular forces which are at the origin of such deformations. The notion of microscopic stress arises as necessary, to recapitulate mechanical information contained in large-scale MD trajectories of nonuniform systems, as well as to connect molecular details with continuum physics at larger scales. The stress tensor is the appropriate quantity to describe how a set of external forces is redistributed, according to the reaction of the internal structure of the material. The proper definition of mechanical stress is that of momentum of a force flowing across an infinitesimal surface, for any arbitrary orientation of its boundary. As such, the molecular stress must include both a kinetic and a potential contribution, [START_REF] Tsai | The virial theorem and stress calculation in molecular dynamics[END_REF][START_REF] Giordano | Brittle fracture: from elasticity theory to atomistic simulations[END_REF] and can in principle be obtained from atomic positions and velocities extracted from MD trajectories. We already demonstrated the insight that can be obtained with mechanical stress calculations in a previous work, dealing with isolated mono-nucleosomes [START_REF] Cleri | Mechanical evolution of DNA double-strand breaks in the nucleosome[END_REF] (see Suppl.Mat. Appendix III for details).

Figure 8 displays some components and integrals of the stress tensor for a simulation of the T183 system at 125 pN constant force. All plots are shown as continuous lines obtained by smoothing the 13/23 noisy discrete atomic values (see Methods). Also, all data are plotted as a function of the pathlength 1 < s < 36 along the DNA central axis, since the z coordinates of the base-pair centers are continuously changing in the deformed configurations.

The plots in panel (a) show the σ zz component of the stress (parallel to the direction of the applied force, between the center of mass of the C and B groups) at time t=1 ns, which in these very early stages of the simulation practically corresponds with the total applied stress. The green-L linker is all under a compressive load, with a positive value of σ zz everywhere, except the last fragment s >27, where it fluctu- ates around zero, indicating that the transmission of the compressive force from the C (top) nucleosome is not yet complete along the DNA structure, at this time. The blue-R linker is also generally under compression, however it shows several spots with a slightly tensile condition (σ zz < 0); this should indicate a role of the under-twist (see again Fig. 4, upper-right) that starts interfering with the pure bending deformation, around s=7-10, 21-23, 28-29, already at such early stages. At later stages of increasing deformation, the simple σ zz component will no longer be an appropriate quantity to characterize the deformation; instead, the deviatoric (Von Mises) stress, which combines differences among the diagonal components and the squared-average of the off-diagonal components (see eq.( 10) in Suppl.Mat. Appendix III), will provide a more synthetic information about the state of deformation.

In panel (b) we report the elastic energy density ε(r) (see eq.(!1) in Suppl.Mat. Appendix III), for

14/23

both linkers at the end of compression stage (t=100 ns). This is a measure of the mechanical elastic energy accumulated in the system, upon the action of the external force. At the beginning, both linkers start storing some elastic energy, for a total of 4.8 (blue) and 6.4 MPa (green) respectively. (Note that such totals are obtained by summing the discrete values for each atom, not by a numerical integration of the continuous line.) Upon increasing load the green-R linker develops the kink, which allows to concentrate the deformation energy around s=17-18, while the two parts of DNA before and after the kink tend to remain rather close to straight, with just a shallow curvature. The blue-R linker instead, displays two maxima close to the two ends, where the linker starts to join with the two nucleosome core particles, and a broad central maximum, the whole plot lying however below the green one. This indicates that twist-bend coupling response is characterized by a lower accumulation of elastic energy, compared to the kinked green linker at the same MD time. By looking at panel (a) of Suppl.Mat. Fig. S3, it can be seen that the kink deformation (cyan line) develops in the early stages of compression, after which the two DNA sections above and below the kink settle back to lower RMSD, whereas the twist-bending (blue line) continues to accumulate during the whole compression: overall, the twist-bend coupling seems to do a better job at redistributing the elastic energy, instead of concentrating at one highly deformed site.

On the other hand, when looking at Fig. 8 panels (c) and (d), where are respectively reported the average curvature, and the deviatoric stress density ε dev (r), we get a different picture, complementary to that of panel (b). According to both measures, the blue linker demonstrates a larger degree of deformation, and an increased local curvature (averages taken over groups of 6 consecutive base-pairs), compared to the green linker. The green DNA shows a shallow curvature of about κ ∼0.05 nm -1 (that is, average radius of 20 nm), except the two small humps immediately close to the kink site at s=17, where the curvature is slightly higher (κ ∼0.08 nm -1 ). Curvature in the blue linker, instead, has three peaks with κ ∼0.15 nm -1 (bending radius ∼6-7 nm), broadly matching the deviatoric stress maxima in panel (d).

The overall picture that these stress measures are providing, is that the (blue linker) twist-bend coupling brings a large structural deformation and extreme local bending, with curvature radii as compact as 6-7 nm (not far from the curvature of DNA wrapping around the nucleosome core); the formation of angled kinks (green linker) instead maintains a rather straight, less deformed and less stressed DNA. However, the sharply kinked structure can locally concentrate a much higher density of elastic energy, compared to the case of a continuous twist-bend, which better redistributes the elastic energy along the whole length.

DISCUSSION

In summary, in this work we focused on the structural effects of mechanical compression on the smallest significative chromatin unit, the trinucleosome. External forces were applied to a central nucleosome sandwiched between two "anchoring", adjacent nucleosomes. Very-large-scale, all-atoms molecular dynamics simulations of the whole system under external forces, embedded in a molecular water bath with physiological ion concentrations, gave a comprehensive picture of the force redistribution in the deformed samples, leading to the following, main conclusions:

1) The elastic/plastic deformation is concentrated in the DNA linkers, while the nucleosome core particles remain practically compact and unaffected, in the explored range of compressive forces (40-400 pN). Monitoring the helical DNA coordinates provides precursor signatures of the sites and modes of maximum deformation.

2) Even very short DNA linkers (7.5 nm, much shorter than the DNA persistence length) can experience extreme elastic and plastic deformations, when severely constrained.

3) Longer DNA linkers display a smaller force threshold to initiate deformation, starting from a 15/23 straight DNA conformation; in the limit of a length approaching the DNA persistence length (λ P ∼ 50 nm), even thermal fluctuations could induce deformations. 4) Kinks can form along the linker, either by a narrowing of the cross-section mediated by a molecular defect (base-flipping), or by a strong interaction with external agents, such as histone tails. These originate a kind of Brazier instability, with near-straight segments of DNA split by a constriction. 5) Large amount of under-twisting can be induced at different sites along the DNA contour length; twist-bend coupling can result in supercoiling of the linker, which seems to be partly irreversible.

6) Structural details (e.g. symmetry of the linker arrangement, initial position of the histone tails) determine the plastic response and the "final state" configuration. It is possible that different sequences of A-T vs G-C pairs could affect the local deformability; this is left for further study.

7) Upon application of external forces, the internal reactions redistribute the deformation in the molecular structure. Calculation of the stress field components provides an essential link between the observed deformation, and the underlying atomic forces.

All such observations are conducive to a deeper understanding of the force transmission, from the outer parts of the cell to the inside of the nucleus, during movements entailing strong compression and shearing of the cell, e.g., migration, nucleokinesis, diapedesis, tissue invasion by cancer cells and so on, with potentially important consequences in downstream signaling pathways and regulation of gene expression.

A final word may be in order concerning the question of the DNA-histone adhesion. One early motivation of our simulations was to investigate the so-called "breathing" movements of DNA about the histone core. Differently from the case of isolated nucleosomes, which often made the subject of experimental and simulation studies, in the continuous chromatin fiber any localized nucleosome deformations in which DNA would detach from the core particle must be subject to mechanical constraints. In our simulations, carried out with pulling forces quite larger than the forces required to unwrap DNA from nucleosomes (as measured by optical tweezers [START_REF] Bennink | Unfolding individual nucleosomes by stretching single chromatin fibers with optical tweezers[END_REF][START_REF] Pope | Single chromatin fiber stretching reveals physically distinct populations of disassembly events[END_REF], of the order of 5-25 pN), we could never observe DNA unwrapping events, apart from minor fluctuations of a few Å amplitude for the terminal base-pairs of each nucleosome (see panel (b) in Suppl.Mat. Fig. S3, where the terminal bp of the nucleosomewrapped DNA display but a minor fluctuation in excess of the average). Blossey and Schiessel [START_REF] Blossey | The dynamics of the nucleosome: thermal effects, external forces and ATP[END_REF] offered well-grounded comparisons between the thermal forces and the forces needed to unwrap, twist and slide DNA about the nucleosome protein octamer, showing that all such processes are energetically too costly, to occur spontaneously at any practical rate. To this, we can add that such considerations become even more stringent when regarded in the context of nucleosomes linked together in chromatin strands, as in the present simulations. For example, an eventual "breathing" mode must apply a tensile stress on the DNA linker between two nucleosomes, and the corresponding force required to drag around a whole nucleosome against the hydrodynamic resistance of the nucleoplasm turns out to be too large, to be simply provided by thermal fluctuations. By taking Stokes' drag on a spherical particle of size ∼ 10 nm, moving in a fluid of viscosity ∼ 2 mPa•s (typical value for the nucleoplasm), the ratio F/v is of the order of 2 • 10 -7 pN•s/nm. If we imagine one nucleosome moving opposite to the other at ∼1 nm/µs, because of the DNA breathing, the mechanical force applied on the symmetric linker terminals would be ∼0.1 pN, barely sufficient for displacing a DNA end.

MATERIALS AND METHODS

Molecular dynamics simulations

We performed all-atom molecular dynamics (MD) simulations of the compressive deformation of two trinucleosome systems, in the following defined as T183 and T169, with linker lengths of 36 and 22 base-16/23 pairs (bp), respectively. By using steered-MD, an external force was applied to the central nucleosome, which is pushed against the other two. As it is seen, in all cases the largest amount of mechanical deformation is shared and taken up by the DNA linkers, notably with an important role of the histone tails.

In each of the simulated systems, the 147 bp of DNA of each nucleosome follow the 601-core sequence [START_REF] Lowary | New DNA sequence rules for high affinity binding to histone octamer and sequence-directed nucleosome positioning[END_REF], while the linkers have a generic sequence. The T183 was reconstructed by assembling three copies of the 197 nucleosome (PDB file 5NL0, [START_REF] Bednar | Structure and dynamics of a 197 bp nucleosome in complex with linker histone h1[END_REF]), dropping the H1 histone and manually adjusting the linker lengths in CHIMERA [START_REF] Pettersen | UCSF-chimera, a visualization system for exploratory research and analysis[END_REF]; (it is worth noting that the experimental sequence is not exactly a pristine Widom-601, but a modified, palindromic-601).The T169 was instead directly taken fro the PDB entry 6L4A [START_REF] Takizawa | Cryo-EM structures of centromeric tri-nucleosomes containing a central CENP-A nucleosome[END_REF]. To maintain the representation of a random piece of chromatin, we avoided using the linker lengths multiple of 147+10n, which are known to force the zig-zag conformation of multiple nucleosomes because of the rotational constraints imposed from exact DNA turns. The DNA sequences of the two trinucleosomes are detailed in the Suppl.Mat. Table 2. In both systems, the protein octamers were replaced by copies of the canonical human H3/H4 and H2A/H2B complete with full histone tails, from PDB file 1KX5 [START_REF] Davey | Solvent mediated interactions in the structure of the nucleosome core particle at 1.9 A resolution[END_REF]. Since this hybrid experimental structure is built with human DNA wrapped around Xenopus proteins, we manually substituted the human protein sequence in CHIMERA ("swapaa" command). A comparison of the human vs. frog sequences is given in Suppl.Mat. Table 3, showing that only very minor differences exists between the two. Notably, the largest sequence difference in the more mobile histone tails is observed for the H2A histone. However, the H2A dimer is arranged parallel to the dyad axis, the N-terminal tails remaining always close to the superhelical positions +4/4, and the C-tails close to the dyad. Given the negligible deformations observed in the nucleosome core, such differences are expected to have a minor dynamical effect, if any. Where needed, the primary structure of the histones was completed for the missing aminoacids with the Swiss-Model utility of SwissProt. [START_REF] Duvaud | the Swiss Bioinformatics Resource Portal, as designed by its users[END_REF] All molecular dynamics (MD) simulations were run with the GROMACS-2020 package. [START_REF] Berendsen | GROMACS: A message-passing parallel molecular dynamics implementation[END_REF][START_REF] Lindahl | GROMACS 3.0: a package for molecular simulation and trajectory analysis[END_REF] Proteins were described with the AMBER14 force field database, [START_REF] Lindahl | 14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB[END_REF] including the PARMBSC1 extension for nucleic acids. [START_REF] Pérez | Frontiers in molecular dynamics simulations of DNA[END_REF] For the preliminary thermal stability {NVT} and {NPT} simulations, the ensembles of the complete trinucleosomes were solvated in water boxes of size 20OE20OE35 nm 3 for the T169, and 32OE32OE40 nm 3 for the T183, with periodic boundary conditions in the three directions, respectively containing about 430,000 and 1,225,000 TIP3P water molecules (about 1.4 and 3.8 million atoms), plus enough Na + and Cl -ions to ensure neutralization of the phosphate backbone charge, at a concentration of 0.1 M NaCl. Similar conditions were also used for the umbrella sampling and force-driven simulations. All the production MD runs were carried out at the temperature of 310K and pressure of 1 atm.

Coulomb forces were summed with particle-mesh Ewald sum (PME), using a real-space cutoff of 1.2 nm, equal to the cut-off radius of shifted Van der Waals potentials. We used rigid bonds for the water molecules with a time step of 2 fs for the thermal equilibration phases, and 1 fs for production and forcepulling runs. For the thermal equilibration, preparatory {NPT} runs with temperatures increasing in steps of 100K from T=10 to T=310K lasted 20 ns each; followed by thermal stability {NVT} simulations which extended to 100-200 ns for each configuration.

For the steered molecular dynamics simulations we identified four pull-groups in either the T169 and T183 systems (see schematic in Fig. 1). The central nucleosome (C, red) was forced against the two anchoring nucleosomes (A1 orange, A2 yellow). Holding the A1, A2 nucleosomes at their average position is assumed to mimic the effect of the embedding in the large background chromatin mass. Either constant-force or an umbrella sampling potential were applied to displace the center of mass of C along the negative z direction, towards the center of mass of a subset of atoms (B, blue) approximately covering the central region shared between A1 and A2. We used values of the umbrella spring constant between 5-10 kJ/mol/nm 2 , and displacement rates of 10 -4 -10 -6 nm/ps. Constant-force simulations were run 17/23 with forces ranging between 40-400 pN. In all simulations, the A1 and A2 nucleosomes were held still about their average position by applying a fictitious constant force with zero rate between their respective centers of mass, along the z direction. Further positional restraints with harmonic springs along z were applied to the heavy atoms in the center part of each nucleosomes A1, A2, about 2,400 heavy atoms within a radius of 2,5 nm from each nucleosome center. According to the displacement rate, typical stress cycles required ∼10-100 ns to compress the nucleosome C towards the baseline of A1-A2, followed (in a few simulations) by 150-200 ns of free motion, during which C could relax back, to get somewhat close to its initial configuration; in relaxation, all atoms are free to move and any restraints are removed.

Typical runs on the SGI-8600 JEAN-ZAY supercomputer of the IDRIS French National Center in Orsay used clusters of 400 to 800 Cascade Lake-6248 2.5 GHz cores, and 40 to 120 NVIDIA 100 SXM2 GPU cards, according to the system size, for a total of about 1M CPU hours. A summary of the production runs is given in the Suppl.Mat. Table 1. Some PDBs of the main atomic configurations listed in the table are freely available from the web repository FigShare under the doi 10.6084/m9.figshare.23634396.

Mechanical stress calculation

Continuous stress fields were extracted from MD trajectories by the so-called covariant central-force decomposition scheme (CCFD) for the intra-and intermolecular forces, which ensures conservation of linear and angular momentum of the molecular systems under very general conditions. [START_REF] Torres-Sanchez | Examining the mechanical equilibrium of microscopic stresses in molecular simulations[END_REF] The method is implemented in a special-purpose patch to the old GROMACS 4.6 version, which reads (all or part of) a MD trajectory for the selected subset of atoms for which stress is to be computed, and performs the entire analysis. Since the patch (called GROMACS-LS) constrains the code to run in serial rather than in parallel, care must be taken to define properly the atomic subset of interest, in order to avoid prohibitive computing times. Comparing stress fields from different MD runs requires an extra care, since the structures need to share exactly the same box size and center, to avoid numerical artefacts; when averaging a sequence of frames from a simulation it is important to subtract any translation or rotation with respect to the first frame. To ensure compatibility between the different versions, we firstly extract the required frames for the subsystem considered (typically a few thousand atoms) from the binary trajectory file into a text pdb file (GROMACS-2000 utility TRJCONV, two successive runs, firstly with option "pbc nojump", and the second with option "fit progressive"). Then, a pseudo-trajectory file is rebuilt from the pdb with the TRJCONV from GROMACS 4.6 (this is necessary to produce a readable binary trajectory for GROMACS-LS). According to the CCFD scheme, stress fields are calculated by GROMACS-LS on a continuous grid of 'voxels' superposed on the molecular structure. In all calculations we used cubic voxels of 1x1x1 Å 3 . Stress components and individual force contributions (pair, angle, dihedral, etc.) can be projected back on the atom sites by defining a conventional (but non unique) volume around each atom, including one or more voxels; we used to average over blocks of m 3 Å 3 with a variable m=5 to 10, including the few nearest neighbors of each given atom. Each stress calculation is averaged over an "instant" of 50 ps of the trajectory, including 10 frames spaced at 5 ps.

series of midpoints along s, projected in the best-fitting common plane. Note that the direction vector of R does not necessarily coincide with b, which instead varies from point to point: this is the very reason to introduce a "global", rather than "local" notion of curvature. The numerical procedure is described in the following Appendix II. In practice, the geometric curvature κ = 1/R is estimated by computing R for a series of equal lengths S n = ∑ k+n i=k s i , the starting point k spanning the whole DNA length S 0 ; each set overlaps with the next one by ±1 midpoint, therefore the values of κ(s) behave quite smoothly for a sufficiently large n. In the calculations shown in Suppl.Fig. S1 we used n=5-8, such that each set spans a length approximately equal to one half-turn of the double helix. However, by using larger values for averaging, large bending movements can be captured. Suppl.Fig. S1c shows an example of curvature calculated with n=25 points.

Appendix II. Numerical approximation of global curvature

We describe the method to fit a circle to a cluster of points in 3D space, used to estimate the global curvature along the DNA backbone. Consider a set of n points {P 1 , .., P n }, where P i = (x i , y i , z i ) T ∈ R 3 , obtained by extracting the centers of the pairs of P atoms of each bp along the DNA backbone. For a DNA segment containing N base-pairs, we consider subsets of n < N consecutive points, typically n ≃ 6 -12, covering the whole DNA length, possibly with some overlap between adjacent subsets. We want to find a circle that fits as close as possible to each subset of points. The circle fitting method can be split into three steps:

1. use the least-squares method to find the best fitting plane to the set of points;

2. project all the points perpendicularly onto the fitting plane in 2D; 3. use again the least-square method, to fit a circle in the 2D plane to the set of projected points, and obtain the circle center and radius.

Eventually, the circle center can be transformed back to 3D coordinates, if we want to collect the ensemble of curvatures in the common reference frame other than just the scalar value of κ. We will use two different implementations of the least-squares method to obtain the radius of the best fitting circle.

Step1. Given the subset of n points in 3D, the best plane can be found by a least-squares method that minimizes the distance of each point to the plane. The equation for a generic plane is: ax + by + c = z. Hence, we build an eigenvalue problem like Ax = B, or:

    x 0 y 0 1 x 1 y 1 1 .. .. .. x n y n 1       a b c   =     z 0 z 1 .. z n     (3)
Then, solve for the x vector of coefficients. However, since each subset is made up by more than just 3 points, the system is over-determined. Therefore, we use the left pseudo-inverse matrix: A + = (A T A) -1 A T . Finally, the coefficients of the plane are found as:

  a b c   = (A T A) -1 A T B ( 4
)

25/37

Note that such a method is not entirely general. A more robust method should be to subtract out the centroid of the subset of points, to form a n × 3 matrix with the resulting coordinates, and calculate its singular value decomposition. As a result, the normal vector of the best-fitting plane should be found as the left singular vector corresponding to the least singular value. However, we preferred the method of least squares because of its simpler implementation in this case.

Step2. Once the best fitting plane is obtained, the 3D coordinates of the subset of points are projected onto this plane according to their respective perpendicular distances, thus obtaining a 2D representation of the same subset in the new coordinates.

Step3. The implicit equation in 2D for a circle with radius R and center (x c , y c ) T can be arranged as:

(x -x c ) 2 + (y -y c ) 2 = R 2 , 2x c x + 2y c y + (R 2 -x 2 c -y 2 c ) = x 2 + y 2 , ( 5
) c 1 x + c 2 y + c 0 = x 2 + y 2
with c = (c 1 , c 2 , c 0 ) T the vector of unknown parameters. Then, by applying this definition to all the input points P i , it yields to a system of linear equations:

Ac = b (6)
with:

A =   x 1 y 1 1 .. .. .. x n y n 1   , b =   x 2 1 + y 2 1 .. x 2 n + y 2 n   , (7)
Since there are more equations than unknowns, an approximate solution is obtained by the method of least-squares, which minimizes the squared sum of the residuals ||b -Ac|| 2 ; we implemented for this purpose the family of subroutines QR_SOLVE from LINPACK. Therefore, the center of the fitting circle is x c = c 1 /2, y c = c 2 /2, and the radius is R = (c 0 + x 2 c + y 2 c ) 1/2 . Then, the line curvature for the subset {P 1 , .., P n } is just κ = 1/R.

Appendix III. Stress calculation for atomistic systems

In classical macroscopic continuum elasticity, stress is constructed as a continuous field at each point r in a homogeneous domain where a distribution of forces exists, namely a 3x3 tensor, σ αβ (r), with the dimensions of a force per unit surface. In practice, for a force vector f and a boundary S with normal unit vector n, the Cartesian components of force and stress are linked as:

f α = ∑ β S (σ αβ • n β )dS, α, β = x, y, z (8)
In atomistic simulations, however, we deal with material points exerting forces across empty space. Hence, the definition of the analogous of stress at the atomic level is complicated by several issues, such as: discretization of a continuum field; what is the volume around each atom; how to precisely define a boundary surface; last but not least, the fact that atoms have a velocity [START_REF] Cleri | Representation of mechanical loads in molecular dynamics simulations[END_REF][START_REF] Shi | A perspective on the microscopic pressure (stress) tensor: History, current understanding, and future challenges[END_REF]. With a much simplified notation for the sake of clarity, Hardy's definition of a continuous stress that can be mapped onto atomic 26/37 positions, velocities and forces (that is, accelerations) of a system of N point particles ("atoms") can be written as [START_REF] Hardy | Formulas for determining local properties in molecular dynamics simulations: shock waves[END_REF]:

σ αβ (r) = - N ∑ i< j=2 (r α i j • f β i j )b i j (r) + 1 2 N ∑ i=1 m i (v α i • v β i )g i j (r) (9)
where r i j = |r i -r j | is the scalar distance between atoms i, j, and f i j is the force on atom i from any other atom j; the b i j and g i j are spatial weighting functions (usually derived by variational arguments). The first term in the equation represents the potential (or "virial") contribution, and the second term is the kinetic contribution to the molecular stress.

In this work we use the recently developed covariant central-force decomposition scheme (CCFD) for the intra-and intermolecular forces [START_REF] Torres-Sanchez | Examining the mechanical equilibrium of microscopic stresses in molecular simulations[END_REF], which is based on thermodynamical arguments rather than on momentum conservation. We aim to extract such information from our simulations and couple it to the state of deformation, along the lines of our first application of this method to a single nucleosome [START_REF] Cleri | Mechanical evolution of DNA double-strand breaks in the nucleosome[END_REF].

In the condition of uniaxial loading used in this work, the component of the stress parallel to the loading condition, σ zz (r), may be a first, interesting quantity, at least at the beginning of the simulation; in the early stage of loading, the compressive force (directed along the line joining the centers of mass of the groups C and B, in Fig. 1a) is parallel to the z-direction. Therefore, the σ zz (r) component is practically coincident with the σ 3 principal component of the diagonal stress.

At later times, however, the configuration quickly becomes so much deformed that other components of the stress tensor describe the redistribution of the internal forces in response to the compression. Next to the purely compressive component, also mixed transverse and shearing components arise. A quantity that can be useful to characterise the complex deformations occurring, instead of looking at the different components one by one, is the distribution of deviatoric, or Von Mises, stress:

ε dev (r) = 1 + ν 3E 1 2 ∑ i< j
(σ ii -σ j j ) 2 + 3 ∑ i< j σ 2 i j [START_REF] Pajerowski | Physical plasticity of the nucleus in stem cell differentiation[END_REF] where all the σ i j = σ i j (r), and E=300 MPa and ν=0.4 are typical values for the Young's modulus and Poisson's ratio for straight DNA, which fix the energy scale. Note that by definition this quantity is always positive, although the various stress components may have both negative and positive values. At a more global scale, another important quantity that can be derived from the local stress tensor is the density of elastic energy at each point in space:

ε el (r) = 1 E ∑ i, j
](σ ii ) 2 -2ν(σ ii σ j j ) + 1 + ν 2 (σ i̸ = j + σ j̸ =i ) 2 [START_REF] Bustamante | Ten years of tension: single-molecule DNA mechanics[END_REF] By looking at the time evolution of this elastic energy at different positions in space, it can be observed how the elastic energy, progressively stored by the constant loading force, is redistributed. 1. Summary of the steered-MD simulations presented in this work. The upper and lower blocks correspond to the T183 and T169, respectively. For the latter, the '*' indicates configurations produced with alternate H3 tails initial configurations. The name of the pdb file containing the final configuration of the run, and available from the repository, is indicated whenever available. Note that water and ions are not included in the pdb.

Figure 1 .

 1 Figure 1. (a) Schematic of the trinucleosome. The center of mass of the central nucleosome C is pulled by a constant force toward the center of mass of the reference group B; the supporting nucleosomes A1 and A2 are held fixed about their average positions, by applying zero-rate force to the respective centers of mass. The ideal line C-B defines the z-axis, positive toward the top of the page. (b) Displacement plots at constant force for the C→B nucleosome compression. Thick full/dashed lines for the T169 (see text), thin lines for the T183 system. The relative C-B distance is normalised to the initial distance L 0 . Constant force values (pN) indicated by colors: 40/green, 85/blue, 125/red, 200/black, 400/purple. (c) Snapshots from the T183 simulation at 125 pN (see thin-red line in (b)).

23 Figure 2 .

 232 Figure 2. (a) Mechanical instabilities in a axially compressed rod. (b) Progressive deformation of two 36-bp DNA linkers into a Brazier-kink (cyan, left-to-right), and a twist-bend mode (blue, right-to-left); both linkers extracted from the same T183 simulation at 125 pN.

Figure 3 .

 3 Figure 3. Time-plots for the buckle (top left) and opening (top right) between each two bases in a bp, and tilt (bottom left) and roll (bottom right) inter-bp helical DNA parameters, computed along the L (cyan) DNA linker of Fig.1c. The vertical axis indicates the position 1-34 of the base pair; the horizontal axis is simulation time. Note that the color scales next to each plot have different units and ranges.

  (left panel), displays a shifting in the relative position of two bases 7/23

Figure 4 .

 4 Figure 4. Time-plots for the local curvature (top left) and twist helical parameter (top right), and for the stretch (bottom left) and stagger (bottom right) intra-bp helical parameters, for the R (blue) DNA linker of Fig.2b. See Fig3 for nomenclature.

Figure 5 .

 5 Figure 5. MD compression simulation of the T183 system. (a) Sequence of molecular events leading to a Brazier kinking instability via extrahelical base flipping. (b) Sequence of molecular events leading to twist-bending. In each panel, the left figure represents the initial state and the right figure a snapshot of the evolved conformation. The time-plots below each panel depict the evolution of the hydrogen bonds in the respective base-pairs. The red and blue vertical bars indicate the approximate time to which the two snapshots correspond.

Figure 6 .

 6 Figure 6. (a-d) Snapshots at subsequent times of two MD simulations of the T169 system, starting with the same DNA configuration and different initial arrangements (blue vs. red) of the H3 histone tail from the "upper" (C in Fig.1) nucleosome; initial configurations superimposed in (a) at t=0; frames (b,c,d) at t=1,5,10 ns (only the first 20 amino-acids of the H3 tail are depicted). (e) Conformation of the G-C base pair indicated in the grey rectangle of (d, red): above, beginning of simulation; below, end of simulation. (f) Time-plot of the hydrogen bonds in the G-C pairs: thick line C82.O2-G82.H2, thin line C82.N3-G82.H1, thin-dashed line C82.H4-G82O6.

Figure 7 .

 7 Figure 7. Sequence of snapshots from a MD simulation of the T169 system. It is shown one DNA linker (cyan) "cradled" between two H3 histone tails, respectively coming from the adjacent nucleosomes, the upper, moving C (red tail) and the lower, fixed A (blue tail). The figures below each snapshot give the contact surface in Å 2 and the ∆G of adhesion of the contact surfaces, in kcal/mol; line above (blue) for the blue H3 tail, line below (red) for the red H3 tail.

Figure 8 .

 8 Figure 8. Sample stress plots for the T183 trinucleosome system. Blue/cyan lines refer to the R or L linker (see also Fig.2). In all panels the abscissa is the position along the contour length s of the DNA linker. (a) Component σ zz of the stress parallel to the compression direction, at the beginning of the simulation (t=1 ns). (b) Elastic energy density (see Suppl.Mat. Appendix III, eq.(11) at time t=100 ns (final stage of the compression). (c) Average line curvature at time t=100 ns (see Suppl.Mat. Appendix II for definitions). (d) Deviatoric stress (see Suppl.Mat. Appendix III, eq.(10) at time t=100 ns.

Table 2 .

 2 DNA sequences of the two trinucleosome systems employed in this work; the original Widom-601 sequence is also reported for comparison. Red letters indicate the DNA linker regions; blue letters indicate the hanging DNA leads; boxed letters indicate the dyad.

12/23

Acknowledgements

We acknowledge useful discussions with E. Carlon and his group in Leuven. The authors acknowledge funding from the French ANR Project ANR-21-CE45-0032 "Dyprosome". We thank generous computing time allocation on the JEAN-ZAY supercomputer of IDRIS CNRS in Orsay, under project GENCI A0130712986.

Author contributions statement F.C. and R.B. conceived the computer experiments. F.C. conducted the molecular dynamics simulations. All authors analysed the results and reviewed the manuscript.

Additional information

Atomic configurations in PDB format from the molecular dynamics simulations in this work are available at the repository FigShare under doi 10.6084/m9.figshare.23634396. The authors declare no competing interests.

Supplementary Materials

Appendix I. DNA reference frames and conventional coordinates Analysis of stress fields and molecular deformation was carried out with our in-home NA-BAKE code, which allows to extract the principal components of the stress, compute invariants and other derived quantities, compare stress fields from different simulations, and write the outputs in the portable Gaussiancube format for visualization. The code also includes the set of subroutines CURVES+, [START_REF] Strick | Conformational analysis of nucleic acids revisited: Curves+[END_REF] for a detailed analysis of the DNA helical parameters and coordinates. These latter are defined according to the socalled Tsukuba conventional reference frame (Suppl.Fig. S1a), with the three orthonormal basis vectors {e 1 , e 2 , e 3 } such that e 3 is the tangent unit vector along the center line defining the DNA helical symmetry axis, and e 1 is the unit vector pointing at the major groove.

It is however worth noting that the base-pair reference frame defined in Ref. [START_REF] Olson | A standard reference frame for the description of nucleic acid base-pair geometry[END_REF] and used by CURVES+ can become extremely distorted, given the large deformations between bases observed in the present work. For this reason, NA-BAKE also uses a second, "smoother" Cartesian reference frame, centered at the midpoint of each pair of P atoms in a base-pair (Suppl.Fig. S1b), the phosphate backbone being more rigid than the stacked base-pairs. At each midpoint, the three basis vectors are defined as follows: the tangent vector τ lies on the direction joining two midpoints in adjacent base-pairs; the normal vector n points at the center of the major groove; and the binormal vector is b = n × τ. (This can be seen as the approximation of a Frenet-Serret continuous frame for the broken polyline connecting the midpoints of the base-pairs.) The origin of this reference frame does not coincide with the {e 1 , e 2 , e 3 }, which depends on the type of base-pair, whereas in the P-centered frame this is uniquely defined as the mid-point along the P-P ideal segment; the origin is also shifted along τ||e 3 by about half a repeat distance ("rise") with respect to the Tsukuba base-pair plane.

The ensemble of midpoints ordered according to the bp numbering, defines a pathlength 0 ≤ s ≤ S 0 along the DNA contour length. For a strand of N bp, with midpoints defined by the set of position vectors {r 1 , ..r N }, s is a sum of discrete segments (a polyline):

The pathlength vectors s i are parallel to the local tangent τ. The variation of the three vectors n(s), b(s), τ(s) along the DNA pathlength can be used to identify the local bending and torsion of the line. From the discrete variation we compute three Euler angles, defining the 3D rotation of the local reference frame between two bp i and i + 1; in NA-BAKE this is done by using quaternion algebra.

Because of the double-helix geometry there are two bending modes with different elastic constants A and B (stiffness, usually given in units of nm -1 ), i.e. the DNA polyline has an effective thickness, and is anisotropic. By looking at the definitions of the helical parameters, it turns out that the variation of n(s) corresponds to the "tilt" θ , that of b(s) corresponds to the "roll" ρ, and that of τ(s) corresponds to the "twist". The first two parameters also give the local curvature between consecutive bp along the line, as

However, such a definition of "strictly local" curvature may miss the ample bending deformations that extend over lengths longer than just 1-2 bp (see Suppl.Fig. S1c). For this reason, we introduced a more geometric-minded, "global" notion of curvature, by calculating the radius R of the best fitting circle to a

36/37

Supplementary Movie 1 -Slow compression (100 ns) of the T183 trinucleosome, followed by free relaxation (150 ns), at a force of 125 pN. The two linkers are represented in cyan and blue, as in the main text. The red nucleotide on the cyan linker is the thymine-95, which flips out in extrahelical position some time after the beginning of compression, and induces the kinking (Brazier-like instability). Water molecules and ions not shown. Supplementary Movie 2 -Comparison of the rapid compression (10 ns) of the T169 trinucleosome, followed by free relaxation (200 ns), at a force of 125 pN. The left and right panels correspond to two identical simulations, starting with different initial configurations of the H3 histone tails. Only the H3 histone tails implicated in the contacts are shown as colored surfaces (red, orange and yellow for the C,A1,A2 nucleosome) plus the H4 tail (grey) from C nucleosome, while the rest of the proteins is shadowed; water molecules and ions are not shown.

37/37

Applied Compression

DNA linker supercoiling compressive stress

DNA curvature

Graphical Abstract (for review)

Click here to access/download;Graphical Abstract (for review);graph_abstr.pdf -Understanding force transmission to the molecular cons4tuents of chroma4n during cell and nucleus compression, e.g. in cell migra4on, 4ssue invasion by cancer cells, and so on, has important consequences in downstream signaling pathways, regula4on of gene expression, DNA repair.

-During simulated compression and shearing of a small chroma4n element, the trinucleosome, molecular deforma4on is concentrated in the DNA linkers, while the nucleosome core par4cles remain prac4cally compact and unaffected.

-Kinks, under-twis4ng, extreme bending and supercoiling of the DNA linkers are observed, leaving behind partly irreversible defects aDer compression is released.