Accessible Region Conformation Capture (ARC-C) gives high-resolution insights into genome architecture and regulation

Worm strains and culture

C. elegans strains were maintained at 20°C as previously described (Brenner 1974). The following strains were used: Bristol N2 (wild type), YJ55 blmp-1(tm548) (Huang et al. 2014), and GW638 met-2(n4256) set-25(n5021) (Towbin et al. 2012).

Worm growth

Strains for ARC-C and ChIP-seq were grown in liquid culture at 20°C using standard S-basal medium with HB101 bacteria. Animals were first grown to the adult stage, bleached to obtain embryos, and the embryos hatched without food in M9 buffer for 24 h at 20°C to obtain synchronized starved L1 larvae. L1 larvae were grown in a further liquid culture at 20°C then harvested at the L3 stage. Worms were collected, washed in M9 buffer, floated on sucrose, washed again in M9, then frozen into small pellets by dripping worm slurry into liquid nitrogen, which was stored at −80°C until use.

ChIP-seq

JMJD-2, SCC-1, ZFP-1, BLMP-1, TOP-2, and SIN-3 chromatin immunoprecipitations in L3 larvae and library preparations were conducted as in McMurchy et al. (2017). Antibodies used were JMJD-2 (antibody Q3951; this study, raised against amino acids 50-149), SCC-1 (Novus antibody 29510002, lot Q0835; raised against amino acids 421-520); ZFP-1 (antibody Q2059; this study, raised against amino acids 1-100), BLMP-1 (antibody Q2919; this study, raised against amino acids 43-142), TOP-2 (antibody Q5515; this study, raised against aa317-437), and SIN-3 (antibody Q6013; Beurton et al. 2019).

Differential expression analysis of met-2 set-25 mutant

Wild-type and met-2 set-25 mutant worms were grown at 20°C on NGM plates, harvested at the L3 stage, and flash frozen in liquid nitrogen. RNA was extracted from frozen worms using TriPure (Roche) and purified with the Zymo Research RNA Clean and Concentrator kit (R1013) after DNase I digestion. Libraries were prepared with the TruSeq RNA Library Prep Kit (Illumina). Reads were aligned to genome assembly ce10 with gene annotation WS235 using STAR (Dobin et al. 2013). Read counts per gene were generated by HTSeq and genes differentially expressed in met-2 set-25 relative to wild type identified using DESeq2 (Love et al. 2014), requiring adjusted P < 0.05 (Supplemental Table S1). Genes were annotated as being marked by H3K9me2 or H3K9me3 if their average signal was 1.5× above the median, using L3 ChIP-seq data sets available under the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) accession numbers GSE49728 for H3K9me2 and GSE49720 for H3K9me3 (Ho et al. 2014).

Differential expression analysis of blmp-1 mutant

Raw RNA-seq data of wild-type and blmp-1(tm548) mutants at L3 stage were obtained from GEO (GSE55225) (Horn et al. 2014). Reads were aligned to genome assembly ce10 with gene annotation WS235 using STAR (Dobin et al. 2013). Read counts per gene were generated by featureCounts, and genes differentially expressed in blmp-1 relative to wild type were identified using DESeq2 (Love et al. 2014), requiring adjusted P < 0.05 and absolute fold change >1.5. Genes significantly up-regulated or down-regulated in blmp-1 mutants compared to wild type and for which a BLMP-1 ChIP-seq peak overlapped one of its assigned regulatory elements were defined as up-regulated “up” or down-regulated “down” BLMP-1 targets (Supplemental Table S6).

ARC-C library preparation

A key element of ARC-C is to digest nuclei with DNase I in situ under conditions in which cutting results in maximal enrichment of informative read pairs at accessible chromatin. This concentration should be empirically determined for the sample type of interest. For C. elegans nuclei, we found that 50–200 units/mL DNase I gives optimal recovery of informative read pairs using the digestion conditions indicated below (Supplemental Fig. S1).

Frozen worm pellets were ground into a fine powder in which worms were broken into approximately 10 fragments. Then 1 mL of worm powder was fixed in 10 mL of 1% formaldehyde in PBS for 10 min at room temperature (RT) with gentle shaking then quenched for 5 min with a final concentration of 125 mM glycine. Fixed worm fragments were then washed with Buffer A (340 mM sucrose, 15 mM Tris-HCl at pH 7.5, 2 mM MgCl₂, 0.5 mM spermidine, 0.15 mM spermine, 1 mM DTT, protease inhibitors), resuspended in 7 mL Buffer A, then the material dounced 20 strokes in a 7 mL stainless steel tissue grinder (VWR 432-5005). The dounced material was spun at 100g for 5 min; the supernatant, which contains nuclei, was transferred to a new tube. The pellet was resuspended in Buffer A and again dounced with 20 strokes. After spinning, the two supernatants containing nuclei were pooled.

Aliquots of 10 million fixed nuclei were spun down at 1000g and resuspended in 200 µL of 1× DNase buffer (Roche), and chromatin was digested with 50 and/or 100 units/mL DNase I for 10 min at 25°C. The reactions were then quenched with a final concentration of 25 mM EDTA and 5 mM Tris at pH 7.5. Nuclei were washed twice with 1 mL of ice-cold Nuclear Washing Buffer (340 mM sucrose, 15 mM Tris-HCl at pH 7.5, 25 mM EDTA, 0.5 mM spermidine, 0.15 mM spermine, 1 mM DTT, protease inhibitors). Nuclei were then resuspended in 100 µL of end repair master-mix (10 µL 10× NEB End Repair Buffer, 5 µL NEB End Repair Enzyme Mix, 85 µL H₂O) and incubated for 30 min at 20°C with rotation. Thereafter, 400 µL of ligation master-mix was added (40 µL 10× ligation buffer, 5 µL T4 DNA ligase [400,000 units/mL], 355 µL H₂O) and the mixture was incubated overnight at 4°C with rotation. Nuclei were pelleted, resuspended in 50 µL of tagmentation master-mix (22.5 µL H₂O, 25 µL 2× Nextera TD buffer, 2.5 µL Nextera Tn5 transposase), and incubated for 30 min at 37°C. Next, 5 µL of 1% SDS and 2 µL of NEB Proteinase K [800 units/mL] was added and left for 15 min at 65°C before DNA from the mixture was purified using Qiagen MinElute columns. Large fragments (>500 bp) were removed from the purified DNA using two rounds of 0.6vol AMPure XP beads. Before PCR amplification, a test qPCR amplification using 1/20th of the input DNA was performed, and the cycle number for PCR was determined by the midpoint of the exponential phase of the amplification curve. The resultant DNA was amplified with NEBNext Ultra II Q5 Master-mix under the following PCR conditions: for 5 min at 72°C, for 30 sec at 98°C, and cycling for 10 sec at 98°C, for 30 sec at 63°C, and for 1 min at 72°C using the determined cycle number. The library was size-selected with AMPure XP beads to a final range of 200–700 bp (insert size: ∼70–570 bp). Before sequencing, libraries were quality tested by qPCR at the gap-3 promoter. DNase I overdigestion or underdigestion led to poor enrichment of regulatory elements. Good libraries generally had gap-3 promoter enrichment values of greater than fourfold. Libraries were sequenced on the Illumina platform paired-end (100 bp). Four libraries were sequenced from wild-type L3 larvae (three biological replicates and two technical replicates; separate ARC-C libraries made from the same biological material: N2-rep1, N2-rep2a, N2-rep2b, and N2-rep3 in Supplemental Fig. S2A), and two biological replicate libraries were sequenced for both blmp-1 and met-2 set-25 L3 larvae.

Processing ARC-C data

ARC-C libraries were sequenced 62–150 bp (100 bp for most samples) from both ends. Supplemental Figure S2 lists libraries analyzed in the paper and sequencing statistics. Adapter sequence was trimmed by cutadapt (Martin 2011), and sequences with >20 bp remaining were removed. Each sequenced end was aligned independently to the ce10 reference genome using BWA-MEM (Li 2013), which allows split-read alignment using the default parameters. The two aligned ends (or the 5′ segment for split alignments) were then paired. We required both ends of a pair to align uniquely and with high confidence (mapping quality ≥ 30 and number of mismatches ≤ 2) to the nuclear genome and outside modENCODE backlisted-regions. PCR duplicates were next removed by sambamba markdup. The remaining read pairs were regarded as valid read pairs. Valid read pairs mapping to different chromosomes or >600 bp apart on the same chromosome were regarded as trans- or cis-informative read pairs, respectively. The 600-bp threshold was established by comparing the proportions of the four possible end alignment orientation configurations (forward–forward, forward–reverse, reverse–forward, and reverse–reverse) as a function of mapping distance. The vast majority of pairs mapping <500 bp apart were in the forward–reverse configuration (nonligated fragments), whereas >600 bp the proportions were stably at ∼25% each (Supplemental Fig. S12).

Contact maps were made from informative read pairs by binning the genome into fixed-width (1, 5, 10, or 50 kb) nonoverlapping bins and counting the number of read pairs between each pair of bins. The maps were then normalized by matrix balancing using the Knight–Ruiz algorithm (Knight and Ruiz 2013). Concordance between replicate maps was assessed by GenomeDISCO (Ursu et al. 2018). The lower resolution 50-kb maps were used for whole-chromosome visualization. The 1-, 5-, and 10-kb maps were further corrected for distance-dependent background contact frequency by dividing the spline-smoothed average contact frequency given the distance to the diagonal. These higher resolution maps were used for aggregated contact analysis (see below) of nuclear factor binding sites and chromatin domain/compartment, respectively.

Processing Hi-C data

Raw FASTQ files of wild-type mixed embryo Hi-C data (Crane et al. 2015) were downloaded from the NCBI Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra) accession number SRX77040 and processed using HiCUP v0.5.9 (Wingett et al. 2015), which filters for same fragment (circularized, dangling ends, internal), religation, wrong size, contiguous sequence, and removes duplicate read pairs. We additionally required mapQ ≥ 30 and number of mismatches ≤2 from both reads, that none of the reads overlapped modENCODE blacklisted regions, and a minimum distance of 600 bp between the two read pairs to be consistent with the processing of ARC-C data. In the end, 25,460,294 read pairs passed all filters, of which 17,100,808 have both reads mapping to the same chromosome.

Calling significant interactions

We first segmented the genome into bins of ∼500 bp using the following procedure. We first took annotated regulatory elements (Jänes et al. 2018) (n = 42,245) and expanded them to 500 bp or until neighboring intervals began to touch; a small number of elements that were within 100 bp were merged first. The rest of the genome was covered with evenly placed 500-bp nonoverlapping fixed-width intervals; hence, the entire genome was covered by a combined set of 192,257 intervals of average size (494 bp). We used this procedure instead of generating fixed nonoverlapping bins to avoid individual Res being split into two bins.

Before assessing the significance of chromatin interactions in Hi-C or other chromatin interaction data sets, inherent bias resulting from uneven coverage and physical proximity need to be accounted for. A difference in ARC-C data compared to that of Hi-C is the specific overrepresentation of open chromatin. The observed enrichments at these regions are a compound effect of enriching the regions themselves and of their being enriched by virtue of linkage to another region of open chromatin (true signals). To address this and avoid normalizing out interaction signals, we developed an approach similar to that used for Capture Hi-C, which has similar coverage biases and used “off-peak” interaction frequency to account for coverage (Cairns et al. 2016). We modeled the number of read pairs linking two intervals as a random variable following a binomial distribution parameterized by an expected contact frequency determined by unevenness of coverage and distance between the interval (Ay et al. 2014), as described next.

For every interval i∈I, the number of cis-informative read pairs c_cis,i were counted. Intervals in the top 10% of the coverage distribution were regarded as peaks, and intervals in the bottom 10% were removed. An off-peak cis-informative coverage c_offpeak,i was calculated for every kept interval, counting the number of contacts not involving peak intervals. We calculated a scaling factor for the interval's representation/visibility as v_i = (c_offpeak,i/median(c_offpeak,.))^0.87 (see below for derivation of the exponent). Chromosome-wide average distance-dependent contact frequency F(d) in the distance range of 1 kb to 1 Mb was modeled by fitting a spline function in a two-pass process (following Fit-Hi-C) (Ay et al. 2014). For every pair of intervals with a distance between 1 kb and 1 Mb, an expected contact frequency was calculated given the distance and the visibility of each interval as f_i,j = v_i v_j F(d_i,j). Given the total number of cis-informative contacts (N) of the chromosome, we considered a null distribution in the form of a binomial, where the observed number of contacts, n_i,j ∼ binomial(N, f_i,j). Significant interactions were called at an FDR level of 0.05 and were post-filtered requiring support by more than five read pairs.

The appropriate correction factors for adjusting the representation bias, X, should be able to transform the unnormalized contact matrix A into a normalized contact matrix B by B = diag(X)*A*diag(X), such that each row or column of B sums up to the same value, thus eliminating the unevenness of representation across different bins. It has been well-established that a correct set of factors can be found using the method “matrix balancing” (MB) (Imakaev et al. 2012; Rao et al. 2014), and an efficient algorithm has been developed (Knight and Ruiz 2013; Rao et al. 2014). Our goal is to find correction factors that enable transformation into a normalized contact matrix in which each bin has the same off-peak coverage. Because an efficient algorithm for solving this problem at high resolution is not yet available, we aimed to find an approximate solution. Two measures, namely, the reciprocal of coverage (c) and square root coverage (c^0.5), have been proposed for use in the place of MB-derived correction factors (Rao et al. 2014). It was reported that the former can overcorrect, whereas the latter gives good approximation to MB-derived correction factors (Rao et al. 2014). We examined correlation between MB-derived correction factors and the reciprocal of different exponents of coverage (Supplemental Fig. S13) and found that the reciprocal of c indeed overcorrects, but the reciprocal of square root coverage undercorrects. The most accurate approximation is achieved at around an exponent of 0.87. Therefore, we used the reciprocal of (c_offpeak)^0.87 as the correction factor.

Aggregated contact analyses

A contact is defined as the region in the contact map that connects a pair of genomic locations. Aggregated contact analysis is a method of visualizing the average contact frequency of a group of many contacts together with local contact frequency (Rao et al. 2014). We applied this method to both small genomic regions, such as nuclear factor ChIP-seq binding sites (NFBS), and to larger intervals, such as chromatin domains. We normalized contact maps using matrix balancing (Imakaev et al. 2012) to account for coverage bias and removed distance-dependent background. In NFBS analysis, we used normalized and background-frequency-corrected contact maps of 1-kb resolution. For each NF, up to 50,000 contacts were randomly sampled from all possible cis-contacts among its binding sites within a distance range from 20 kb to 1 Mb, and local maps of 21 × 21 bins centered at the contacts were extracted and aggregated. For the case of BLMP-1 regulated targets, all possible cis-contacts between BLMP-1 binding sites that involves a BLMP-1 regulated target were aggregated (i.e., at least one end of the contact is at a BLMP-1 regulated target). The log₂ fold change of the central point over the mean of the rest of the points in the aggregated map was calculated to measure the relative increase in contact frequency over local background. To assess statistical significance while controlling for accessibility and the distance between the pair of NFBSs, 1000 sets of random contacts were generated, each containing the same number of contacts with matching accessibility and distance as the NFBS contacts. Each of the 1000 random sets was aggregated and a relative increase in contact frequency was calculated in the same way as the NFBS set, forming a distribution of values against which the NFBS value was compared and a P-value generated, which was corrected for multiple testing using the FDR method. The rex–rex APA analyses used normalized and background-frequency-corrected contact maps of 10-kb resolution and a distance range of 100 kb to 4 Mb. Data sets (Kranz et al. 2013; Araya et al. 2014; Ho et al. 2014; Latorre et al. 2015; Wiesenfahrt et al. 2016; McMurchy et al. 2017; Kudron et al. 2018) are listed in Supplemental Table S7. Processed and curated modENCODE ChIP-seq peaks (Araya et al. 2014; Kudron et al. 2018) were obtained from Jänes et al. (2018). For ChIP-seq performed in this paper (JMJD-2, SCC-1, ZFP-1, BLMP-1, TOP-2, and SIN-3), peaks were called for each replicate separately as for accessible sites in Jänes et al. (2018) using YAPC (https://github.com/jurgjn/yapc), except for BLMP-1, which used MACS2 (Feng et al. 2012). ChIP-seq peaks were combined by IDR (Li et al. 2011) with a P-value threshold of 0.01. We removed intervals that are highly occupied target (“HOT”) because such binding events are thought to represent non-sequence-specific TF binding or ChIP artifacts (Gerstein et al. 2010; Kudron et al. 2018). This was defined as the top 20% of peak intervals ranked by the number factors in which the interval is called (effectively removing intervals called in 12 or more factors). We only considered data sets having at least 300 peaks following filtering.

For domain analyses, we used normalized and background-frequency-corrected contact maps of 5-kb resolution. L3 stage chromatin domains are from Evans et al. (2016). Regions annotated as active and border were merged to generate the active domains used here; the H3K27me3 domains are those termed “regulated.” Informative reads mapped at similar levels to active and H3K27me3 domains: 10.6 million (44.4%) of all informative reads mapped to active domains (39.0 Mb, 38.9% of the genome), whereas 8.1 million (33.8%) mapped to H3K27me3 domains (42.3 Mb, 42.1% of the genome). The rest of the reads mapped to Chromosome X, for which active and H3K27me3 domains were not mapped. The small (1.42-fold) bias in coverage was normalized by the matrix balancing that was performed on the contact frequency matrix to normalize the differences in sequencing coverage across all genomic bins.

We tested for TADs by assessing intra-domain contacts. For each type of domain (active or H3K27me3), maps of each contact region containing a domain of at least 5 kb together with up to 25 kb of the flanking domains of opposite type were extracted and aggregated. Where neighboring domains were <25 kb, only the domain was extracted. The aggregated contact was scaled to a square of 9 × 9 bins, and the flanking intervals were scaled to five bins wide. The log₂ fold change of the mean of the central 9 × 9 square over the mean of the four neighboring 5 × 10 rectangles on top, bottom, left, and right was calculated to measure the relative strength in contact frequency with P-values generated by t-test. We tested for compartments using the same approach, by assessing all possible pairs of inter-domain contacts in the range from 50 kb to 2 Mb. To compare domain and compartment strength in met-2 set-25 mutants relative to wild type, the fold change in percentage difference between central block and side blocks across strains were calculated.

To test the effect of local H3K9 methylation on TAD and compartment strength, we separated active and H3K27me3 domains into highly and lowly marked sets. For each active and H3K27me3 domain, the average coverage of H3K9me1, H3K9me2, and H3K9me3 was calculated using wild-type L3 ChIP-seq data (Ho et al. 2014). For each modification, active and H3K27me3 domains in the top 25% of signal were defined as H3K9me1, H3K9me2, or H3K9me3 high domains, and the remainder as low domains. Aggregated contact analysis of TADs and compartments was performed using wild-type and met-2 set-25 ARC-C data for each of the domain sets as described above.

Aggregated contact analysis of blmp-1 mutant ARC-C data

We used the procedure described above for wild type to measure the relative contact frequency between NFBSs in blmp-1 ARC-C data. To normalize overall open chromatin measurements between blmp-1 and wild-type maps, we first fitted blmp-1 ACA measures as a linear function of the respective measures in wild type, that is, log₂(blmp1_FC) = a + b * log₂(N2_FC). The residuals were then used to measure the difference in contact frequency in blmp-1 relative to wild type. Statistical significance was assessed by 10,000-time bootstrapping the distribution of residuals, and P-values were adjusted by FDR (Supplemental Table S5).

Expression correlation between interacting promoters

For genes linked by promoter–promoter interactions, we calculated Pearson correlation coefficients of gene expression across cell types using data from Cao et al. (2017). For pairs involving bidirectional promoters, the gene pair with highest correlation was chosen. For each gene, we also calculated a coefficient of variation of gene expression (CV) across the cell types (Cao et al. 2017) as a measure of tissue-biased expression. Genes with similar expression across cell types have low gene expression CV values and those with tissue-biased expression have high CV values. In Figure 3C, we assessed expression correlation between all linked pairs of genes (n = 5124), linked genes for which both were in the bottom 30% of all CV values (n = 2224), and linked genes for which both were in the top 70% of CV values (n = 879). To assess the statistical significance of the correlations, we generated control sets by sampling the same number of random pairs using genes from the respective linked gene sets while requiring the distribution of distance between the random pairs to match that of the observed set. Difference in correlation between the observed and the random control set was tested using a t-test: all (P = 7.77 × 10⁻²⁷), bottom 30% (P = 4.38 × 10⁻⁰⁶), and top 70% (P = 5.69 × 10⁻²⁰).