Suv39h-catalyzed H3K9me3 is critical for euchromatic genome organization and the maintenance of gene transcription

Mice

Suv39h1 and Suv39h2 null mice (a generous gift from Thomas Jenuwein, Max Planck Institute) were backcrossed to the C57BL/6J strain on which background the Suv39 double-deficient (DKO) genotype was found to be embryonic-lethal (Keenan et al. 2020). We therefore used 8- to 12-wk-old chimeric mice for this study as previously described (Keenan et al. 2020). Briefly, chimeric mice were generated using Suv39h1 and Suv39h2 double-deficient (DKO) CD45.2⁺ fetal liver or bone marrow donor cells, which were injected into the tail vein of lethally irradiated (2 × 550 Gy) congenic (CD45.1⁺) recipient mice. These DKO cells are compared with littermate control cells (designated CON) of the Suv39h1^+/ySuv39h2^+/− genotype, as we used a mating strategy unable to produce wild-type cells to increase the frequency of the submendelian DKO genotype. Our previous studies have shown no defects in the development, function, or heterochromatin structure of Suv39h2^−/− mice, making the Suv39h2^+/− a reasonable surrogate for a true wild type. All mice were maintained at The Walter and Eliza Hall Institute animal facility under specific-pathogen-free conditions. All animal experiments were approved in advance by the Walter and Eliza Hall Institute animal ethics committee and conducted in accordance with published guidelines.

Flow cytometric sorting

Single-cell suspensions were generated from mouse thymus by mechanical homogenization following by red blood cell lysis (156 mM NH₄Cl, 11.9 mM NaHCO₃, 97 μM EDTA). Fluorochrome-conjugated antibodies against the following mouse antigens were then used for sorting by flow cytometry: CD45.2-FITC (clone 104) and CD19-PECy7 from BD Pharmingen and CD4-PE (clone GK1.5) and CD8α-APC (clone 53.6.7), which were generated internally. Surface staining was performed for 30 min at 4°C. Live (propidium iodide [Thermo Fisher Scientific] negative), CD45.2⁺CD19⁻CD4⁺CD8⁺ cells were sorted using a BD FACSAria cell sorter (Supplemental Fig. S1).

RNA-seq

RNA was extracted using an RNeasy plus mini kit (Qiagen) and quantified in a TapeStation 2200 using RNA ScreenTape (Agilent). Libraries for sequencing were prepared using the TruSeq RNA sample-preparation kit (Illumina) from 500 ng RNA, amplified with KAPA HiFi HotStart ReadyMix (Kapa Biosystems), size-selected to 200–400 bp, and cleaned up with AMPure XP magnetic beads (Beckman). Final libraries were quantified with TapeStation 2200 using D1000 ScreenTape for sequencing on the Illumina NextSeq 500 platform to produce 81-bp paired-end reads. Around 35 million read pairs were generated per sample.

All read pairs were aligned to the mouse genome, build mm10, using align from the Rsubread package (Liao et al. 2019) v1.24.0 in R (R Core Team 2023). Over 91% of all read pairs were successfully mapped for each sample. The number of read-pair-overlapping genes were summarized into counts using Rsubreads featureCounts (Liao et al. 2014). Of those successfully mapped, an average of 89% of these read pairs were assigned to a gene. Genes were identified using GENCODE annotation to the mm10 genome, version M22. Differential expression analyses were then undertaken using the limma (Ritchie et al. 2015) and edgeR (Robinson et al. 2010) software packages, versions 3.42.2 and 3.28.1, respectively.

Before analysis, all genes with no symbol, non-protein-coding, and Riken genes were removed. Expression-based filtering was also performed. All genes were required to achieve a count per million reads (CPM) greater than 0.4 in at least three samples to be retained. Following filtering 11,036 genes remained. Compositional differences between samples were then normalized using the trimmed mean of M-values (TMM) method (Robinson and Oshlack 2010). All counts were transformed to log₂-CPM with associated precision weights using voom (Law et al. 2014). Differential expression between the knockout and wild-type samples was evaluated using linear models and robust empirical Bayes moderated t-statistics relative to a fold-change threshold of 1.2 (McCarthy and Smyth 2009; Phipson et al. 2016). P-values were adjusted to control the false-discovery rate (FDR) <5% using the Benjamini and Hochberg method.

Analyses of the GO terms (The Gene Ontology Consortium 2019) was performed using limma's goana function. Gene set testing of the aging signature and corresponding barcode plots were generated using limma's roast and barcodeplot functions. The aging signature was evaluated using both directional and nondirectional roast tests with 9999 rotations. The heatmap was generated using the pheatmap package available through CRAN.

Lamin B1 ChIP-seq

Cells (>5 M) for Lamin B1 ChIP-seq were fixed in 1% formaldehyde for 20 min at room temperature with mixing and then quenched by the addition of glycine to a final concentration of 125 mM. Fixed cells were washed twice with ice-cold PBS and then the cell membrane was lysed by incubation in ChIP buffer (150 mM NaCl, 50 mM Tris, 5 mM EDTA, 0.5% (v/v) NP-40, 1% (v/v) Triton X-100) containing 1× protease inhibitors (Roche) on ice for at least 10 min before Dounce homogenization. Chromatin within intact nuclei was then digested for 5 min with micrococcal nuclease (NEB M0247, 500 U/1 M cells) before nuclei lysis by Covaris sonication. Nuclear lysate was then cleared by centrifugation (>15,000g for 1 min); a small aliquot was taken as a whole-cell extract (WCE) control; and Lamin B1 antibody (Abcam ab16048, 3 μg) was added to the remaining cleared lysate and incubated overnight at 4°C with rotation. Immunoprecipitation was performed using Protein G Dynabeads (Thermo Fisher Scientific), and then chromatin was eluted from beads using elution buffer (1% (v/v) SDS, 0.1 M NaHCO₃). ChIP and WCE samples were then adjusted to 200 mM NaCl, digested with RNase A (20 μg/sample) overnight at 65°C, and then digested with Proteinase K (20 μg/sample) for a further 1 h at 65°C. DNA was then extracted using a Zymo ChIP DNA clean and concentrator kit (D5205), and sequencing libraries were prepared using an Illumina TruSeq DNA sample prep kit using half volumes of all reagents. Libraries were sequenced on an Illumina NextSeq 500 to produce 81-base-length paired-end reads. Independent biological duplicates were prepared for each condition and sequenced to a depth of about 50 million reads.

Lamin B1 ChIP-seq processing and peak calling

Biological replicate FASTQ files were combined for ChIP and WCE samples. FASTQ files were aligned to the mouse genome (mm10) with Rsubread package v1.28 (Liao et al. 2019). BAM files were sorted with Rsamtools, and duplicate reads were removed with MarkDuplicates (Picard suite v1.117) (https://broadinstitute.github.io/picard/). WCE samples were down-sampled to the size of the corresponding ChIP libraries (about 53 million reads) using DownsampleSam from Picard tools with strategy chained and accuracy 0.001. For viewing and plotting, bedGraph files of the log of the ratio of the ChIP coverage to the WCE coverage were created with bamCoverage from deepTools (Ramírez et al. 2016) and R (R Core Team 2023). Tracks were visualized in IGV (Robinson et al. 2011). Enriched domain detector (EDD) (Lund et al. 2014) was used to call LADs on the control and DKO samples separately using the WCE as input for the corresponding ChIP. EDD was run 10 times without specifying the bin size and gap penalty. The mean of the estimated bin size and gap penalty from the 10 trials was used to call LADs (wild-type: gap penalty 4.79 and bin size 6 kbp; KO: gap penalty 9.1 and bin size 8 kbp). LADs called on the Y and scaffold chromosomes were excluded from subsequent analyses. LADs were also defined with methods based on circular binary segmentation (CBS) and LADetector (Harr et al. 2015; Luperchio et al. 2017). With the down-sampled BAM, we used featureCounts to create the log₂ ratio of the ChIP to WCE in 10-kbp bins. Bins overlapping blacklisted regions or gaps in the genome (e.g., telomere, centromere) were excluded. We partitioned this ratio using the DNAcopy package in R (Seshan and Olshen 2023) with alpha = 0.001. LADs were classified as regions > 100 kb in size of positive signal (seg.mean), allowing for smaller regions of negative signal < 10 kb in size on autosomes and Chr X; iLADs were classified the same but with negative signal. After defining LADs in the CON sample, we calculated the median seg.mean in each LAD for DKO and CON. LADs with a median seg.mean decreasing by more than 0.1 were classified as “disrupted,” and all other LADs were classified as “nondisrupted.”

Histone mark ChIP-seq

ChIP-seq was performed as previously described (Lee et al. 2006) with amendments. Briefly, 2 × 10⁶–5 × 10⁶ cells per IP were cross-linked by adding 1/10 volume of fresh 11% formaldehyde solution (50 mM HEPES-KOH at pH 7.4, 100 mM NaCl, 1 mM EDTA, 1 mM EGTA, 11% formaldehyde) to a 1 × 10⁶ cells/mL cell suspension for 20 min at room temperature with mixing. Then, 1/12.5 volumes of 2.5 M glycine was added to quench formaldehyde, and cells were washed twice in ice-cold PBS. Cross-linked cells were lysed in lysis buffer 1 (50 mM HEPES-KOH at pH 7.5, 140 mM NaCl, 1 mM EDTA, 10% [v/v] glycerol, 0.5% [v/v] IGEPAL CA-630, 0.25% [v/v] Triton X-100, 1× Roche cOmplete, EDTA-free protease inhibitor cocktail) for 10 min at 4°C with rocking, pelleted at 1350g for 5 min at 4°C and then with lysis buffer 2 (10 mM Tris-HCl at pH 8.0, 200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 1× Roche cOmplete, EDTA-free protease inhibitor cocktail) for 10 min at room temperature with rocking, pelleted again at 1350g for 5 min at 4°C, then resuspended in lysis buffer 3 (10 mM Tris-HCl at pH 8.0, 100 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 0.1% [w/v] sodiam deoxycholate, 0.5% [v/v] N-lauroylsarcosine, 1× Roche cOmplete, EDTA-free protease inhibitor cocktail), and transferred to Covaris miniTUBEs. Cell lysates were sonicated using a Covaris S220 focused-ultrasonicator with the following settings to shear DNA to a size of 300–500 bp suitable for high-throughput sequencing: fill level 10, duty cycle 10, PIP 175, cycles/burst 200 for 60 sec. Then 1/10 volume of 10% (v/v) Triton X-100 was added to the sonicated lysate, and cellular debris was pelleted at 20,000g for 10 min at 4°C. Immunoprecipitation was then performed overnight at 4°C with rotation using ProteinG Dynabeads (50 μL per IP) precoupled with antibody (5 μL per IP) in 0.5% BSA/PBS for >5 h. Antibodies used were H3K9me3 Rabbit pAb (Abcam 8898 lot GR3217836-1), H3K27ac (Abcam ab4729 lot 3357415-3), H3K27me3 (Cell Signaling Technology 9733S lot 16), H3K4me3 (EMD Millipore 07-472, lot 3746), and rabbit IgG pAb (Abcam ab46540). Antibody specificity was empirically confirmed by MODified histone peptide array (Active Motif) (Supplemental Fig. S3). Following immunoprecipitation, beads were collected using a magnet and washed five times with wash buffer (50 mM HEPES-KOH at pKa 7.55, 500 mM LiCl, 1 mM EDTA, 1% (v/v) IGEPAL CA-630, 0.7% (w/v) sodium deoxycholate) and then one time with Tris-EDTA-NaCl buffer (10 mM Tris-HCl at pH 8.0, 1 mM EDTA, 50 mM NaCl). DNA was then eluted from beads by incubation in elution buffer (50 mM Tris-HCl at pH 8.0, 10 mM EDTA, 1% (w/v) SDS) for 30 min at 65°C with shaking. Bead-free supernatant was reverse crosslinked by incubation for 6–15 h at 65°C, diluted with equal volume of TE buffer (10 mM Tris-HCl at pH 8.0, 1 mM EDTA), and then digested sequentially with RNase A (0.2 mg/mL) for 2 h at 37°C and Proteinase K (0.2 μg/mL) for 2 h at 55°C. Final immunoprecipitated DNA was then purified using a Zymo ChIP DNA clean and concentrator kit and prepared for sequencing using an Illumina TruSeq ½ reaction library preparation kit. Final libraries were quantified using an Agilent TapeStation using a D1000 screentape. Libraries for H3K9me3 were sequenced on an Illumina NextSeq 500, generating 76-bp paired-end reads, and all other libraries were sequenced on an Illumina NextSeq 2000, generating 66-bp paired-end reads.

Histone mark ChIP-seq processing and peak calling

FASTQ files of H3K27ac, H3K27me3, H3K4me3, H3K9me3, and IgG from CON and H3K27ac, H3K27me3, H3K4me3, and IgG from DKO-aligned were with Rsubread with the mm10 reference genome. Duplicate reads were removed with Picard tools. The HOMER pipeline was used to call peaks (Heinz et al. 2010). Tag directories were created for each library with the genome specified. For ChIP-seq samples, the findPeaks function was used to find peaks with style as histone, size as auto, and the corresponding IgG sample as input.

ATAC-seq

ATAC-seq was performed using the Omni-ATAC method (Corces et al. 2017). Briefly, 50,000 cells were pelleted at 500g for 5 min at 4°C and then lysed on ice for 3 min in 50 μL ATAC-RSB buffer (10 mM Tris-HCl at pH 7.4, 10 mM NaCl, 3 mM MgCl₂) also containing 0.1% (v/v) IGEPAL CA-630, 0.1% (v/v) Tween-20, and 0.01% (w/v) digitonin. After lysis, 1 mL of ATAC-RSB containing only 0.1% (v/v) Tween-20 but no IGEPAL CA-630 or digitonin was added, and nuclei were pelleted at 500g for 10 min at 4°C. All supernatant was aspirated and discarded, and each pellet of nuclei was resuspended in 50 μL of transposition mix (25 μL 2× TD buffer, 2.5 μL Tn5 transposase [100 nM final], 16.5 μL PBS, 0.1 μL 5% digitonin, 0.5 μL 10% Tween-20, 5 μL H₂O) and incubated for 30 min exactly at 37°C in a thermomixer with 1000 RPM mixing. DNA was then purified using a Zymo DNA clean and concentrator-5 kit to generate 20 μL of product, which was preamplified for five cycles using NEBNext 2X master mix (E7649A) and indexed primers as previously described (Buenrostro et al. 2013). Five microliters of preamplified DNA was then used in a 15 μL qPCR reaction using Promega GoTaq qPCR master mix to determine the required number of additional cycles to yield minimally amplified DNA for optimal library diversity for massively parallel sequencing. Amplified DNA was then purified using NucleoMag NGS clean-up and size select beads (Macherey-Nagel) using a 0.5× ratio of beads to DNA to initially clear large DNA fragments from the library and then a 1.5× ratio of beads to supernatant to adsorb DNA fragments of appropriate size for sequencing. DNA-bound beads were then washed twice with 80% EtOH, and size-selected DNA was eluted in nuclease-free water. Final libraries were quantified using an Agilent TapeStation using a D5000 screentape and sequenced using an Illumina NextSeq 2000 to generate about 50 million paired-end 165-bp reads per library.

ATAC-seq data were analyzed as previously described (Milevskiy et al. 2023). Briefly, Nextera adapter sequences were trimmed with trim_galore (https://github.com/FelixKrueger/TrimGalore), and then reads were mapped to the mm10 genome with Bowtie 2 with arguments ‐‐maxins 2000 ‐‐minins 38 ‐‐very-sensitive -p 48 -k 10. Duplicates were marked with sambamba (Tarasov et al. 2015). SAMtools (Danecek et al. 2021) was used to remove reads on Chr M and reads with FLAG 1804 (including duplicates and multi-mapped reads). BAM files of the biological replicates were merged with SAMtools. The HOMER pipeline was used to call peaks with the findPeaks function with style as factor and size as auto.

DNA fluorescence in situ hybridization

DNA fluorescence in situ hybridization (DNA-FISH) was performed on control and Suv39DKO primary DP thymocytes as previously described (Tapia Del Fierro et al. 2023) with 40- to 48-h hybridization time. BAC DNA (BACPAC resources) was isolated using a large-construct kit (Qiagen) and converted to fluorescently labeled DNA-FISH probes by a nick translation reaction (Abbott Laboratories). BACs used in this study are RP24-263I20 (Chr 12: 19,359,827–19,451,909), RP23-277F6 (Chr 2: 31,162,711–31,356,593), RP24-284I15 (Chr 1: 99,904,338–100,058,959), RP23-305P17 (Chr 6: 30,032,262–30,187,634), and RP24-391L20 (Chr 3: 72,874,595–73,009,096). Stained cells were visualized on an LSM 880 (Zeiss) microscope with Airyscan processing.

Microscopy data were analyzed automatically with a custom written macro in FIJI (Schindelin et al. 2012). Briefly, nuclei were segmented from the background with a 3D autothreshold and labeled on a scaled-down version of the DAPI channel. FISH foci were detected and assigned to a location within the nuclei via the 3D maxima detector and the EVF calculator in the 3D ImageJ suite (Ollion et al. 2013). The radial position was computed as the fractional distance of all detected foci between the nuclear periphery and center, with a maximum of two foci per nucleus (based on intensity) used in analysis. The P-value between groups was determined with a Wilcoxon test.

Hi-C

In situ Hi-C was performed on fixed cells as previously described (Johanson et al. 2018). Libraries were sequenced on an Illumina NextSeq 500 to produce 81-base-length paired-end reads. Independent biological duplicates were prepared for each condition and sequenced to a depth of about 200 million reads (Supplemental Table 1).

Hi-C data processing

Data preprocessing and analysis were performed as previously described with changes in the parameters (Johanson et al. 2018). Samples were aligned to the mm10 genome using the diffHic package v1.14.0 (Lun and Smyth 2015), which uses cutadapt v0.9.5 (Martin 2011) and Bowtie 2 v2.2.5 (Langmead and Salzberg 2012) for alignment and Picard suite v1.117 (https://broadinstitute.github.io/picard/) for fixing mate pair information and marking duplicate reads. Dangling ends and self-circling artifacts were identified and removed if the pairs of inward-facing or outward-facing reads were separated by <1000 bp (inward-facing) or 2000 bp (outward-facing). Read pairs with fragment sizes >1000 bp were removed. The proportion of reads removed through each processing step is tabulated in Supplemental Table 1.

The HOMER Hi-C pipeline (Heinz et al. 2010, 2018) was also used for Hi-C analysis. After processing with the diffHic pipeline, libraries in HDF5 format were converted to the Hi-C summary format with R. Then input-tag directions were created for each library with the makeTagDirectory function, with the genome (mm10) and restriction-enzyme-recognition site (GATC) specified. Summed biological-replicate tag directories for each cell type were also created.

To confirm the reproducibility of the libraries, the HiCRep R package was used to quantify the similarity between all libraries with the stratum adjusted correlation coefficient (SCC) (Yang et al. 2017). For every combination of libraries, the raw contact matrices (50-kbp resolution) of individual chromosomes for each replicate were used to compute the SCC with a smoothing parameter of three and a maximum distance considered at 5 Mbp. For each pairwise comparison, a median SCC was calculated across all chromosomes (Supplemental Fig. S7).

Plaid plots were constructed using the contact matrices and the plotHic function from the Sushi R package (Phanstiel et al. 2014). The inferno color palette from the viridis package (Garnier et al. 2024; https://sjmgarnier.github.io/viridis/) was used. The range of color intensities in each plot was scaled according to the library size of the sample to facilitate comparisons between plots from different samples. DI arcs were plotted with the plotBedpe function of the Sushi package. The z-score shown on the vertical access was calculated as −log₁₀ (P-value).

Detection of A/B compartments

A/B compartments were identified at a resolution of 100 kbp with the method outlined by Lieberman-Aiden et al. (2009) using the HOMER Hi-C pipeline (Lin et al. 2012). With the summed biological replicate tag directories, the runHiCpca.pl function was used on each library with −res 100,000 and the genome (mm10) specified to perform principal component analysis to identify compartments. To identify changes in A/B compartments between libraries, the getHiCcorrDiff.pl function was used to directly calculate the difference in correlation profiles. A region was considered to flip compartment if the difference in correlation profile was <−0.1 and was determined from control cell A/B compartments. Chr Y was excluded from the analysis.

Detecting DIs

DIs between the different libraries were detected using the diffHic package v1.16.0 (Lun and Smyth 2015). Read pairs were counted into 50-kbp bin pairs for all chromosomes. Bins were discarded if they had a count of fewer than 10, contained blacklisted genomic regions as defined by ENCODE for mm10 (The ENCODE Project Consortium 2012), or were within a centromeric or telomeric region. Filtering of bin pairs was performed using the filterDirect function, in which bin pairs were only retained if they had average interaction intensities more than threefold higher than the background ligation frequency. The ligation frequency was estimated from the interchromosomal bin pairs from a 1-Mbp bin pair count matrix. Diagonal bin pairs were also removed. The counts were normalized between libraries using a loess-based approach with bin pairs <100 kbp from the diagonal normalized separately from other bin pairs. Tests for DIs were performed using the QL framework (Lund et al. 2012; Chen et al. 2016) of the edgeR package v3.26.5 (Robinson et al. 2010). A design matrix was constructed using a one-way layout that specified the group to which each library belonged. A mean-dependent trend was fitted to the negative binomial dispersions with the estimateDisp function. A generalized linear model (GLM) was fitted to the counts for each bin pair (McCarthy et al. 2012), and the QL dispersion was estimated from the GLM deviance with the glmQLFit function. The QL dispersions were then squeezed toward a second mean-dependent trend, using a robust empirical Bayes strategy (Phipson et al. 2016). A P-value was computed against the null hypothesis for each bin pair using the QL F-test. P-values were adjusted for multiple testing using the Benjamini–Hochberg method. A DI was defined as a bin pair with a FDR <5%. DIs adjacent in the interaction space were aggregated into clusters using the diClusters function to produce clustered DIs. DIs were merged into a cluster if they overlapped in the interaction space. The significance threshold for each bin pair was defined such that the cluster-level FDR was controlled at 5%. Cluster statistics were computed using the csaw package v1.18.0 (Lun and Smyth 2016). Overlaps between unclustered bin pairs and genomic intervals were performed using the InteractionSet package (Lun et al. 2016). Differential interactivity for DE genes was calculated from 50-kbp interactions with P-value < 0.05 before clustering and multiple testing. For each DE gene, the mean logFC from overlapping interactions was calculated.

Detecting TAD boundaries

TAD breakpoints were detected in each Hi-C library with TADbit v0.2.0.5 (Serra et al. 2017). Read pairs were counted into 50-kbp bin pairs (with bin boundaries rounded up to the nearest MboI restriction site). The TADbit tool find_tad was run on the raw counts specifying normalized = false. Only TAD boundaries with a score of seven or higher were included in the results.

Detecting differential TADs

Differential TAD boundaries (DTBs) between groups were detected with the diffHic and edgeR packages (Lun and Smyth 2015) using the approach described in chapter 8 of the diffHic user's guide. This approach adapts the statistical strategy recently described for differential methylation (Chen et al. 2017) to identify DTBs. The strength of a putative TAD boundary was assessed based on the upstream versus downstream intensity contrast at genomic loci, defined as the ratio of upstream to downstream Hi-C reads anchored at that genomic region. edgeR was used to test whether the ratio at the loci significantly increased or decreased in absolute size between groups. This method directly assesses differential boundary strength relative to biological variation without needing to make TADs calls in individual samples. Upstream and downstream read counts were determined for 200-kbp genomic regions and for 2 Mbp up- and downstream. Low-abundance regions with average log₂-counts per million below two were removed along with regions in blacklisted regions, telomeres, or centromeres. Tests for DTBs were performed using the QL framework of edgeR. Regions with FDR below 0.05 were considered to be DTB.

Detecting differentially interacting promoters

Differentially interacting promoters (DIPs) were detected across all libraries with the diffHic package (Lun and Smyth 2015) and the method described previously (Chan et al. 2021) with alterations in parameters. Gene TSSs were defined with the mouse GENCODE gene set annotation (v. M24). Promoter regions (upstream 10 kbp and downstream 5 kbp) were defined with the promoter function from the GenomicFeatures package v1.36.2. Interactions between the promoters and the genome in 15-kbp bins were counted with the connectCounts function. Interchromosomal interactions or interactions in blacklisted regions or with an anchor >20 kbp or in the diagonal were excluded. A loess curve (span 0.1) was fitted to the average log counts per million (logCPM) for all libraries as a function of distance ^ 0.25 (loessFit function from the limma package v3.40.2) (Ritchie et al. 2015). An interaction was required to have an abundance larger than the fitted curved plus two times the mean of the absolute values of the residuals from the loess fit. Interactions for each promoter were then aggregated to produce a count matrix. Low-abundance promoters were filtered using edgeR's filterByExpr function. Only promoters for protein-coding genes were retained. Obsolete Entrez gene indentifications were removed, as were mitochondrial, Riken, and olfactory genes. The counts were normalized between libraries using a loess-based approach with the normOffsets function from the csaw package v1.18.0 (Supplemental Fig. S7).

DIPs were detected with the QL framework (Lund et al. 2012; Chen et al. 2016) of the edgeR package. A design matrix was constructed with a one-way layout that specified the cell type. Using the promoter counts, the estimateDisp function was used to maximize the negative binomial likelihood to estimate the common dispersion across all promoters with trend = none and robust = TRUE (Chen et al. 2014). A GLM was fitted to the counts (McCarthy et al. 2012), and the QL dispersion was estimated from the GLM deviance with the glmQLFit function with robust = TRUE and trend.method = none. The QL dispersions were then squeezed toward a second mean-dependent trend, with a robust empirical Bayes strategy (Phipson et al. 2016) to share information between genes. A P-value was calculated for each promoter using a moderated t-test with glmQLFTest. The Benjamini–Hochberg method was used to control the FDR <5%. Heatmaps of the filtered and normalized logCPM value were plotted with the coolmap function from the limma package. Gene set enrichment was tested using limma's fry function and visualized using limma's barcode enrichment plot.

Coverage plots and coverage of genomic features

For viewing and plotting, bedGraph files were created for ChIP-seq samples with bamCoverage from deepTools with arguments: normalizeUsing RPGC, bs 10, -e, and effectiveGenomeSize 2652783500. For bedGraph files from ATAC-seq samples, the arguments used were as follows: normalizeUsing RPGC, bs 10, and effectiveGenomeSize 2652783500. Coverage plots were created with computerMatrix and plotHeatmap from deepTools. For genes, the arguments used were as follows: scale-regions -S, -a 3000 -b 3000, -bs 10 –, and regionBodyLength 8000. For LAD boundaries with ChIP, the arguments used were as follows: reference-point -S, -a 500000 -b 500000, -bs 10000, and ‐‐averageTypeSummaryPlot median. For LAD boundaries with RNA-seq the argument -bs was changed to 20,000. Blacklisted regions were excluded.

ChIP-seq and ATAC-seq read coverage (log₂ of the read per kilobase per million) of different genomic features (gene promoters, LADs, iLADs, and TAD boundaries) was determined with featureCounts from Rsubread with allowMultiOverlap = TRUE and the rpkm function from R package edgeR with log = TRUE. The read coverage of genes was determined from regions that are 2 kbp upstream of and 500 kbp downstream from the TSS. Nondifferentially expressed genes survived filtering in the DE analysis and adjusted P-value > 0.05. For LADs/iLADs, DIs and differential TAD boundary regions the genome regions excluded blacklisted regions. The P-value between groups was determined with a Wilcoxon test.

Overlaps between genomic features

Overlaps between genomic regions were identified with the overlapsAny function of the IRanges package v2.20.2 (Lawrence et al. 2013). Distances between nearest genomic features was calculated with the distanceToNearest function from the GenomicRanges package (Lawrence et al. 2013).

Motif enrichment

Motif enrichment was performed with the HOMER pipeline on genomic regions with the findMotifsGenome.pl function and specifying the genome as mm10 and size given and for genes with the findMotifs.pl function with mouse specified as the genome. For DIs, enrichment was performed on the unique anchor regions of unclustered DIs relative to a background of the same number of not-significant DI anchors (FDR > 0.05).

Identification of TAD cliques and modeling with Chrom3D

The Chrom3D pipeline was used to analyze the clique structure and model the 3D chromosome structure of the wild-type and DKO cells (Paulsen et al. 2017, 2018) with some modifications. For creating the input GTrack file we used the following: intrachromosomal interactions at 100 kbp, interchromosomal interactions at 1 Mbp, the LADs as determined by EDD in the wild-type and DKO libraries separately, and a diploid cell. Beads were created from TADbit segmentation at 50 kbp of summed biological replicates of the wild-type libraries excluding the Y Chromosome with the boundary rounded down to the nearest 50 kbp. As there were no overlapping TADs, we did not merge beads. For calculating the significance of the intrachromosomal interactions, we used fdr_bh with a cutoff of 1.5 and threshold of 0.05; for interchromosomal interactions, a cutoff of two and threshold of 0.05. TAD cliques were identified as previously described (Paulsen et al. 2019). For each sample type, we simulated 100 independent Chrom3D models, each with a separate seed value. We ran each simulation for 5 million iterations with a 5-μm radius, the nucleus flag, and 0.3 scale. Bead measurements were calculated as the median radius of beads of each feature type for a given model. Distributions of these measurements were compared by Welch's unequal variances t-test.