Early Critical Phase Transitions of Pericentromere-Associated Domains in MCF-7 Breast Cancer Cells Committed to Differentiation by Heregulin

Finding out how cells with the same genome change fates in differentiation commitment is a challenge of biology. We used MCF-7 breast cancer cells treated with the ErbB2 ligand heregulin (HRG), which induces differentiation, to address if and how the constitutive pericentromere-associated domains (PADs) may be involved in this process. PAD-specific repressive heterochromatin (H3K9me3) and active euchromatin (H3K4me3) marking, centromere (CENPA) labelling, qPCR, acridine-orange-DNA structural test, and microscopic image analysis were applied. We found a two-step DNA unfolding, at 15-20 min and 60 min after HRG treatment, coinciding with bi-phasic activation of the early response genes (c-FOS family) and two steps of critical phase transition which were revealed in transcriptome studies. In control, the distribution of PAD number and size displays a power-law scaling with a boundary at the nucleolus. PADs’ clustering correlates with centromere numbers. 15 min after HRG treatment, the unravelling of PADs occurs, coinciding with the first step of euchromatin unfolding. The second step is associated with transcription of long-non-coding-RNA from satellite III DNA. We hypothesize that splitting of the PAD clusters under the critical size threshold of the silencing domain abrupts position effect variegation. It allows the first genome transcription avalanche to occur, starting differentiation commitment.


Introduction
In human cells and tissues, the genomes represented by the DNA in the cell nuclei are identical. However, differentiated cells are functioning in their own highly specialized way and this means that the genes being available in total must show a different expression depending on the cell/tissue functional requirements. In order to find answers to this challenge, biophysics can offer holistic approaches looking at the genome of a cellular system in its entirety. Apparently, it requires a non-linear thermodynamic approach [1]. The well-established model of luminal breast cancer MCF-7 cells [2][3][4] was used for differentiation induction by the treatment with type I neuregulin 1 (also named heregulin (HRG)), a ligand of HER2-positive breast cancer [5]. HRG induces cooperative trans-activation of ErbB2-receptors through their dimerization (ErbB2/ErbB3), internalization, and nuclear transport of epidermal growth factor receptor (EGFR) [6][7][8][9]. Similarly, neuregulin 1 causes the differentiation of embryonic stem cells into cardiac myocytes [10]. In MCF-7 breast cancer cells, 10 min after HRG application, it stimulates the down-stream signalling of Akt and ERK kinase pathway, coupled with bi-phasic induction of the immediate response genes, the earliest of c-FOS peaking at 30 min [11]. The sustainability of the ERK network activation lasting for 1-1.5 h is dependent on the sequential induction by c-FOS of FOSL-1 (FRA-1) and then FHL2. Interaction of the last two causes a negative feed-back suppression of c-FOS with gradual downregulation of ERK whose activation is abated by 3 h. This bi-phasic circuit of the early gene response was shown indispensable for the irreversible determination of cell fate change by HRG, with the specification of differentiation (production of lipid droplets) displayed several days later [11][12][13]. Analysis of the kinetics of the transcriptome data supplied by this group (22277 transcripts) evaluating the change of the transcription rate allowed Tsuchiya and colleagues to reveal the early critical acceleration of transcription between 15-20 min involving the low-expressed genes and operating with the ~1 Mbp of the topology-associated domains (TAD)-sized units [14,15]. This can be interpreted as a phase transition of the whole genome superstructure. The turbulence of transcription ended after 3 h. In turn, on the same model, Chen and colleagues [16] revealed a critical coordinated upregulation of the promoting group of 104 genes at 1-1.5 h. These studies highlighted two points of the self-organized criticality (SOC) in the applied MCF-7-HRG model at 15-20 min and 1-1.5 h after treatment.
The earlier ideas of the self-organisation ordering of the genome as a whole [17][18][19] were more recently developed by nucleome studies using the chromatin high-structure capture and image analysis. This and similar approach disclosed cooperative assembly and phase separation of the transcriptionally active and inactive chromatin [20][21][22][23][24][25] largely acting on the base of simple physical principles [26]. Studies of cell differentiation, in general, showed an increase of the proportion and perinuclear assembly of the repressed facultative heterochromatin; these studies were set on mature differentiated cells for a few days after receiving the corresponding stimuli [25]. The role of the constitutive heterochromatin in the pericentromere-associated domains (PADs) composed of the 234 bp satellite repeats [27] has been also investigated in similar sets. It was shown that PADs being rather dynamic in differentiation mostly cluster and co-segregate with repressive genomic regions and determine the repression of juxta-posed genes by position variegation effect [27][28][29]. PADs also react by increase of the long non-coding RNA (LncRNA) synthesis to stress conditions [30]. However, the behaviour of PADs has not been explored at the early stage of differentiation commitment described above in transcription studies. Here, we set to address the topological behaviour of PADs in the MCF-7 cells grown, stimulated, and fixed in situ on the chamber slides for the analysis at the early terms after treatment with HRG, when critical phase transitions of the genome transcription have been revealed.

Cell line and treatments
Cells of the human breast cancer cell line MCF-7 (ATCC) were cultured in Dulbecco's modified Eagle's media (DMEM) supplemented with 10% fetal bovine serum (FBS, Sigma-Aldrich, St. Louis, MO, USA). Cells were grown on chamber slides without antibiotics in 5% CO2 incubator at 37°C. For experimental studies, cells were serum-starved for 48 h and treated with 50 nM heregulin (HRG, 396-HB, R&D Systems Inc., Minneapolis, MN, USA or ab50227, Abcam, Cambridge, UK) for different times (10 min -3 h). For ErbB-2 receptor clustering and internalization experiments proving the activation of ErbB2 receptors, after 48 h serum starvation, additionally reduced serum starvation (0,5% FBS) for 3 h was performed. Then the medium was replaced by serum-free medium containing 50nM HRG for 60 min and finally changed for 10% FBS DMEM for another 1 h before fixation.
RNAse A (DNAse, protease-free, 10 mg/ml, Fermentas, Vilnius, Lithuania) in PBS for 30 min at 37°C and rinsed in PBS. Then hydrolysis for 40 s or 5 min at 60°C in 1N HCl was performed. After hydrolysis, the slides were washed in three changes of distilled water (1 min each), passed through PBS for 2 min and 0.1 M acetic buffer pH 4.1 for 5 min. Staining with the high purity Acridine Orange (AO) (Polyscience Inc.) at a concentration of 10 -4 M was performed in the same buffer for 15 min (in the dark) following three changes of AO diluted to a concentration of 10 -6 M in the same acetic buffer, sealed in the latter solution under a glass coverslip with nail polish and immediately examined (AO starts fading after two hours).

Immunofluorescence (IF)
Cells on chamber slides were washed with warm PBS and fixed in 4% paraformaldehyde/PBS for 15 min at RT. Slides were then washed thrice in PBS 0.1 M glycine for 5 min, once in PBS 0.2% Triton X-100 for 5 min and once in PBS for 5 min. Slides were subsequently blocked for 15 min in TBS, 0.025% Tween 20, 1% BSA at room temperature. Specimens were covered with a blocking solution containing the primary antibody and incubated in a humidified chamber at 4°C overnight. Specimens were then washed thrice in TBS 0,01% Tween 20 and covered with appropriate secondary antibodies (Goat anti-mouse IgG Alexa Fluor 488, Goat anti-Rabbit-IgG Alexa Fluor 594, Invitrogen, Carlsbad, CA, USA) and incubated for 40 min at room temperature in the dark. Slides were washed thrice for 5 min with TBS and then counterstained with 0.25 μg/ml DAPI for 2 min, and finally embedded in Prolong Gold (Invitrogen, Carlsbad, CA, USA).
For further illustration, super-resolution single-molecule localization microscopy was used described in detail in [9,31]. The instrument was equipped with a 100x / NA 1.46 oil immersion objective lens. A solid-state laser (491 nm) was used for inducing molecular blinking and 2,000 image frames were detected by an iXon Ultra Andor EMCCD grayscale camera.
AO stained slides were analysed and images were acquired using an epifluorescence microscope Leica DM LS2 equipped with a Sony DXC-S500 colour video camera and using a Leica 40x / NA 0.65 objective. The fluorescence filter cube Leica I3 was used (excitation filter BP420-490 nm, dichromatic mirror 510 nm, suppression filter LP 515 nm).
Image analysis and quantification of the red/green index of AO staining was performed using the Image-Pro Plus 4.2 software (Media Cybernetics, Rockville, MD, USA).

Real-time PCR (qPCR)
The total RNA was extracted from 10 6 MCF-7 cells by using TRIZOL (Invitrogen). First-strand cDNA was synthesized using 2 µ g of RNA, random hexamers and RevertAidTM M-MuLV Reverse Transcriptase (Fermentas, Lithuania) according to the manufacturer's protocols, and subsequently diluted with nuclease-free water ten times. qPCR was run on a ViiA7 (Applied Biosystems). Amplification mixtures (25 µ l) for all amplicons contained 2 µ l template cDNA, 2x SYBR Green Master Mix buffer (12.5 µ l) (Thermo Fisher Scientific, Vilnius, Lithuania) and 100 nM forward and reverse primer.
The cycling conditions for the c-FOS, FOSL-1, FHL2, c-MYC, GABDH and B2M comprised 10 min polymerase activation at 95 o C and 40 cycles at 95 o C for 15 s and 60 o C for 60 s. A melt curve was also performed after the assay to check for the specificity of the reaction. This consists of 20 s at 72 o C followed by a ramp-up of 1 o C steps with 5 s hold at each step. Every cDNA sample was normalized against two housekeeping genes: GAPDH and B2M using giNorm software [32]. The calculated gene expression stability coefficient M was applied to qPCR results.
The cycling conditions for the GAPDH-2, HS3-9 and HS3-1 comprised 10 min polymerase activation at 95 o C and 40 cycles at 95 o C for 15 s, 57 o C for 20 s and 60 o C for 50 s. A melt curve was also performed after the assay to check for the specificity of the reaction using the same conditions as above. Expression was normalized to the GAPDH using the formulas: ΔCT = CT (HS3-9 or HS3-1) -CT(GAPDH), ΔΔCT = ΔCT -ΔCT, RQ = 2 -ΔΔCT.

Statistical analysis
Comparison between groups was performed by two-tailed t-test analysis. For activation and visualization of ErbB2 or ErbB2/ErbB3 dimerization, respectively, the procedures and the antibody were used as described in [9]. After treatment with 50 nM HRG tested at 60 min, we obtained clustering of ErbB2 receptors, their transition from the cell surface and aggregation in the cytoplasm in the cells (Fig. 1A,B). We next checked by real-time PCR (qPCR) the activation and dynamics of the transcription of the key genes of the early gene response at 20 min, 60 min and 3 h. As presented in Fig. 2A, we could reproduce the bi-phasic kinetics reported by Saeki and colleagues [11] showing a very high up-regulation of c-FOS by 20 min, with the following peak of FOSL-1 at 60 min and both downregulation at 3 h, with the opposing dynamics of the third gene of this group, the FHL2 gene. c-MYC, a master inducer of bivalent genes, was reaching an activation peak at 60 min after 50 nM HRG treatment ( Fig. 2A). So, we have validated the MCF-7-HRG model acting in our hands as comparable with the studies carried out on it before. As in these earlier studies, there were demonstrated critical phase transitions of the transcribed genome -at 15-20 min by Tchuchiya et al. [14,15] and at 60-75 min by Chen et al [16] -we decided to probe the DNA conformation at both terms using the acridine orange in situ test.

3.2.
A two-step DNA unfolding induced by HRG is revealed by the DNA structural test with acridine orange (AO), at 20 min and at 60 min.
AO is a metachromatic fluorophore which binds to nucleic acids in two forms -intercalating monomers (fluorescence absorption maximum at 502 nm and emission maximum at 530 nm) and co-parallel dimers attached to the phosphate residues of the DNA backbone (fluorescence absorption maximum at 460 nm and emission maximum at 640 nm). So, after RNA extraction, the conformation changes correlated to DNA double/single-strand transitions (and the interplay between a secondary and tertiary structure can be detected by changes of fluorescence spectra of the dye [35]). This test has been used for chromatin conformation evaluation in many cases [36][37][38][39], particularly for evaluating sperm chromatin in male infertility [40,41]. The red/green index (~intensity ration 640/530 nm) was suggested as a physical parameter of the DNA conformation [40].
The structural DNA test by AO, applied after RNA removal by RNAse, at equilibrium conditions and pH 4.1 is preceded by short provocation with acid pretreatment. It releases more loosely bound histones from DNA phosphate groups, thus allowing to stain them by forming the AO dimers -redder, differently from the chromatin with more stable DNA-histone interactions remaining more green due to preferential intercalation mechanism of the staining in the latter [38,39]. For fine-probing of the DNA structure, the duration of the acid pretreatment can be varied; also air-drying of cells before fixation (elevating the concentration of small counter-ions) can reveal the differences of the various chromatin domains if compared with the undried specimens [37]. We applied this methodology and tested the red/green AO index by imaging the freshly AO-stained cells by random scanning of 250-400 cells with the same colour camera exposure for each time point. As can be seen in Fig. 3A, under these conditions AO was well staining the whole nuclear chromatin, i.e., both, euchromatin and heterochromatin The images were composed of two components, red and green. By removal of the green component of the fluorescence images in Image-Pro Plus, the red one has become visible ( Fig. 3B) and vice versa. The images were digitally separated and converted into grayscale pixel values. Thus, the integrated fluorescence intensity (IFI) was measured and the red/green index was determined as the ratio of the two IFI. Using the minimal selected acid pretreatment (40 s, 1N HCl at 60°C), the DNA conformation indicated by red/green IFI was found to be shifted in the whole nucleus at 20 min and rather similarly at 60 min after HRG treatment (Fig. 4A). However, by introducing into the protocol a 5 min air-drying step before fixation, we were able to better separate both peaks: 20 min peak was discriminated by a smaller shift from the 60 min spectral shift (Fig. 4B). Using this composition of air-drying and 40 s HCl exposure, we could also separate a peak at 10 min after HRG treatment (Fig. 4C), which differed from the control with low confidence (p = 0.03) only but was clearly different from the more shifted 20 min peak (p < 0.01). Two 60 min peaks, the first practically coinciding with the 20 min one and another more shifted towards red fluorescence are also clearly seen (Fig. 4C). However, when we applied a drastic longer treatment by acid (5 min 1N HCl at 60°C), we could find a spectacular shift of the 20 min peak from the control peak, unequivocally indicating the conformational influence of the HRG action on DNA at this time point (Fig. 4D). The second peak of the more unfolded 60 min fraction of the transcriptionally active DNA became simply extracted (Fig. 4D). Contrary to it, the first 60 min peak showed its resistance to acid extraction. With this sharp provocation, the 10 min HRG fraction still remained not differentiated from the control (p = 0.24). So, our DNA probing approach (repeated in independent experiments) revealed that the cells underwent two clear DNA conformation shifts from unimodal to bimodal profiles after HRG action: the first bifurcation occurred between 10 to 20 min and the second between 20 to 60 min. The highly unfolded DNA fraction of the second peak after 60 min of HRG underwent extraction by prolonged acid pretreatment. The extraction of the active euchromatin by prolonged acid treatment (being here half-time of depurination) is known in DNA histochemistry and could be quite expectable [42]. As to splitting of the 60 min spectrum into two fractions, the first less sensitive to acid provocation than the second and more resistant even to harsh acid treatment, most likely reveals the subpopulation of cells, which were not involved in the second (after 20 min) phase of the HRG stimulation of transcription. Our findings well correspond to the qPCR data of the early gene response, with c-FOS activation peak at 20 min and sequential FOSL-1 and c-Myc peaking at 1 h ( Fig. 2A) and with two peaks of self-organizing criticality (SOC) in the whole transcriptome of HRG-treated MCF-7 cells revealed at 15-20 min by Tsuchiya et al [14,15] and at 60-75 min, by Chen and colleagues [16].
In literature, the stress-induced activation of the whole genome by transient long non-coding RNA (lncRNA) transcription from constitutive heterochromatin releasing silencing [43] and modulating the chromatin by lncRNA [44] were reported.
Therefore, due to our interest in PADs, we decided to check the synthesis of pericentromere lncRNA at the two-time points revealed by AO DNA conformation shifts.

The qPCR of the lncRNA synthesis from the human centromere satellite 3 DNA shows difference at 20 and 60 min after treatment with HRG
Testing with two probes for human centromere satellite 3 of the chromosomes 1 and 9 transcripts by qPCR, we could find a slight down-regulation at 20 min (p<0.05) and ~ 1.5-fold upregulation at 60 min (p<0.02) after treatment with HRG. These dynamics are displayed similarly with two chromosome-specific probes (Fig.2B). Together with the AO-test results, this data is another evidence of the reality and peculiarity of the chromatin SOC transition occurring between 15 to 20 min, found earlier by Tsuchiya with colleagues [14,15]. We then decided to study the possible change of the nuclear topology of pericentromere heterochromatin domains (PADs) just around this time point.
By using immunofluorescence (IF) and image analysis, we addressed if and how PADs size and area were changed after HRG treatment in comparison with the starving control.

The chromocentres revealed by H3K9me3 epigenetic mark contain centromeres identifying pericentromere-associated domains (PADs), which are highly heterogeneous in the cells of serum-starving (ST) control
The MCF-7 cells were stained by IF with an antibody against the histone H3 methylation site H3K9me3, a PAD-specific repressive mark. The typical fluorescence microscopy results in control cells (ST) are presented in Fig. 3C. For comparison, the pointillist pattern of the H3K9me3 distribution obtained by Single-Molecule Localization Microscopy is presented (Fig. 3D) indicating high-density regions compatible with PAD chromocenters. In the confocal section shown in Fig. 3E, H3K9me3 and the centromere protein CENPA were co-stained. The signals are clearly colocalized. These PADs show both, a tendency to densely arrange around the nucleoli and to separate in many small clusters with sometimes a few but often only one or two centromeres (Fig. 3E).
For the evaluation of PADs heterogeneity, we had to know the cell cycle of the (ST) control. Specimens stained for PADs with antibodies against H3K9me3 (for PAD counts) and counterstained with DAPI for DNA by measuring the nuclear integral fluorescence of the DAPI channel were analysed. The representative distribution performed with the aid of epifluorescence microscopy (X1000) and Image-Pro Plus analysis is shown in Fig. 5A. Most cells contain 2C DNA content (determined in ana-telophases of the control, not shown) which means that they are in the G0-G1 phase. About 5% show a pre-apoptotic aggregation of the chromocenters with some DNA loss (left box), while about 9-10% of the cells are in S-G2-phase since these cells being mostly delayed in the G0-G1-phase by serum starvation never fully stopped the cell cycle. Thus, the majority of cells (about 85%) reside in the 2C fraction and their chromocenter (PAD) numbers and areas are not a subject of DNA replication. Nevertheless, one can see following the distribution along the Y-axis of the major 2C fraction that these cells contain a very heterogeneous average number of centromeres per nucleus, the majority (middle-box) varying between 6 to 16 PADs per cell.
We further deliberately selected in three experiments in the population of ST cells those which possessed mostly small PADs and those with mostly large ones (excluding the frankly preapoptotic cells). In Fig. 5B, the relationship between the inverted number of PADs (Y-axis) and their perimeter derivative (X-axis) is presented. The high correlation (r = 0.76) of the two values shows that small and large chromocenters could arise only by splitting and fusing of the same material, corresponding to the constitutive, constant proportion of the pericentric heterochromatin.

The PAD sizes and numbers in the population and in individual cells are displaying the scale-free powerlaw distribution
The ST cell populations were further examined relating the average PAD number to their average area per cell in the epifluorescent microscopic study of 406 randomly selected cells (Fig.  6A) and in the confocal microscope by evaluating the same parameters of 567 individual PADs from 22 ST cells (Fig. 6B) using the software Image-Pro Plus. In both cases, we obtained a very similar power-law distribution. In general, the power-law (characterized by "a bursting synchronization") is displayed when the objects of different size transform one into another (e.g., by splitting-folding), organize networks, or/and in any case are driven by a common dynamical process like pulsing [45]. The typical power-law scale-free character of the relationship between size and number of PADs is clearly seen from the confocal microscopy result which resolved the minute ~0.4 µ m 2 PADs in addition to those larger ones obtained by epifluorescent microscopy (Fig. 6B, compare with Fig. 6A). The transition between the smallest PADs (~ 0.4 µ m 2 ) to the next (~1 µ m 2 ) and the following order (~ 2 µ m 2 ) is very steep ('bursting"); the final step (~ 4 μm 2 ) still shows a small and flatter located cell cluster. However, the larger ( > 4µ m 2 ) relatively scarce chromocenters are already on the plateau of both distributions.
The cluster analysis performed on the maximum projection images of confocal sections of ST cells estimating the number of centromeres by centromere identifier -antibody to CENPA protein [46] per PAD area stained for H3K9me3 (by Image-Pro Plus in interactive mode, with visual control) is presented in Fig. 6C. It shows that the size of the PADs is generally proportional to the number of the centromeres within (r = 0.84). Fig. 6C suggests that small PADs including 1 to 4 centromeres undergo dynamic transitions from one to another, with the kinetics of the power-all-or none-law, while a deterministic boundary creates the 'stopping condition' at 6 and more centromeres. This boundary undoubtedly belongs to the nucleolar organizers (NORs) of the acrocentric chromosomes stuck around the mostly centred nucleoli as can be seen from Fig. 3E. So, the dynamic chromocenters likely tend to radially unravel (pulse) from the nucleolus-associated heterochromatin. The distribution in Fig. 6C, in fact, corresponds to the distributions seen in Fig. 6A and B, obeying the power law. It is worth noting that induced by HRG in this model, the acceleration of the transcription speed for the gene avalanche group found by the analysis of the whole genome (22 277 genes), also displays a very similar distribution pattern with a critical point determined between 15 to 20 min (Fig. 6D, republished from [15]).

Unravelling of PADs after HRG stimulation occurs by a critical phase transition
Our next question was -how this heterogeneity of PADs in ST cells is changed after HRG application, in the dynamics of these first 15-30 mins, where we also have recorded the DNA conformational shift by AO test. The results of three experiments recording the PAD number and area (averaged per cell) after 15 and 30 min HRG application, showed the accumulation of the treated cells in the population compartment with small PADs -up to 2 µ m 2 (Fig. 7A), with reduction of the heterogeneity and correlation values (red trendline to be compared with the blue one). So, the larger PADs originating from the perinucleolar rim have unravelled into smaller units during this short period. The boxplots of both distributions of ST and HRG treated cells presented in Fig. 7B also shows the decrease of the average PAD area with much smaller variations as compared to the ST control. In Fig. 7C, where the parameters are presented for PAD diameters in ST and in HRG treated cells at 15 min and 30 min separately, we see that the comfort attractor (encircled) for the cells with small PADs is achieved already at 15 min HRG treatment and remains such at 30 min. It is likely dynamically reached through pair-wise fluctuations of the PAD units (see the boxed zone above). As judged from Fig. 6C, this level of dynamic PADs organisation likely enrols 4-2-2-1 centromeres.

Unfolding of the active H3K4me3-positive chromatin around the unravelled PADs
We also applied immunostaining of the two histone H3 modification markers H3K9me3 and H3K4me3 for the visualisation of the pericentric heterochromatin and active euchromatin, correspondingly (exemplified in Fig. 3F). We wanted to see whether the smaller size of the PADs

Discussion
In line with the early transcription change at two critical points, at 15-20 min found by Tsuchiya et al [14,15] and at 60-75 min revealed by Chen et al [16] in the MCF-7 breast cancer cell model, we found two DNA conformation changes between 10 to 20 min and at 1 h after HRG treatment expressed by two shifts from green to red fluorescence in the established AO-DNA structural test [36][37][38][39]. Our AO test showed that at 60 min the DNA conformation was more unfolded than that at 20 min HRG treatment. Our qPCR studies showed that the 20 min shift was coinciding with the peak of c-FOS activation, while the 60 min peak coincided with a peak of FOSL-1. The expression of FHL2 was opposite bi-phasic to FOSL-1; the annihilation of this early gene response occurred by 3h as established earlier [11,16]. Activation of c-Myc was also found at the 60 min time point (in one technical replicate of four in two experiments, it was more activated already at 20 min after HRG exposure). c-Myc is a master activator of the early response of bivalent genes and can contribute to the chromatin expansion by acetylation of histone H3 tails in the second euchromatin unfolding peak [48,49]. We also revealed that the 20 min peak of the DNA conformation change was not associated with activation of lncRNA transcription from centromeric human satellite 3 DNA, while the 60 min peak was.
We were interested, what really happens with PADs at the first critical phase transition involving a sandpile acceleration of transcription of several hundred genes [14,15,47] at 15-20 min after HRG treatment. Firstly, we have studied the behaviour of the constitutive PADs characterized by a repressive histone H3K9me3 mark in the control ST cells (mostly accumulated in the 2C G0+G1 phase) and found that the distribution of these domains according to size and number was obeying a scale-free, power-law. Moreover, the same power-law distribution was obtained when PADs were measured and averaged per cell in the whole cell population, or when individual PADs were measured in single cells; as well it was irrespective of the epifluorescent or confocal levels of microscopic resolution for small and smallest PADs. The power-law distributions and measurements of the perimeters of the large and small PADs indicated that these domains did not appear de-novo but convert one into another (splitting-fusing or rather folding and unravelling). With the same power, the area of PADs was proportional to the number of centromere clusters, while the larger clusters around nucleoli containing six and more centromeres put a boundary to their dynamics. In fact, this data corresponds to the earlier findings on the dynamism of the chromocenters of G1 cells [50].
When studying the HRG-stimulated samples, we verified the coincidence of the critical phase transition in the state of transcriptionally active chromatin at 15-20 min revealed by Tsuchiya et al [14,15]. The AO test applied here confirms the data by a critical decrease of PAD clusters to a size ≤ 2 µ m 2 (dominated by two or one centromeres). Interestingly, after HRG treatment the cells with smaller PADs formed a clear attractor at 15 min, which was sustained at 30 min. The transient zone of the "smeared" values seen in Fig. 7C suggests fluctuations between a pair-wise fusing-splitting of small PADs before lodging in the attractor. In this attractor, correlation and variability between the PADs is low (a zone of comfort) in contrast to the initial situation with the starving control, which is characterized by the high variability and high power-law correlation of the PAD size and number. As measured by the doubled smoothening of the active euchromatin around these unravelled PADs and the DNA-AO shift, this PAD attractor corresponds to the physical unfolding (phase transition) of the transcriptionally activated euchromatin.
Both events, a sandpile-like change of transcription rate with the unfolding of the euchromatin DNA and the splitting (unravelling) of PADs coincided and occurred by the same scale-free powerlaw with the all-or-none effect developed very steadily within 5 min [14,15,47]. Our observations are thus consistent with the classics of the self-organized criticality in dynamical systems with spatial degrees of freedom naturally evolving into a self-organized critical point [51]. Formation of the early HRG-induced attractor by PADs composed of smaller units, with low variation and correlation, corresponds with the physical principle of minimum energy in phase transitions [52] and observations of Gorban and colleagues on the system comforting in crisis [53]. Dynamic small chromocenters have also been described in ESC and in stress-response of cells [30]. Because of the restricted nuclear volume, DNA cooperative transitions were predicted by Yoshikawa to occur in a stepwise manner [19]. The experiments described here really confirmed this stepwise transitions since a pair-wise unravelling of PADs and two sequential conformational shifts of transcribing DNA are strongly correlated to cooperative DNA activities.
In addition, our finding by the centromere cluster analysis of the tentative NORs boundary of PADs unravelling (radial pulsing) is intrigue in accord with the studies on the relationship between NORs and chromosome territories as a global interactive network for cell nucleus functional organisation [54,55] involving intracellular phase transitions [56].

Conclusions
The effect of PADs co-segregating with repressive genomic regions and determining gene silencing in the proximity of centromeres (by position effect variegation) is long known and well established [27,28]. We hypothesize from the data obtained here that PADs participate in the early commitment of differentiation at 15-20 min after HRG stimulation by the physical splitting/unravelling the centromere clusters under the threshold of their negative position effect on transcription. Apparently, a certain critical mass of the repressive PAD clusters may be needed for the cooperative organisation of the epigenetically silenced domains (attracting the polycomb-group proteins and fusing with facultative H3K27me3-marked heterochromatin). Coming under a threshold by a pair-wise splitting of PAD clusters may cause the PAD-mass-dependent phase