Ultrastructural Visualization of 3D Chromatin Folding Using Serial Block-Face Scanning Electron Microscopy and In Situ Hybridization (3D-EMISH)

The human genome is extensively folded into 3-dimensional organization, yet the detailed 3D chromatin folding structures have not been fully visualized due to the lack of robust and ultra- resolution imaging capability. Here, we report the development of a novel electron microscopy method that combines serial block-face scanning electron microscopy with in situ hybridization (3D-EMISH) to visualize 3D chromatin folding at targeted genomic regions with ultra-resolution (5×5×30 nm in xyz dimensions, respectively). We applied 3D-EMISH to human lymphoblastoid cells at a 1.7 Mb segment of the genome and visualized a large number of distinctive 3D chromatin folding structures in high ultra-resolution. We further quantitatively characterized the reconstituted chromatin folding structures by identifying sub-domains, and uncovered a high level of heterogeneity in chromatin folding ultrastructures, suggestive of extensive dynamic fluidity in 3D chromatin states.

the techniques revealed smaller subunits of chromosome territories -topologically associated domains (TADs by Hi-C) [7][8][9] and chromatin contact domains (CCDs by ChIA-PET) 10 , demarcated by multiple binding sites for CTCF protein 11 . However, even though the contact probabilities are clearly displayed, the mapping data lack physical scale (e.g. in micrometers). To overcome these limitations, microscopy can be used to visualize the actual genome architecture and metric scale in various spatiotemporal resolutions in individual nuclei with different DNA staining methods. Standard fluorescence microscopy has an optical resolution of about 250 nm 12 and has been applied to image mammalian nuclei, successfully establishing the concept of chromosome territory 13 showing chromosomal morphologies in different cell cycle phases, and is also routinely used in cytogenetics to study abnormal chromosomes 14 . Recently, a super-resolution fluorescence light microscopy (20x20 nm in xy dimensions, 50 nm in z) with sequence-specific DNA binding probes was applied to visualize specific chromatin folding structures for a 10 -500 Kb target genomic region in Drosophila cells 15,16 and 1.2 -2.5 Mb target genomic region in human cells 17 . However, this method has a limit in z-dimension image depth (up to 3μm) 17 , and thus larger 3D chromatin structures beyond this limit would be truncated or lost. Electron microscopic in situ hybridization (EMISH) using biotin-labeled DNA probes coupled with diaminobenzidine staining has been used to image chromosomal DNA in the nuclei 18 . When standard FISH protocols are used in electron microscopy (EM) study, severe artefacts have been attributed to formamide and high temperature treatment 19  chromatin structure is EM tomography with even higher resolution (1x1x1 nm), in which photo-oxidized label is used to mark all DNA and to visualize the overall chromatin structures for 3D imaging in the nucleus 20 , but within a limited depth at z-axis (250nm). Nonetheless, our ability to visualize 3D chromatin folding structures in high quality ultra-resolution remains inadequate, hampered in particular by the lack of sequence specificity and low depth of imaging in the z-dimension.
To achieve ultra-resolution visualization of sequence-specific 3D chromatin folding structures, we developed 3D-EMISH, which combines advanced in situ hybridization using biotinylated DNA probes 21 with silver staining and serial block-face scanning electron microscopy (SBF-SEM) 22,23 . The serial z-stack EM images can be assembled computationally into a 3D conformation of the targeted genomic region with a resolution of 5x5 nm in the xy-plane and 30 nm in the z-axis. We applied 3D-EMISH to visualize a specific chromatin folding location in the human genome, and analyzed more than 200 distinctive 3D chromatin structures derived from individual nuclei of human lymphoblastoid cells. To our knowledge, this is the first time when SBF-SEM is applied for imaging specific chromatin folding structures.

The 3D-EMISH methodology
Our method includes 3D preservation of the nuclei, in situ hybridization of specific DNA probes, serial scanning by EM, and imaging data analysis to reconstruct the spatial models of chromatin folding structures (Fig. 1a).
In the first step of the 3D-EMISH protocol, cells were fixed with 4% freshly-made paraformaldehyde, and then embedded in a thrombin-fibrinogen clot, postfixed, cryoprotected in 30% sucrose, and frozen in dry ice-cooled isopentane. The embedded clots (less than 5mm in size) were cut into 40 µm-thick sections in the freezing microtome.
The next step in 3D-EMISH is in situ hybridization (ISH). The free-floating 40 µm-thick sections were incubated with the biotinylated DNA probes targeting a specific genomic region of interest. The signals of the biotinylated DNA probes for EM detection were processed with the use of 1.4nm-thick streptavidin-conjugated fluoronanogold particles and subsequent silver enhancement. The procedure of silver enhancement was experimentally optimized to obtain electron dense 4-5 nm diameter granules as the smallest items ( Supplementary Fig. 1a-c). To optimize chromatin ultrastructure preservation, we performed a set of control in situ hybridization (ISH) reactions. Using transmission electron microscopy (TEM), we assessed in details the influence of potentially harmful factors associated with the ISH procedure. We found, that the permeabilization with Triton X-100 causes negligible changes in chromatin structure (Compare the Supplementary Fig. 1d and 1e) . Formamide and high temperature did evoke some damage in the form of small, empty spaces in the nucleoplasm and a significant diminution of the condensed chromatin layer beneath the nuclear envelope ( Supplementary Fig. 1f). However, surprisingly, we found that the factor most significantly disturbing the chromatin ultrastructure was the incubation of the sample with dextran sulfate ( Supplementary Fig. 1f Fig. 1hi). However, it led to somewhat increased background (especially at the surface of the sliced section). The 40 µm-thick sections were stained with uranyl acetate and tannic acid, followed by dehydrating and flat-embedding in an epoxy resin. Because uranyl acetate stains nucleic acids and proteins broadly, its inclusion in 3D-EMISH enabled us to visualize the overall nuclear space in relation to the target chromatin structure (Fig. 1a).
One block of the embedded specimens (usually less than 1mm in width) was placed in the chamber of a SBF-SEM (Supplementary methods, Zeiss Sigma VP) equipped with a built-in ultramicrotome (Gatan 3View). The specimen block was consecutively sliced one layer at a time at 30-50nm intervals, and the exposed surfaces of the specimen were serially scanned in a field of 8192x8192 pixels (approximately 1700 µm 2 , when a pixel size is 5 nm) to obtain volumetric data including the specific 3D-EMISH signals. The specimen was sliced again for the next round of signal acquisition and so forth. This cycling process of slicingscanning was performed hundreds of times to generate z-stack images in a 3D-EMISH experiment (Fig. 1a).
Since all in situ hybridization reactions result in non-specific background 21 , the 3D-EMISH signals in each slice also included punctate patterns of staining that we deemed likely to be noise. We reasoned that the real signals from the targeted chromatin would be in continuation in multiple layers of the Z-stack images, whereas the background signals would not be connected in the Z-stacks (Fig. 1b). Therefore, we developed an image processing algorithm based on multi-layer connectivity to retain true chromatin positive signals and remove false-positive ones. After filtering out the background noise, the chromatin signals were assembled to reconstruct the 3D chromatin folding structure in 3D EM resolution for the targeted genomic region (Fig. 1c-d).

3D EM resolution of chromatin folding structures captured by 3D-EMISH
We applied 3D-EMISH to investigate chromatin folding structures in the genome of the To confirm the presence of the chromatin loop domains inferred by these sequencing-based approaches, we designed and produced biotinylated DNA probes derived from 11 BAC clones to cover this region for 3D-EMISH experiments (Fig. 2a). To test the hybridization efficiency of the probes to the target genomic region, we also generated two-color fluorescence probes from the BAC clones that could help distinguish the three domains ( Fig.   2a). Subsequently, we performed 3D FISH imaging analysis using confocal microscopy, and validated that these probes were efficient and highly specific (Fig. 2b).
Using the biotinylated DNA probes, we performed 3D-EMISH experiments in GM12878 cells at exponential growth phase. We performed two independent experiments using two different EM settings, producing image voxel size of 7x7x50 nm in replicate 1 and 5x5x30 nm in replicate 2. The sample cubical blocks were analyzed by a series of 300 (replicate 1) or 600 (replicate 2) slicing and scanning cycles, rendering consecutive 50 nm (replicate 1) or 30 nm (replicate 2) sections, respectively. One specimen, typically about 1000x1000x40 µm, included tens of thousands of cells. However, the cubic volumes that were scanned by EM microscopy were much smaller and they were 3200 µm 2 in xy and 15 µm in z for replicate 1, and 1700 µm 2 in xy and 18 µm in z for replicate 2. Each of the Z-stack images could capture a cross section of about 10 -30 cell nuclei, intact or truncated (Fig. 2c).
After signal de-noising and assembling, the outline of a nuclear framework and the specific chromatin objects could be visualized respectively based on the uranyl acetate and the silver staining signals (Fig. 2d). From two independent 3D-EMISH experiments, a total of 166 nucleus image-stacks were obtained. Most of them (140) were truncated nuclei and contained from 1 to 4 specific target signals, whereas 26 nuclei were intact with 2 or 4 specific chromatin targets (Fig. 2e), suggesting that the cells were possibly in G1 phase (2 copies of target region) or S-G2 phase of the cell cycle.
To compare 3D-EMISH with super-resolution microscopy, we applied iPALM 24 and visualized the same targeted chromatin structures using iPALM-specific probes tagged with two-color "blinking" fluorophores for the same genomic region in GM12878 cells ( Supplementary Fig. 2). The iPALM 3D image structure showed about 4-fold lower resolution (20x20 nm in xy) than the 3D-EMISH images. In addition, the depth of an iPALM image at z-dimension was limited at 750nm 25 , thus capturing only incomplete chromatin structures, even though isolated nuclei were used for iPALM imaging to reduce the depth of cell specimen. Thus, 3D-EMISH not only provided higher resolution of chromatin images at xy plane than iPALM, but also provided greater depth at the z-axis than iPALM, enabling visualization of complete 3D chromatin folding structures. However, the use of two different fluorophores in iPALM allowed determination of the orientation of the chromatin folding structure along the linear DNA strand. This representation allows clear dissection of the chromatin signals in more detailed folding structures (Fig. 3a). Next, we calculated the local density centers in each chromatin structure, followed by diffusing the density from the centers to demarcate the boundaries and to identify distinctive domains in the 3D-EMISH structures (Fig. 3a). For example, in the three different 3D-EMISH structures (sID50, sID12, sID42), one, two and three local density centers were identified, respectively, corresponding to morphological domains ( Fig. 3b and Supplementary   Fig. 4). The sub-domains in each 3D-EMISH structure were demarcated by different colors (Fig. 3c).
Only two structures were identified as five domains, and we merged them with four domains structure group for our statistical analysis. Remarkably, the structures with multiple sub-domains accounted for a combined 75% of all the structures, which is approximately in line with the ChIA-PET mapping data (Fig. 2a). To further characterize these chromatin folding structures, we measured the volume and the surface area for each. These measurements showed that the chromatin structures with one domain 1 1 showed increased values along with the numbers of domains (Fig. 4c-d). We also calculated the form factor (volume/surface; see the detailed formula in Online Methods) for each structure, which showed the same trend as the other measurements, i.e. single-domain structures had the lowest form factor value (Fig. 4e).
Taken together, our imaging analyses suggest that the small percentage of one-

Discussion
We have demonstrated 3D-EMISH as an effective imaging technique with ultra-resolution Moreover, studies of transcription reveal expression heterogeneity (cit). Also in the same cells that we studied -GM12878, heterogeneity was shown using single cell sequencing approaches -scATAC-Seq and scRNA-Seq when compared to bulk studies. It suggests that variability observed by us and others could interplay with transcriptional activity.
Nevertheless, the relation between dynamic changes of chromatin structure and stochasticity in gene expression is not yet fully understood (cit).

DNA probes used in in situ hybridizations
The

Confocal microscopy imaging and post-processing
The samples were imaged with Zeiss LSM 780 microscope, using a 63x oil immersion

Two-color DNA-FISH experiment -iPALM
DNA probe signals were captured by one-step system detection of biotin-labelled probe with streptavidin-CAGE590 (Abberior), and one-step system detection of digoxigenin-labelled probe with anti-digoxigenin-Alexa Fluor 647 (Jackson ImmunoResearch).

iPALM imaging and post-processing
The samples for iPALM were imaged in STORM buffer 28 , using coverslips with fiducial markers, prepared according to Shtengel et al. 25 . 3D-FISH samples were imaged using noncommercial iPALM system 24 available at Advanced Imaging Center-HHMI's Janelia Research Campus, with excitation: 647nm laser (5.5 kW/cm 2 ) and 561nm laser (4.5 kW/cm 2 ) and activation: 405nm laser (100W/cm 2 ). The number of frames is 25000, and the exposure time is about 50ms. iPALM data was processed to extract single molecule 3D coordinates and localization precisions, and subsequently visualized using the PeakSelector software (Janelia Research Campus) 25 . FISH signals acquired by iPALM were segmented based on the threshold with custom-made software PartSeg (https://4dnucleome.cent.uw.edu.pl/PartSeg/). Then, 3D reconstructions were obtained using Imaris (Bitplane).

Cell preparation for 3D-EMISH
The experiments were performed on cultured GM12878 cells. After centrifugation, fixation of the cell pellet (~7x10 7 cells) with 4% paraformaldehyde in PBS overnight at 4 C and washing with PBS, the artificial tissue with embedded cells was formed. To the fixed and washed pellet, 0.2 ml the fibrinogen solution (100 mg fibrinogen (Sigma; F3879), 53.33 mg sodium citrate, 283.33 mg sodium chloride, 33.33 ml H 2 0) was added and stirred. After 1 8 centrifugation, 0.2 ml the thrombin solution (Sigma; T7009-100UN + 1ml H 2 0) was added 29 .
The clot was postfixed with 4% paraformaldehyde in PBS (10 min; RT) and washed with PBS.
Then, the formulated artificial tissue was soaked with 30% sacharose in PBS overnight at 4 C and frozen in Tissue Freezing Medium (Leica). After overnight storage at -80 C, the clot was cut into 40 µm slices using cryostat (Leica).

In situ hybridization for 3D-EMISH
In situ hybridization was performed on the free floating sections of the clot according to

3D-EMISH staining
Sections after passing the quality control using confocal microscopy were washed with ddH 2 0 with silver conductive resin (CircuitWorks), and incubated 2 hours at 65 C. Using ultramicrotome (Leica) and diamond knife (DiATOME) smooth surface of sample to SBF-SEM was obtained. Then, to enhance conductivity, sides of the sample were painted with silver paint (Ted Pella).

EM imaging
To image 3D-EMISH samples, Zeiss Sigma electron microscope with backscatter electron detector and ultamicrotome with diamond knife inside the chamber was used (3View2 from Gatan). All images were collected using variable pressure mode with EHT 4-6kV, aperture 30 µm and resolution 8192x8192 pixels. Replicate 1 was collected with pixel size 7 nm, sliced each 50 nm, 300 times (voxel size: 7x7x50 nm). Replicate 2 was collected with pixel size 5 nm, sliced each 30 nm, 600 times (voxel size: 5x5x30 nm). To examine potential harmful effect of reagents on nuclear ultrastructure, transmission electron microscope JEM 1400 (JEOL) was used.

Image processing and quantitative analysis
All 3D electron microscopy stacks were preprocessed with ImageJ plugin -Linear Stack Alignment with SIFT 31 . The images were manually inspected, and then a cuboid containing the region of interest (ROI) was cut. The mask defining the ROI was determined by connecting component analysis of objects inside the cuboid. We took the largest connected component using custom made software, where we used Maximum Entropy algorithm (ImageJ plugin) to optimize the threshold value to segment the structure. The resulting mask 1 was overlaid on the image to segment the structure. Subsequently, 3D images were resampled in order to obtain the same voxel size (5 nm) in x, y, and z. The reslicing was performed by upsampling with linear interpolation between adjacent planes. We also applied a gaussian filter of size 1 pixel (i.e. 5 nm) in the x-y planes in order to remove possible pixel noise. The isotropic scale was required by the plugin used to produce the movies (with unequal scales the structures appear in the movies unnaturally flattened), and also it was needed by the algorithm for morphological domains separation, which operates on a cubical grid. The upsampling is not obviously visible unless we rotate the structure through an angle perpendicular to the axial direction, in the supplementary movies the structures were rotated through an axis perpendicular to its longest (principle) axis, which in general is randomly oriented with respect to the slicing direction Then, Gaussian fliter (5 nm cut-off size) was applied to each z-section to remove the detector noise. The volume (V) and surface (S) of the segmented structure were calculated by 3D Object Counter plugin (ImageJ). The form factor (ff) was defined as where the value is 1 for a spherical object and increase as the structure is less compact than a sphere. For illustrative purpose, the segmented structures were aligned accordingly to their principal axis 32 . The structures were classified into 5 different morphological groups: 1 compact sub-domain group (one spherical ball shape), 2 distinctive sub-domain group (dumbbell shape), and 3 -5 distinctive sub-domain groups by applying the following subdomain identification algorithm.
First, the structure was smoothed out by applying the "maximum filter" with 135 nm cut-off size. This parameter was experimentally tested and determined to prevent from oversegmenting sub-domains of a structure , yet preserving the overall morphological feature.
Second, the local maximal density centers were calculated per each identified distinctive group. Third, these local maxima were subsequently used as seed points in computing of the diffusion process, where the image density was used as a local diffusion coefficient. The diffusion process was simulated by solving numerically an uncoupled set of diffusion equations, denotes the density of the diffusing material, and D denotes the EMISH signal density, i=1...N (the number of seed points) , r Ԧ are the pixel coordinates in 3D space.
Initially, we assume that all EMISH signals of each sub-domain was concentrated in each seed point. We represented it as normalized gaussian functions centered at the seed points.
Image voxels, diffused from the same seed point, belong to the same sub-domain group. The 3D image region of low EMISH signal density will diffuse slower than that of high EMISH signal density. The boundary voxels had diffused signals from two or more different seed points, and these boundary voxels belong to the sub-domain group which is the highest diffused signal. The diffusion process was completed when all voxels were filled by a nonzero density and the structure was divided into separate sub-domains.

Statistical analysis
Statistical analysis was performed using IBM SPSS Statistics software. In all cases, error bars indicate SEM. The morphological feature statistics of 3D-EMISH structures for different sub-domain groups was performed using Kruskal-Wallis test with Bonferroni post-hoc test.
For the three-dimensional ROI, we created the following binary masks: specific signal mask (SP), non-specific background signal mask (NS) and no-signal background mask (BN). The specific and non-specific signal volume, was obtained by summing all mask elements for SP and NS masks respectively, and multiplying it by the voxel volume. The content of the unspecific signal in the background volume was obtained as a ratio of NS mask volume, and the total background volume (sum of NS and BN mask volumes).

Code and Data availability
3D-EMISH image processing code and data files are available at the following public repository server: