Cohesin prevents cross-domain gene coactivation (2024)

This research complies with all relevant ethical regulations by Environment, Health and Safety from Howard Hughes Medical Institute (HHMI)/Janelia Research Campus.

Chemicals

The plant auxin analog indole-3-acetic acid sodium salt (Millipore Sigma, I5148) was dissolved in double-distilled water to a stock concentration of 500 mM, aliquoted and stored at −20 °C and used at a final concentration of 500 µM (acute depletion) or 5 µM for 6 h. α-Amanitin (Tocris, 4025) was dissolved at 1 mg ml^–1 in double-distilled water, aliquoted and stored at −20 °C and used at a final concentration of 100 µg ml^–1.

Cell culture

JM8.N4 mouse ES cells from the C57BL/6N strain and their genome-edited derivatives were routinely cultured in 60-mm plates coated with 0.1% gelatin without feeders at 37 °C and 5% CO₂. The mouse ES cell culture medium was composed of optimized knockout DMEM for mouse ES cells (Thermo Fisher Scientific, 10829-018), 15% ES cell-qualified fetal bovine serum (ATCC, SCRR-30-2020), 1,000 U of leukemia inhibitory factor (self-purified), 1 mM GlutaMAX (Thermo Fisher Scientific, 35050-061), 0.1 mM MEM non-essential amino acids (Thermo Fisher Scientific, 11140-50), 0.1 mM β-mercaptoethanol (Thermo Fisher Scientific, 21985-023) and antibiotic–antimycotic (Thermo Fisher Scientific, 15240-062). 2i inhibitors were also added to the medium at final concentrations of 1 μM PD0325901 (Millipore Sigma, PZ0162) and 3 μM CHIR99021 (STEMCELL Technologies, 72052), respectively.

Generation of stable cell lines and genetic knock-in cells

We constructed plasmids PiggyBac-EF1-HaloTag-OCT4, PiggyBac-EF1-H2A.Z-HaloTag and PiggyBac-EF1-HaloTag-mouse TBP (mTBP) based on an available PiggyBac-EF1-HaloTag backbone vector. We used them together with the PiggyBac supertransposase to generate stable cell lines that expressed HaloTag–OCT4, H2A.Z–HaloTag and HaloTag–mTBP fusion proteins for single-particle tracking experiments. For H2B single-particle tracking experiments, we adopted a previously constructed PiggyBac-EF1-H2B-HaloTag plasmid⁵⁵. Maps of plasmids used in this study are available upon request.

To generate stable HaloTag–OCT4- and HaloTag–mTBP-expressing mouse ES cells, 8 µg of PiggyBac-EF1-HaloTag-OCT4 or PiggyBac-EF1-HaloTag-mTBP was nucleofected with 8 µg of PiggyBac supertransposase into ~1 × 10⁶ mouse ES cells. The derived mouse ES cells were selected with 500 µg ml^–1 G418 (Thermo Fisher Scientific, 10131035) for ~5 days, stained with JF549 dye and sorted by fluorescence-activated cell sorting (FACS) for cell populations with intermediately expressed genes.

We have applied previously described methods to generate single guide RNA (sgRNA) constructs and to construct HaloTag–MED1, MED6–HaloTag and HaloTag–BRD4 donor plasmids for generating knock-in cells by CRISPR–Cas9 genome editing¹⁸. The corresponding Emerald versions of MED1/MED6 donor constructs were generated by replacing the HaloTag DNA coding sequence with Emerald coding sequence by Gibson assembly (New England Biolabs, E5520S). We generated HaloTag–MED1, MED6–HaloTag and HaloTag–BRD4 genetic knock-in cells based on a previously generated Rad21 mAID degron ES cell line¹⁸. Briefly, 1.0 µg µl^–1 SpCas9-sgRNA-PGK-Venus construct and 1 µg µl^–1 donor construct were nucleofected into ~3 × 10⁶ Tir1 mouse ES cells using Amaxa 4D-Nucleofector and P3 Primary Cell 4D-Nucleofector X kits (Lonza, V4XP-3024) following the manufacturer’s protocols. Three days following nucleofection, cells were then stained with 50 nM JF549 HaloTag ligand (JF549) for 30 min, washed three times in 1× PBS and then in mouse ES cell medium for 15 min and subjected to sorting by FACS. JF549⁺ cells were plated sparsely in 60-mm tissue culture plates and grown for another 5~7 days. Single colonies were picked for genotyping by designing PCR primers outside of the hom*ology arms. Biallelic knock-in mouse ES cell clones were verified by PCR genotyping, Sanger sequencing and western blotting.

Western blotting

Mouse ES cells were lysed in 1× SDS sampling buffer (200 mM Tris-HCl (pH 7), 10% glycerol, 2% SDS, 4% β-mercaptoethanol, 400 mM DTT and 0.4% bromophenol blue). Lysates were sonicated and incubated on ice for 30 min, mixed with 2× loading buffer and denatured at 95 °C for 5 min. Lysed proteins were resolved by SDS–PAGE using Mini-PROTEAN TGX Precast Gels (Bio-Rad). The following primary antibodies were used: rabbit polyclonal anti-RAD21 (D213; Cell Signaling, 4321, 1:1,000 dilution), rabbit polyclonal anti-MED1 (CRSP1/TRAP220; Bethyl Laboratories, A300-793A, 1:1,000 dilution), rabbit polyclonal anti-MED6 (Abcam, ab220110, 1:1,000 dilution), rabbit monoclonal anti-BRD4 (E8V7I; Cell Signaling, 54615, 1:1,000 dilution), rabbit polyclonal anti-TBP (Cell Signaling, 8515, 1:1,000 dilution), rabbit polyclonal anti-OCT4 (Abcam, ab19857, 1:2,000 dilution), anti-histone H2A.Z (EPR18090; Abcam, ab188314, 1:2,000 dilution) and rabbit monoclonal anti-α-tubulin (11H10; Cell Signaling, 2125, 1:2,000 dilution). We used horseradish peroxidase-linked anti-rabbit IgG (Cell Signaling, 7074) secondary antibodies at a dilution of 1:1,000. Western blots were exposed by using the 20× LumiGLO chemiluminescent detection system (Cell Signaling, 7003) and imaged by using a Bio-Rad ChemiDoc MP detection system.

Cell cycle analysis

Cell cycle analysis was performed by using propidium iodide staining following the protocols from the propidium iodide flow cytometry kit (Abcam, ab139418). Briefly, cells were trypsinized into single-cell suspensions and fixed with 67% ice-cold ethanol in PBS overnight at 4 °C. The next day, the cells were rehydrated with PBS, stained with propidium iodide (final concentration of 50 μg ml^–1) and treated with RNaseA (final concentration of 50 μg ml^–1) for 30 min at 37 °C before flow cytometry analysis. All samples were acquired on a Beckman Coulter CytoFLEX S with four lasers (405 nm, 488 nm, 561 nm and 638 nm) and operated by CytExpert Software v2.3 (Beckman Coulter) using 488-nm FSC-H (488-FSC-H) as the threshold parameter (threshold automatic setting). The detector’s gain for fluorescence (CD4-FITC) and FSC/SSC detection were optimized by using control cells without treatment. The following gains were used for detectors with different spectral filters: 39, 15, 28 and 10 arbitrary units for 488-FSC, 488-SSC, 561-585/42 and 561-610/20, respectively. SSC-A versus FSC-W was used for initial gating of singlet cells, followed by SSC-A versus FSC-A to further define cellular events. An event count versus 561-610/20 histogram plot was used to determine the percentage of cells in G1, S and G2/M phases of the cell cycle. All samples were acquired for 5 min at a sampling rate of 30 μl min^–1 or up to 15,000 cells. FlowJo v.10.7.1 (Flowjo) was used for analysis of the flow cytometry data.

Smart-SCRB data acquisition

Mouse ES cells with an engineered AID system treated with or without auxin were resuspended in culture medium without Phenol red and sorted by FACS into 96-well PCR plates containing 3 μl of mild lysis buffer (nuclease-free water with 0.2% Triton X-100 + 0.1 U μl^–1 RNase inhibitor (Lucigen, 30281-2)). A total of ~400 cells were sorted and collected for each sample. The PCR plates with collected cells were briefly centrifuged, immediately frozen and stored at −80 °C until cDNA synthesis.

One microliter of harsh lysis buffer (50 mM Tris (pH 8.0), 5 mM EDTA (pH 8.0), 10 mM DTT, 1% Tween 20, 1% Triton X-100, 0.1 g l^–1 proteinase K, 2.5 mM dNTPs and ERCC Mix (10⁷-fold dilution)) and 1 μl of 10 mM barcoded RT primer was added to each well. Plates were incubated at 50 °C for 5 min to lyse cells, and proteinase K was heat inactivated by subsequently incubating at 80 °C for 20 min. To minimize contamination across wells, heavy-duty plate seals and quantitative PCR (qPCR) compression pads (Thermo Fisher Scientific, 4312639) were used to seal the plates. The lysis reaction was mixed with 2 μl of reverse transcription master mix 5× buffer (Thermo Fisher Scientific, 11756500), 2 μl of 5 M betaine (Sigma-Aldrich, B0300-1VL), 0.2 μl of 50 mM E5V6NEXT template switch oligonucleotide (Integrated DNA Technologies), 0.1 μl of 200 U μl^–1 Maxima H-RT (Thermo Fisher Scientific, EP0751), 0.1 μl of 40 U μl^–1 NxGen RNase Inhibitor and 0.6 μl of nuclease-free water (Thermo Fisher Scientific, AM9932). The reaction system was then incubated at 42 °C for 1.5 h, followed by 10 min at 75 °C to inactivate reverse transcriptase. PCR was performed by adding 10 μl of 2× HiFi PCR mix (Kapa Biosystems, 7958927001) and 0.5 μl of 60 mM SINGV6 primer and running the following program: 98 °C for 3 min, 20 cycles of 98 °C for 20 s, 64 °C for 15 s and 72 °C for 4 min and a final extension step of 5 min at 72 °C.

Single-cell cDNA was pooled by plate to make libraries. cDNA (600 pg) from each sample plate was used in a modified Nextera XT (Illumina, FC-131-1024) library preparation but using the P5NEXTPT5 primer and a tagmentation time of 5 min. The resulting libraries were purified following the Nextera XT protocol (0.6× ratio) and quantified by qPCR using a Kapa Library Quantification kit (Kapa Biosystems, KK4824). Four plates were pooled together on a NextSeq 550 high-output flow cell or NextSeq 2000 P2-100 flow cell with 26 bp in read 1 and 50 bp in read 2. PhiX control library (Illumina) was spiked in at a final concentration of 7.5% to improve color balance in read 1. Read 1 contains the spacer, barcode and unique molecular identifier and read 2 represents a cDNA fragment from the 3′ end of the transcript. The entire experimental procedure was replicated two more times, and ~400 cells were analyzed for each condition (control versus auxin treated).

Alignment and count-based quantification of single-cell data were performed by removing adapters, tagging transcript reads to barcodes and UMIs and aligning the resulting data to the mouse genome mm10. After quantification, 2,600 detected genes (gene_det > 2,600) was set as a threshold to eliminate unreliable sequencing results, and the data matrix (M_RNA) for all remaining cells was used for analysis.

ACD calling

To identify ACDs across the genome, we referred to available bulk ATAC-seq data¹⁸ and used a genomic distance of 200 kb as the binning unit for analysis. The MATLAB function findpeaks() was then called in ‘MinPeakProminence’ mode, and the value of ‘MinPeakProminence’ was set to the average of overall ATAC signals. In total, 776 ACDs were identified by using these criteria. The reason for using a large bin is to make sure that we generate a dense matrix without too many zeros. This is essential for robust downstream coexpression and coaccessibility calculations, as normal statistics break down when dealing with a sparse matrix with a lot of zeros.

To align ACDs with genomic features including A/B compartments, LADs and ChromHMM regions, the mouse genome (mm10) was divided into 500-bp windows using the ‘tileGenome’ function within the GenomicRanges package. These windows were then overlapped with multiple genomic features, including mouse ES cell bulk ATAC-seq peaks¹⁸, LAD regions²⁷, A/B compartments²⁶ and ChromHMM states²⁸, and annotated accordingly. To assess the enrichment of ATAC-seq peaks within each of these annotated groups, we calculated the fold change as the ratio of ATAC-seq peaks in each group compared with the ratio across the entire genome. Groups exhibiting a fold change of greater than 1 were considered to be enriched with ATAC-seq peaks.

Calculation of normalized contact frequency per ACD pair

To calculate the normalized contact frequency for a specific ACD pair, we used published mouse ES cell Hi-C data and calculated the average value over all the contact frequencies within the boundaries of the ACD pair¹⁸. The average value was defined as the normalized contact frequency for that specific ACD pair.

Gene coexpression analysis

To evaluate Spearman’s rank-order correlation among gene pairs, we filtered out low-expressed genes from M_RNA with a threshold of 1 after averaging the RNA-seq counts over the single cells analyzed. For each pair of expressed genes, the Spearman correlation coefficient S_x,y (considering tied ranks) was calculated by

x_i and y_i are the sequencing counts for genes x and y in the ith individual cell, respectively. For each ACD pair (ACD1 and ACD2), the Spearman correlation matrix S was determined by calculating the Spearman correlation coefficient for each gene pair from ACD1 and ACD2, and the average Spearman correlation coefficient \({\bar{S}}_{{\rm{RNA}}}\) was calculated by taking the average value across the whole matrix S for control and cohesin-depleted conditions.

The differential Spearman correlation matrix ΔS_RNA was calculated by

where S_RNA,RAD21– and S_RNA,Ctrl are Spearman correlation matrices for cohesin-depleted and control conditions, respectively. The workflow for the above analyses is illustrated in Extended Data Fig. 2a.

The gene coexpression coefficient ΔC_RNA(i,j) between two ACDs within the same chromatin was estimated by using the following equation:

where n_ΔS+ and n_ΔS– represent the number of positive and negative elements within the differential Spearman correlation matrix ΔS_RNA. The derived differential gene coexpression matrix ΔC_RNA was plotted as a heat map. The workflow for evaluating ΔC_RNA is illustrated in Extended Data Fig. 2b.

To compare the cis (within chromosome) and trans (between chromosomes) effects of gene coexpression per ACD pair, we calculated chromosome-wise differential coexpression coefficients by binning differential coexpression coefficients (ΔC_RNA(i,j)) in cis or within one chromosome pair in trans. The results were plotted in a heat map or bar plot (the bar for trans effects is an averaged value) in Extended Data Fig. 3c,d.

Cell cycle phase classification based on Smart-SCRB data

To dissect the effect of cell cycle stages on gene coexpression analysis, we adopted a computational method described by Scialdone et al. for classifying cells into cell cycle phases based on Smart-SCRB data⁵⁶. Using a reference dataset, the difference in expression between each pair of genes was computed. Pairs with significant changes across cell cycle phases were selected as markers for classification and applied to a test dataset. Cells were then classified into the appropriate phase based on whether the observed change for each marker pair was consistent with one phase or another. This approach was implemented in the cyclone() function from the ‘scran’ package, which contains a pretrained set of marker pairs (mouse_cycle_markers.rds) for mouse cells.

scATAC-seq data acquisition

Mouse ES cells treated with or without auxin were washed and resuspended in 1× PBS with 0.04% bovine serum albumin. Nuclei isolation for scATAC-seq from cell suspensions was performed according to the manufacturer’s demonstrated protocol (CG000169, Rev E, 10x Genomics). Nuclei were counted using a Luna II automated cell counter (Logos Biosystems). Approximately 15,000 nuclei per sample were loaded and subjected to a Chromium NextGem scATAC-seq v2 assay (10x Genomics). The resulting libraries were sequenced on a NextSeq 2000 (Illumina; 50 bp read 1, 8 bp i7 index read, 16 bp i5 index read and 50 bp read 2).

FASTQ files of raw data were processed by using the CellRanger ATAC (10x Genomics, v2.1.0) analysis pipeline. Reads were filtered and aligned to mouse genome mm10 (10x Genomics, refdata-cellranger-arc-mm10-2020-A-2.0.0) using the cellranger-atac count() function with default parameters. The barcoded and aligned fragment files were then loaded by ArchR (version 1.0.1). Low-quality cells with minimum transcription start site enrichment scores of less than 4 and minimum fragment numbers of less than 1,000 were filtered out. Doublets were inferred by the addDoubletScores() function and removed using the filterDoublets() function with default parameters.

Cross-ACD chromatin coaccessibility analysis

To evaluate the chromatin coaccessibility per ACD pair, we binned ATAC counts within each ACD to derive an ACD counts matrix (M_ATAC) for each experimental condition. After adding ACD regions using the addPeakSet() function, counts for each ACD per cell were aggregated together using addPeakMatrix() with a very high ceiling value (1,000). To correct batch effects, we downsampled (50%) the scATAC-seq data to match the distribution of reads per ACD count between experimental conditions, and ten randomly downsampled dataset pairs were used for the calculations below. The ACD counts matrix (M_ATAC) for each condition was normalized to the average value \({\bar{M}}_{{\rm{ATAC}}}\) of the matrix to derive a normalized ACD counts matrix (N_ATAC). To mitigate the influence from poorly sequenced cells (have 0 count for many ACDs), each column j of N_ATAC with mean value \({\bar{N}}_{{\rm{ATAC}}}\left(j\right)\ge {\bar{M}}_{{\rm{ATAC}}}\) was selected and integrated to form a new matrix N_{ATAC-Filtered} for Spearman’s rank-order correlation evaluation. The final results from ten downsampled dataset pairs were pooled and averaged.

For each ACD pair, the Spearman correlation coefficient for chromatin coaccessibility was calculated by using Eq. (1). The Spearman correlation matrix S_ATAC was determined by calculating the Spearman correlation coefficient for each ACD pair across the genome. The differential Spearman correlation matrix ΔS_ATAC was calculated by

where S_{ATAC,RAD21–} and S_ATAC,Ctrl are Spearman correlation matrices for cohesin-depleted and control conditions, respectively. S_{ATAC,RAD21–}, S_ATAC,Ctrl and ΔS_ATAC were plotted as heat maps. The workflow for evaluating the Spearman correlation coefficients for chromatin coaccessibility is illustrated in Extended Data Fig. 4a.

smRNA-FISH

smRNA-FISH probe blends were designed through Stellaris RNA-FISH probe designer (Biosearch Technologies) and include 40~48 serial probes that target only the intron regions of selected genes. The probe pairs designed for gene pairs selected from neighbor TADs were labeled by Quasar 570 and Quasar 670, respectively, for two-color experiments. The commercially synthesized oligonucleotide probe blends were dissolved in 400 μl of TE buffer (10 mM Tris-HCl and 1 mM EDTA, pH 8.0) to create a probe stock of 12.5 μM.

smRNA-FISH experiments were performed according to the Stellaris RNA-FISH protocol for adherent cells provided by Biosearch Technologies. Specifically, cells were grown on 18-mm round number 1 cover glass (Warner Instruments, CS-18R) in a 12-well cell culture plate coated with human recombinant laminin 511 (BioLamina, LN511-0202), fixed in 1× PBS with 3.7% formaldehyde (Millipore Sigma, F8775-25ML) for 10 min and permeabilized in ice-cold 70% (vol/vol) ethanol for 2 h. The coverslips with cells were immersed in 100 µl of hybridization buffer (90 µl of Stellaris RNA-FISH Hybridization Buffer (Biosearch Technologies, SMF-HB1-10) and 10 µl of deionized formamide (Millipore Sigma, S4117)) at a final probe concentration of 125 nM and placed into a humidified chamber. The assembled humidified chamber was incubated overnight at 37 °C before the sample coverslips were washed, co-stained with 5 ng ml^–1 DAPI (Sigma-Aldrich, D8417) and mounted onto slides for future imaging analysis.

Chr 2 intron-FISH probe design, synthesis and amplification

For 208 active genes across mouse Chr 2, we selected non-overlapping 35-nucleotide (nt) probes with several constraints, including a maximum melting temperature of 100 °C, minimum melting temperature of 74 °C, secondary structure temperature of 76 °C, cross-hybridization temperature of 72 °C, 30–90% GC content, no more than six contiguous identical nucleotides and spaces of at least 2 nt between adjacent probes. Primary probes were screened for potential non-specific binding with Bowtie2 (–very-sensitive-local) against mm10 genome sequences. Probes with more than one binding site were filtered out. Thirty qualified probes per gene closest to the transcription start site were selected. Spacers of 3 nt (random sequence) were extended at the 5′ and 3′ ends of the 35-nt probe sets. Two readout sequences (15 nt) separated by a 2-nt spacer (random sequence) were added at both the 5′ and 3′ ends of the probe, respectively, for the potential of performing seqFISH experiments. Universal primer sequences were then attached at the 5′ and 3′ ends. The 5′ primer contains a T7 promoter. The total length of each probe is 147 nt. The oligonucleotide probe pool (6,180 probes) was purchased from Twist Bioscience.

For probe amplification, limited PCR cycles were used to amplify the designated probe sequences from the oligonucleotide complex pool with Kapa HiFi HotStart Polymerase (Roche, KK2502). The amplified PCR products were then purified using a Zymo DNA Clean and Concentrator kit (Zymo Research, D4014) according to the manufacturer’s instructions. The PCR products were used as the template for in vitro transcription (New England Biolabs, E2040S) followed by reverse transcription (Thermo Fisher Scientific, EP7051) with the forward primer. After reverse transcription, the probes were subjected to uracil-specific excision reagent enzyme (New England Biolabs, N5505S) treatment for ~24 h at 37 °C. Probes were then alkaline hydrolyzed with 1 M NaOH at 65 °C for 15 min to degrade the RNA templates, followed by 1 M acetic acid neutralization. Next, to clean up the probes, we used ssDNA/RNA Clean & Concentrator (Zymo Research, D7011) before hybridization.

Chr 2 intron-FISH

Chr2 intron-FISH was performed following a revised seqFISH protocol as described by Shah et al.³². Mouse ES cells were plated on human recombinant laminin 511-coated coverslips (Electron Microscopy Sciences, 72196-25). Cells were then fixed using 4% formaldehyde (Thermo Fisher Scientific, 28908) in 1× PBS diluted in molecular biology-grade water (Corning, 46-000-CM) for 15 min at 20 °C, washed with 1× PBS a few times and incubated in 70% ethanol for about 3 h at room temperature.

The coverslips were then washed twice with 2× SSC. For primary probe hybridization, samples were incubated with primary Chr 2 intron probes for 30 h at 37 °C in 50% Hybridization Buffer (2× SSC, 50% (vol/vol) formamide (Thermo Fisher Scientific, AM9344) and 10% (wt/vol) dextran sulfate (Sigma-Aldrich, D8906) in molecular biology-grade water) and washed in 55% Wash Buffer (2× SSC, 55% (vol/vol) formamide and 0.1% Triton X-100 (Sigma-Aldrich, 93443)) for 30 min at room temperature, followed by washing in 2× SSC. Alexa Flour 647-coupled Imager probes (Integrated DNA Technologies) for the first round of hybridization were incubated for 20 min at 50 nM each at room temperature in 10% EC buffer (10% ethylene carbonate (Sigma-Aldrich, E26258), 2× SSC, 0.1 g ml^–1 dextran sulfate and 0.02 U ml^–1 SUPERase·In RNase Inhibitor (Thermo Fisher Scientific, AM2694)) and washed for 5 min at room temperature in 10% Wash Buffer (2× SSC, 10% (vol/vol) formamide and 0.1% Triton X-100), followed by a 1-min wash in 2× SSC.

Samples were then co-stained with 5 ng ml^–1 DAPI for 15 min and imaged in an antibleaching buffer (50 mM Tris-HCl (pH 8.0), 2× SSC, 3 mM Trolox (Sigma-Aldrich, 238813), 0.8% D-glucose (Sigma-Aldrich, G7528), 100-fold diluted catalase (Sigma-Aldrich, C3155), 0.5 mg ml^–1 glucose oxidase (Sigma-Aldrich, G2133) and 0.02 U ml^–1 SUPERase·In RNase Inhibitor).

Intron-seqFISH data acquisition

Fifteen rounds of imaging were performed with transistor–transistor logic automated imaging and a fluidic system consisting of a Nikon CSU-W1 spinning disk microscope with a CFI Plan Apochromat ×60/1.42-NA oil objective lens and a spatial genomics fluidics pack from Elveflow. After primary probe hybridization, the 40-mm coverslip (number 1.5) with cells was mounted into a closed-top FCS2 chamber (Bioptechs) and subsequently loaded into a custom stage adaptor on the microscope. The flow rate was normally maintained at 100 µl min^–1. For each readout probe hybridization, cells were first washed with the Wash Buffer (2× SCC with DAPI) for 5 min, followed by 20 min of 10% EC buffer injection (10% ethylene carbonate (Sigma-Aldrich, E26258), 2× SSC, 0.1 g ml^–1 dextran sulfate (Sigma-Aldrich, D4911) and 0.02 U ml^–1 SUPERase In RNase Inhibitor (Invitrogen, AM2694); flow rate of 50 µl min^–1) and 10 min of hybridization (flow rate of 0 µl min^–1). Samples were then washed with 10% Wash Buffer (2× SCC, 10% formamide and 0.1% Triton X-100) for 10 min and then with Wash Buffer (2× SCC with DAPI) for 5 min. After these two washes, imaging buffer (50 mM Tris-HCl (pH 8.0), 300 mM NaCl, 2× SSC, 3 mM Trolox (Sigma, 238813), 0.8% d-glucose (Sigma, G7528), 100-fold diluted catalase (Sigma, C3155), 0.5 mg ml^–1 glucose oxidase (Sigma, G2133) and 0.02 U ml^–1 SUPERase·In RNase Inhibitor (Invitrogen, AM2694)) was flowed into the chamber for 10 min, followed by zero flow implementation with the T junction method recommended by Elveflow. Imaging acquisition was performed on a Nikon CSU-W1 equipped with the Uniformizer, five laser lines (405 nm/514 nm/561 nm/594 nm/640 nm) and a Hamamatsu BT fusion camera. The four image stacks were acquired sequentially from the 640-nm channel (Alexa Fluor 647N) to the 561-nm (Cy3), 488-nm (Alexa Fluor 488) and 405-nm (DAPI + blue beads) channels under ultraquiet mode with a fixed framerate (5.1 Hz). A total axial range of 20 µm with z steps of 300 nm were covered. To detect weaker single-molecule signals, four-frame averaging was used for 640-nm, 561-nm and 488-nm channels, and two-frame averaging was used for the DAPI channel. After imaging, cells were washed with probe stripping buffer (2× SCC, 55% formamide and 0.1% Triton X-100) for 10 min before starting the next round of hybridization. The z position of the objective between imaging rounds was maintained by Nikon’s perfect focusing system.

seqFISH cell-level gene decoding and data analysis

Three-dimensional single-molecule localizations were performed using FISH-Quant 2 (refs. ^57,58) and fixed thresholds for each color channel. xyz drift corrections were performed based on localizations of blue beads in the 405-nm channel. Chromatic corrections were calculated based on 3D multicolor bead image stacks acquired using coverslips coated with tetraspeck beads. Localizations from different hybridization rounds and color channels were pooled after applying corresponding drift and chromatic corrections. Gene decoding was performed according to the code book (Supplementary Table 3) with a maximal xyz distance cutoff of 1 pixel between all readout rounds. Three-dimensional nucleus segmentation was performed by using Cellpose 2 (ref. ⁵⁹) with a pretrained specialist’s model optimized for our imaging condition. Decoded genes were parsed into single cells based on gene localizations and 3D masks in the field of view.

The Chr 2 intron-FISH data were screened for number (N) of cells. The number of bursting (N_b) events for each gene was counted over all the cells detected. The number of co-bursting (N_co) events over N cells for each gene pair was counted when a visible and approximated (Euclidean distance ≤ 1 µm) pair of intron-FISH localizations from two channels were identified. The bursting frequency (F) of each gene was computed by

The normalized co-bursting frequency (F_co) for gene pairs (A and B) was computed by normalizing co-bursting frequency (N_co/N) to the product of the bursting frequencies of both genes

F_A and F_B are the bursting frequencies for A and B, respectively.

Airyscan imaging and image analysis

Fluorescence images obtained from smRNA-FISH or Chr 2 intron-FISH were acquired on a Zeiss LSM880 inverted confocal microscope attached to an Airyscan 32 gallium arsenide phosphide-PMT area detector. Before imaging, the beam position was calibrated to center on the 32-detector array. Images were taken under Airyscan super-resolution mode by a Plan Apochromat ×63/1.40-NA oil objective in a lens immersion medium with a refractive index of 1.515. The Airyscan super-resolution technology used a very small pinhole (0.2 AU) at each of its 32 detector elements to increase the signal-to-noise ratio approximately 4- to 8-fold and enable ~1.7-fold improvement of resolution after linear deconvolution in both lateral (x, y) and axial (z) directions. DAPI, mEmerald, Quasar 570 (and JF549) and Quasar 670 signals were illuminated/detected at excitation/emission wavelengths of 405 nm/460 nm, 488 nm/510 nm, 561 nm/594 nm and 633 nm/654 nm, respectively. z-stacks were acquired with a step of 300 nm. After image acquisition, Airyscan images were processed and reconstructed using the provided algorithm from the Zeiss LSM880 platform.

Three-dimensional Airyscan image stacks were processed by using Imaris 7.2.3. As an initial step, we manually inspected the images and removed low-quality images based on the following criteria: (1) non-specific signal outside of the DAPI-stained nuclei, (2) cropped signal at the edge of the images and (3) very faint signal.

To characterize and visualize MED6 and MED1 protein hubs, we used the Surfaces object in Imaris and ran the following algorithms: (1) apply ‘Background Subtraction (Local Contrast)’ mode and set ‘Diameter of Largest Sphere’ as 800 nm, (2) uniformly use a threshold value of 30 for MED6 and 10 for MED1 to segment protein hubs and (3) set the ‘Minimal Number of Voxels’ as 20 for filtering out noise signals. To characterize and visualize smRNA-FISH puncta, we used the Surfaces object from Imaris and and ran the following algorithms: (1) apply ‘Background Subtraction (Local Contrast)’ mode and set ‘Diameter of Largest Sphere’ as 2,000 nm, (2) uniformly use a threshold value of 35 for Quasar 570 and 70 for Quasar 670 to illustrate smRNA-FISH puncta and (3) set the ‘Minimal Number of Voxels’ as 50 for filtering out noise signals. To characterize the colocalization signal between Mediator hubs and smRNA-FISH puncta, we used the ‘Colocalization’ function in Imaris with a common threshold value of 100 to identify the overlapping volumes and reconstruct the 3D isosurfaces.

To measure the 3D distance between puncta of different smRNA-FISH signals, we localized the voxels corresponding to the local maximum of identified RNA-FISH signal using the Imaris ‘Spots’ function module and calculated the Euclidean distance by using the Measurement Points object with a voxel size of 43.6 nm × 43.6 nm × 300 nm. The potential drifts among different imaging channels were estimated by using 100-nm multispectral beads under the same acquisition settings and were considered in the calculation of Euclidean distance between smRNA-FISH foci.

The gene pair smRNA-FISH data were screened for number (N) of cells. The number of bursting (N_b) events for each gene was counted when a visible smRNA-FISH puncta was identified. The number of co-bursting (N_co) events over N cells for each gene pair was counted when a visible and approximated (Euclidean distance of <1 µm) pair of smRNA-FISH puncta from two channels were identified. The bursting frequency (F) of each gene and the normalized co-bursting frequency (F_co) were computed according to Eqs. (5) and (6), respectively.

Three-dimensional ATAC-PALM, Mediator–HaloTag PALM and image analysis

We prepared the reagents for 3D ATAC-PALM experiments as described previously²¹. One day before the experiment, cells were plated on 5-mm coverslips (Warner Instruments, 64-0700) at a confluency of around 70–80% with proper coating. Cells were fixed with 3.7% formaldehyde (Millipore Sigma, F8775-25ML) for 10 min at room temperature. After fixation, cells were washed three times with 1× PBS for 5 min and permeabilized with ATAC lysis buffer (10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl₂ and 0.1% Igepal CA-630) for 10 min at room temperature. After permeabilization, the sample coverslips were washed twice in 1× PBS, and the transposase mixture solution (1× Tagmentation buffer: 10mM Tris-HCl (pH 7.6), 5 mM MgCl₂, 10% dimethylformamide and 100 nM Tn5-PA-JF549) was added to the sample. The coverslips were placed in a humidified chamber and incubated for 30 min at 37 °C. After the transposase reaction, the coverslips were washed three times in 1× PBS containing 0.01% SDS and 50 mM EDTA for 15 min at 55 °C before being mounted onto the lattice light-sheet microscope sample stage for imaging.

Three-dimensional ATAC-PALM data were acquired by lattice light-sheet microscopy at room temperature⁶⁰. The light sheet was generated from the interference of highly parallel beams in a square lattice and dithered to create a uniform excitation sheet. The inner and outer numerical apertures of the excitation sheet were set to be 0.44 and 0.55, respectively. A variable-flow peristaltic pump (Fisher Scientific, 13-876-1) was used to connect a 2-l reservoir with the imaging chamber with 1× PBS circulating through at a constant flow rate. Labeled cells seeded on 5-mm coverslips were placed into the imaging chamber, and each imaging volume took 100~200 image frames, depending on the depth of the field of view. Specifically, spontaneously activated PA-JF549 dye was initially pushed into the fluorescent dark state through repeated photobleaching by scanning the whole imaging volume with a 2-W, 560-nm (or 640-nm) laser (MPB Communications). The samples were then imaged by iteratively photoactivating each plane with very-low-intensity 405-nm light (<0.05-mW power at the rear aperture of the excitation objective and 6 W cm^–2 power at the sample) for 8 ms and by alternatively exciting each plane with a 2-W, 560-nm laser and a 2-W, 640-nm laser at its full power (26-mW power at the rear aperture of the excitation objective and 3,466 W cm^–2 power at the sample) for an exposure time of 20 ms. The specimen was illuminated when laser light went through a custom 0.65-NA excitation objective (Special Optics), and the fluorescence generated within the specimen was collected by a detection objective (CFI Apo LWD water immersion ×25/1.1-NA, Nikon), filtered through a 440/521/607/700-nm BrightLine quad-band bandpass filter (Semrock) and N-BK7 Mounted Plano-Convex Round cylindrical lens (f = 1,000 mm, Ø 1′, Thorlabs, LJ1516RM) and eventually recorded by an ORCA-Flash 4.0 sCMOS camera (Hamamatsu, C13440-20CU). The cells were imaged under sample scanning mode and the dithered light sheet at a step size of 500 nm, thereby capturing a volume of ~25 µm × 51 µm × (27~54) µm, considering a 32.8° angle between the excitation direction and the stage moving plane.

To precisely analyze the 3D ATAC-PALM data, we embedded nano-gold fiducials within the coverslips for drift correction as previously described²¹. ATAC-PALM images were taken to construct a 3D volume when the sample was moving along the ‘s’ axis. Individual volumes per acquisition were automatically stored as TIFF stacks, which were then analyzed by in-house-developed scripts in MATLAB. The cylindrical lens introduced astigmatism in the detection path and recorded each isolated single molecule with its ellipticity, thereby encoding the 3D position of each molecule relative to the microscope focal plane. All processing was performed by converting all dimensions to units of x–y pixels, which were 100 nm × 100 nm after transformation due to the magnification of the detection objective and tube lens. We estimated the localization precision by calculating the standard deviation of all the localization coordinates (x, y and z) after the nano-gold fiducial correction. The localization precision is 26 ± 3 nm and 53 ± 5 nm for xy and z, respectively.

Three-dimensional pair cross-correlation function

The 3D pair cross-correlation function c₀(r) between localizations of molecule A and those of molecule B can be formulated as

M is the total number of localizations for molecule A, and N is the total number of localizations for molecule B. Δr = 50 nm is the binning width used in the analysis. The Dirac Delta function is defined by

where r_ij represents the pairwise Euclidean distance between localization points i and j. The normalized 3D pair autocorrelation function \(C(r)\) was calculated by

\({c}_{r}(r)\) refers to the pair cross-correlation function calculated from uniform distributions with the same localization density in the same volume as real data used.

Single-molecule imaging

Single-molecule imaging experiments were performed as previously described^42,49 on a Nikon Eclipse TiE motorized inverted microscope equipped with a ×100/1.49-NA oil immersion objective lens (Nikon), four laser lines (405/488/561/642 nm), an automatic TIRF illuminator, a perfect focusing system, a tri-cam splitter, three EMCCDs (iXon Ultra 897, Andor) and Tokai Hit environmental control (humidity, 37 °C, 5% CO₂). Proper emission filters (Semrock) were switched in front of the cameras for GFP and JF549 emission, and a band mirror (405/488/561/633-nm BrightLine quad-band bandpass filter, Semrock) was used to reflect the laser into the objective.

To perform single-molecule imaging of transcription factors and cofactors, cells were seeded on 25-mm number 1.5 coverglass precleaned with potassium hydroxide and ethanol and coated with human recombinant laminin 511 according to manufacturer’s instructions. Live-cell imaging experiments were conducted by culturing mouse ES cells in imaging medium composed of FluoroBrite DMEM (Thermo Fisher Scientific, A1896701), 10% fetal bovine serum, 1× GlutaMAX, 1× non-essential amino acids, 0.1 mM β-mercaptoethanol and 1,000 U ml^–1 leukemia inhibitory factor. The TIRF illuminator was adjusted to deliver a highly inclined laser beam to the cover glass with the incident angle smaller than the critical angle. Oblique illumination (HILO) has much less out-of-focus excitation than regular epi-illumination. Transcription factors and cofactors linked to HaloTag were labeled with 5 nM HaloTag ligand-JF549 for 15 min and imaged using a 561-nm laser with an excitation intensity of ~50 W cm^–2. To minimize drift, the imaging experiments were performed in an ultraclear room with a precise temperature control system. The environment control chamber for cell culturing was thermoequilibrated. The imaging system was calibrated with beads to confirm a minimal drift during imaging (x–y drift < 100 nm h^–1).

For fast-molecule tracking and throughout the experiments, we used an imaging acquisition time of 10 ms and took 5,000 continuous frames per imaging view after photobleaching of saturatedly labeled single molecules at the beginning. About 10~15 views were imaged for each labeled transcription factor/cofactor under normal or cohesin-depleted conditions. By contrast, for evaluating the RoC, which captures the motility of stably bound molecules, we used a longer acquisition time of 50 ms and took 5,000 continuous frames per imaging view.

Single-particle tracking analysis

Each imaging view was recorded as a TIFF stack, and single molecules were tracked using SLIMfast, a custom-written MATLAB implementation of the MTT algorithm⁶¹. Frame-to-frame motions are defined by the distance between consecutive positions of the particle and can be potentially related to (1) Brownian (random walk) or confined motions of the molecule or (2) potential artifactual effects such as imperceptible movements of the nucleus. To filter out the effects resulting from point 1, it is thus necessary to define a maximal expected diffusion coefficient (D_Max), which defines the maximal distance (d_m) between two consecutive frames for a particle to be considered as the same object. As in the previous publication⁴², a cutoff was set to 3 d_m to ensure a 99% confidence level, and D_Max was set as 1 μm² s⁻¹.

For each imaging view, SLIMfast generated a .txt output file consisting of a series of successive x/y coordinates and times of detection corresponding to the displacement of each individual molecule. The output SLIMfast.txt files included the following information: x/y coordinates (two-dimensional coordinates of the molecule in micrometers), trajectory index (ID number of the trajectory) and frame number (the index of the frame on which each single molecule was detected). These track files were used as the inputs to perform two-state (for histone subunits) or three-state (for transcriptional regulators) kinetic fitting by using Spot-On⁴⁵ software to compute the biophysical parameter values of single particles.

The RoC represents the circle best encompassing the motion track rather than encompassing it strictly. Thus, the measurement of the RoC is largely independent of the track duration. Tracks with lengths of <5 frames were discarded in the preprocessing step. To quantify the RoC, the mean square displacement (MSD) curves of each track were fitted using the nonlinear least-squares approach in MATLAB with a circle confined diffusion model⁶² as illustrated in the following equation:

The fitting provided values for RoC, the diffusion coefficient at short timescales D and a constant offset value due to the localization precision limit, which is inherent to all the localization-based microscopy methods. To discard fitting errors related to artifacts such as erroneously connected jumps, we have discarded the trajectories with squared norm of the residual higher than 10⁻⁵ and a RoC higher than 500 nm.

Statistics and reproducibility

Unless specified, data are presented as mean ± s.d. We normally applied two-sided Student’s t-tests for measurements of technical replicates among different conditions but applied a non-parametric two-sided Wilcoxon rank-sum test for data that clearly do not follow normal distribution. No statistical method was used to predetermine sample size, but our sample sizes are similar to those reported in previous publications^18,25. No data were excluded from the analyses. The experiments were not randomized, and the investigators were not blinded to allocation during experiments and outcome assessment. For results shown for representative experiments, each measurement was repeated three times independently with similar results.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.