Impact of 3D genome organization, guided by cohesin and CTCF looping, on sex-biased chromatin interactions and gene expression in mouse liver

Matthews, Bryan J.; Waxman, David J.

doi:10.1186/s13072-020-00350-y

Research
Open access
Published: 17 July 2020

Impact of 3D genome organization, guided by cohesin and CTCF looping, on sex-biased chromatin interactions and gene expression in mouse liver

Epigenetics & Chromatin volume 13, Article number: 30 (2020) Cite this article

5088 Accesses
13 Citations
5 Altmetric
Metrics details

Abstract

Several thousand sex-differential distal enhancers have been identified in mouse liver; however, their links to sex-biased genes and the impact of any sex-differences in nuclear organization and chromatin interactions are unknown. To address these issues, we first characterized 1847 mouse liver genomic regions showing significant sex differential occupancy by cohesin and CTCF, two key 3D nuclear organizing factors. These sex-differential binding sites were primarily distal to sex-biased genes but rarely generated sex-differential TAD (topologically associating domain) or intra-TAD loop anchors, and were sometimes found in TADs without sex-biased genes. A substantial subset of sex-biased cohesin-non-CTCF binding sites, but not sex-biased cohesin-and-CTCF binding sites, overlapped sex-biased enhancers. Cohesin depletion reduced the expression of male-biased genes with distal, but not proximal, sex-biased enhancers by >10-fold, implicating cohesin in long-range enhancer interactions regulating sex-biased genes. Using circularized chromosome conformation capture-based sequencing (4C-seq), we showed that sex differences in distal sex-biased enhancer–promoter interactions are common. Intra-TAD loops with sex-independent cohesin-and-CTCF anchors conferred sex specificity to chromatin interactions indirectly, by insulating sex-biased enhancer–promoter contacts and by bringing sex-biased genes into closer proximity to sex-biased enhancers. Furthermore, sex-differential chromatin interactions involving sex-biased gene promoters, enhancers, and lncRNAs were associated with sex-biased binding of cohesin and/or CTCF. These studies elucidate how 3D genome organization impacts sex-biased gene expression in a non-reproductive tissue through both direct and indirect effects of cohesin and CTCF looping on distal enhancer interactions with sex-differentially expressed genes.

Introduction

Sex differences in gene expression are found in several non-reproductive tissues, including the brain [1], immune system [2], kidney [3] and liver [4]. Sex differences in the liver are associated with a higher incidence of aggressive liver cancer in males [5], increased susceptibility to autoimmune hepatitis in females [6], and sex differences in metabolism of diverse pharmaceuticals and environmental chemicals [7]. Transcriptomic and epigenetic sex differences in the transcriptome are best characterized in mouse liver, where more than 1000 genes [8], including many lncRNA genes [9, 10] and miRNAs [11], exhibit sex-biased expression regulated by the sex-differential temporal patterns of pituitary growth hormone secretion [12]. Sex differences in the epigenome are widespread, and frequently are associated with sex differences in gene distal, but not gene proximal, regulatory elements, which show characteristic sex-differential patterns of histone marks and chromatin accessibility (DNase hypersensitive sites, DHS) [13, 14]. Three-dimensional looping is one mechanism that could potentially link the few thousand mostly distal sex-biased enhancers identified to individual sex-biased genes.

CCCTC-binding factor (CTCF) and the multi-protein complex cohesin are two major transcription factors regulating 3D genomic architecture. CTCF has primarily been studied for its role in DNA looping and insulation [15, 16], while cohesin is a molecular motor powering DNA looping via a loop extrusion mechanism [17, 18]. Loss of CTCF or cohesin is lethal in developing mouse embryos [19, 20]. However, when degradation of the cohesin loading factor Nipbl is induced in adult mouse liver, a dose-dependent loss of both cohesin binding and virtually all focal DNA looping is seen without major hepatocyte toxicity [21]. Loss of DNA looping also occurs in other systems following depletion of either cohesin [22] or CTCF [23]. Thus, CTCF and cohesin are both required for DNA looping.

The functional role of cohesin at a given genomic site is largely dependent on its binding partners. Cohesin lacks sequence-specific DNA binding activity, but is loaded and unloaded from chromatin by specific protein complexes [24]. Cohesin can participate in shorter-range looping between enhancers and promoters (‘enhancer–promoter loops’) in association with Mediator and tissue-specific transcriptional regulators [25,26,27]. Genomic sites bound by cohesin but not CTCF, i.e., cohesin-non-CTCF sites (CNC), tend to be highly tissue specific, but frequently show weaker binding than sites where cohesin and CTCF are both bound, i.e., cohesin-and-CTCF (CAC) sites [25, 28, 29]. Cohesin forms insulating loops at CAC sites, which have been variously characterized as Topologically Associating Domains (TADs) [30] and intra-TAD loops [31], loop domains [32] and insulated neighborhoods [26]. A majority of CAC-mediated insulated loops are conserved across tissues and cell types [30,31,32,33] and many are even conserved across mammalian species at syntenic regions [34]. A distinct subset of CAC sites is important for enhancer–promoter interactions [35, 36], and may specifically help structure super-enhancers around key constituents enhancers, known as hubs [37]. The role of CTCF binding in the absence of cohesin, i.e., at Lone CTCF sites, is less clear. CTCF binds specific DNA motifs via its 11 zinc fingers [15], yet the function of Lone CTCF sites is likely also dependent on additional interacting proteins, such as the transcription factors YY1 [38] and STAT5 [39].

CAC-mediated insulating loops are typically anchored by a pair of tightly bound CTCF sites oriented toward each other by the non-palindromic CTCF motif [32]. Supporting this finding, inversion or deletion of CAC-bound anchors results in a loss of specific loops and subsequent rewiring to form alternative loops [17, 40]. This stereotypical pattern of convergently oriented CTCF binding enables the accurate prediction of CAC-mediated loops based on CTCF and cohesin binding activity and motif orientation alone [17, 31, 41]. Computational prediction of loops linking enhancers and promoters, i.e., gene loops, is a more elusive goal, although some recent progress has been made [42, 43]. Chromatin conformation capture technology using formaldehyde crosslinking and restriction enzyme digestion followed by proximity ligation can be employed to identify such loops experimentally, and thereby determine interaction frequencies between different genomic regions [44]. Circularized chromatin conformation capture with sequencing (4C-seq) is one such method that interrogates all potential interactions between a single site of interest and the rest of the genome [45].

Here we take a multi-pronged approach to elucidate the role of architectural proteins and 3D genome organization in regulating the widespread sex differences in gene expression seen in mouse liver. First, we identify sex-biased binding sites for cohesin and CTCF in mouse liver chromatin, a majority of which were found at intergenic sites distal to sex-biased genes. Further, we investigate the effects of a deficiency in cohesin loading in male mouse liver [21], and find that cohesin is specifically required for expression of male-biased genes with distal sex-biased regulatory elements. Finally, we use 4C-seq to directly evaluate sex differences in chromatin interactions involving sex-biased enhancers from five different genomic regions, and demonstrate the importance of loop domains for insulation of enhancer–promoter contacts at sex-biased genes. Overall, our findings highlight how 3-dimensional genome organization contributes to sex differences in liver gene expression in both direct and indirect ways.

Materials and methods

Mouse protocols, extraction of liver nuclei, and chromatin preparation

Male and female CD-1 mice were purchased from Charles River Laboratories [strain # Crl:CD1(ICR)] and housed separately in the Boston University Laboratory Animal Care Facility. All animal protocols were specifically reviewed for ethics and approved by the Boston University Institutional Animal Care and Use Committee. Mice were euthanized by cervical dislocation at 8 weeks of age, livers were excised, and nuclei were purified then crosslinked with 0.8% formaldehyde for 9 min at 30 °C [31]. The crosslinked chromatin was sonicated and stored at − 80 °C after a small aliquot was removed for crosslink reversal (6 h at 65 °C in a thermocycler with a heated lid) to quantify the chromatin by Quanti-iT PicoGreen assay (Invitrogen, cat. # Q33130) and determine the fragment size distribution, as described [31].

Chromatin immunoprecipitation and sequencing (ChIP-seq)

Immunoprecipitation of sonicated mouse liver chromatin was performed as described [14]. Specifically, 5 µl of rabbit polyclonal antibody to either CTCF (Millipore, cat. # 07-729) or to the cohesin subunit Rad21 (Abcam, cat. # ab992) was mixed with 30 µl of Protein A Dynabeads (Invitrogen, cat. # 1002D) and incubated in blocking solution (0.5% bovine serum albumin in PBS) for 3 h at 4 °C. Beads were then washed with blocking solution, followed by incubation overnight with 70 µg of sonicated liver chromatin in 1X RIPA buffer (50 mM Tris–HCl, pH 8.0, 150 mM NaCl, 1% IPEGAL, 0.5% deoxycholic acid) containing 0.1% SDS. After washing with 1X RIPA (containing 0.1% SDS), formaldehyde crosslinks were reversed for 6 h at 65 °C, followed by RNase A digestion (Novagen, cat. # 70856) at 37 °C for 30 min and then proteinase K digestion (Bioline, cat. #37084) for 2 h at 56 °C. The resulting crude DNA extract was purified using a QIAquick Gel Extraction Kit (Qiagen, cat. # 28706) and quantified with a Qubit HS DNA kit (Invitrogen, cat. # Q32854). All samples were processed using the same protocol and conditions. Sequencing was performed for a total of eight CTCF ChIP-seq samples (n = 4 individual male and n = 4 individual female livers) and a total of six Rad21 ChIP-seq samples (n = 3 male, n = 3 female livers). The male liver ChIP-seq samples were those reported previously [31] and are available at GSE102997. Female liver ChIP samples are available at GSE130908.

4C-seq methods

Isolation of liver nuclei and crosslinking were performed as described [31] through the step where crosslinked nuclei were pelleted by centrifugation. Digestion, proximity ligation, and inverse PCR were then carried out as described [31]. Briefly, frozen nuclei were resuspended in Buffer A (15 mM Tris–HCl pH 8.0, 15 mM NaCl, 60 mM KCl, 1 mM EDTA pH 8.0, 0.5 mM EGTA pH 8.0, 0.5 mM spermidine, 0.3 mM spermine) and quantified using a Countess Automated Cell Counter (Invitrogen, cat. # C10227). An aliquot of 10 million nuclei was used for each individual 4C experiment. Nuclei were pelleted at 1000×g for 5 min at 4 °C then resuspended in 450 μl of 1X NEBuffer 3 (NEB, cat. # B7003S). Primary restriction enzyme digestion of intact nuclei was carried out overnight at 37 °C using 50,000 U of DpnII (NEB: #R0543) with agitation at 900 RPM. DpnII was inactivated by adding SDS to a final concentration of 2%. Samples were then diluted fivefold in 1X ligation buffer (Enzymatics, cat. # B6030). Proximity ligation was performed overnight at 16 °C with 200 U of T4 DNA ligase (Enzymatics, cat. # L6030). DNA was then reverse-crosslinked and purified using a standard phenol/chloroform cleanup method following the manufacturer’s protocol (VWR, cat. # VWRV0883). Secondary digestion of purified DNA was performed overnight at 37 °C with 50 U of Csp6I (Thermo Scientific, cat. # ER0211) in 500 µl of 1X Buffer B (Fermentas, cat. # BB5; 10 mM Tris–HCl (pH 7.5), 10 mM MgCl₂, 0.1 mg/ml BSA), followed by heat inactivation at 65 °C for 30 min. Samples were diluted tenfold and secondary ligation was carried out overnight at 16 °C, as described above. The effectiveness of primary digestion, proximity ligation, secondary digestion, and secondary ligation was verified at each step by reverse crosslinking and gel electrophoresis analysis (see gel image in Additional file 1: Figure S5A). The final PCR template was purified by phenol/chloroform clean up, followed by QiaPrep 2.0 column cleanup (Qiagen, cat. # 27115) to yield a standard, circularized 4C inverse PCR template, which was amplified using specific viewpoint primers, as described below.

PCR reactions were performed using inversely oriented primer pair sequences for valid 4C-seq viewpoints, obtained from the 4CSeqpipe primer database [http://www.wisdom.weizmann.ac.il/~atanay/4cseq_pipe/], with the addition of 5′ dangling truncated Illumina adapters (Additional file 2: Table S3A). Candidate viewpoints were selected based on the following criteria. First, we only considered viewpoints that are in the same TAD as at least one protein-coding or lncRNA gene showing > 3-fold sex bias in its expression. Second, the viewpoint must be within 1 kb of the transcription start site (TSS) of a sex-biased gene, or it must overlap a sex-biased enhancer (minimum twofold sex-bias in normalized DHS opening or H3K27ac mark intensity). Third, the non-reading primer (Additional file 2: Table S3A) was required to map to the genome uniquely, while the reading primer was more stringently required to have > 89% unique sequence identity (i.e., no 18-mer within a 20 nt primer sequence that maps elsewhere in the genome). Inverse PCR amplification of 1 µg of each 4C template was performed using Platinum Taq DNA polymerase (Invitrogen, cat. #10966026), as follows: 94 °C for 2 min, 25 amplification cycles (94 °C for 30 s, 55 °C for 30 s, 72 °C for 3 min), then 4 °C hold. To minimize the impact of PCR artifacts, a total of 6 identical inverse PCR reactions were performed for each sample, except as noted. These replicate reactions were processed in parallel and pooled prior to library preparation. For 4C-seq analysis of the Nox4 genomic region viewpoints, single independent inverse PCR reactions were carried out for each of four liver samples (males M1 and M2; females F1, and F2) and sequence libraries were then prepared independently. Sequencing data was pooled at the raw read level. Examples of pooled 4C libraries are shown in Additional file 1: Figure S5B.

4C-seq samples were amplified using barcoded primers, such that each biological replicate had a unique barcode (NEB, cat. # E7335). To minimize over-amplification, each sample was amplified with 5 additional cycles of PCR to fill in the full Illumina adapter sequence needed for sequencing. Due to the low sequence complexity in the first 20 bases of each sequence read, 4C-seq libraries were multiplexed by combining with high complexity sequencing libraries for unrelated samples (e.g., RNA-seq or ChIP-seq libraries), which were sequenced in the same Illumina sequencing lane. In practice, 4C-seq libraries constituted no more than 15% of the total library pool, by molarity. Samples were sequenced on an Illumina HiSeq2500 instrument for 50, 125, or 150 bp paired end reads.

Computational analysis of ChIP-seq datasets

Sequence reads were split by barcode and mapped to mouse genome assembly mm9 using Bowtie2 (v2.2.9). All reads not uniquely mapped to the genome were excluded from downstream analyses. Peak calling was performed using MACS2 (v2.1.1) with default parameters, and peaks that overlapped blacklisted genomic regions (http://www.sites.google.com/site/anshulkundaje/projects/blacklists) were filtered out. Additionally, we removed spurious peaks that exclusively contained PCR duplicated reads, defined as 5 or more identical sequence reads that do not overlap any other reads. All BigWig tracks used to visualize sequencing data in a genome browser were normalized for both sequencing depth and sample quality, expressed as reads in peaks per million mapped reads (RIPM). In practice, the browser y-axis displays the read count from a given sample divided by the total number of reads in peaks (reads that overlap a peak identified in any sample, or the union peak list), per million. Normalization was performed separately for the CTCF and cohesin datasets. This approach is functionally similar to the quality control metric known as Fraction of Reads in Peaks (FRiP) used by the ENCODE consortium [46]. All samples used in this study were judged to be of good quality, with a mean FRiP value of 0.217 and ranging from 0.103 to 0.344. A full listing of samples sequenced and sequencing statistics is provided in Additional file 2: Table S3B.

To identify sex-differential ChIP-seq peaks, diffReps (v1.55.4) [47] was used with default parameters and a window size of 200 bp to identify in-peak differential sites, i.e., diffReps sites that overlap a MACS2 peak, defined below. The diffReps output list of sites was filtered to remove diffReps-identified sites that did not meet the following conditions: overlap with at least one of the peaks from the union peak list for the relevant factor (Additional file 3: Table S1F, G for CTCF; Additional file 3: Table S1H, I for cohesin), contains at least 10 sequence reads, shows > 2-fold sex difference, and has an FDR < 0.05. The resultant sets of ChIP-seq peaks were defined as standard stringency sex-biased peaks, and were used for all analysis, except as noted. A set of lenient stringency sex-biased CTCF and cohesin peaks was defined, as follows. Sequence reads for biological ChIP-seq replicates were combined (merged) to give a single merged sample for each sex and each factor (Male CTCF, Male cohesin, Female CTCF, and Female cohesin). For each transcription factor, the male and female merged samples were then compared, and the merged sample with the higher FRiP was down-sampled, so that the mapped read files for each sex contained the same, normalized number of reads in peaks. In practice, the combined mapped reads for the merged female cohesin samples, and for the merged male CTCF samples, were down-sampled by a factor of 0.979869 and 0.745955, respectively. MAnorm [48] was then used to compare the FRiP-normalized male and female samples to identify a set of sex-differential binding sites for each factor (minimum twofold sex bias; p-adj < 0.01; read count > 15 for the upregulated peak). Binding sites identified by MAnorm that were not on the standard sex-biased peak list were designated lenient stringency sex-biased peak lists. Bedtools (v2.26.0) was used for overlap analysis to determine the distance of ChIP-seq peaks to other genomic features. Genomic coordinates (mm9) for TAD and intra-TAD boundaries were downloaded from [31], where TAD definitions are based on experimental Hi-C analysis for male mouse liver [34]. Given the low resolution of this dataset (20–40 kb bin size), the TAD definitions used in this study are most similar to early Hi-C studies, and may underestimate the total number of domains in liver.

Sex-differential ChIP-seq peaks were analyzed separately for the sets of male and female biological replicates and then divided into four groups, CAC(∆Coh), CAC(∆CTCF), CNC(∆Coh), and Lone CTCF(∆CTCF), based on the following criteria. Sex-differential cohesin peaks (∆Coh) were designated CAC(∆Coh) peaks if they overlapped a CTCF peak identified in at least 3 of the 4 CTCF ChIP-seq biological replicates in the sex showing higher cohesin binding (e.g., male-biased cohesin peaks were compared to male liver CTCF ChIP-seq replicates). Alternatively, they were designated CNC(∆Coh) peaks if they overlapped a CTCF peak found in either 0 or 1 of the four CTCF-seq biological replicates. Those ∆Coh peaks that overlapped 2 of the 4 CTCF replicates were excluded from downstream analyses. Similarly, sex-differential CTCF peaks (∆CTCF) were designated CAC(∆CTCF) peaks if they overlapped a cohesin peak found in 2 or 3 of the three available cohesin ChIP-seq biological replicates. Sex-differential CTCF peaks were designated Lone CTCF(∆CTCF) peaks if they overlapped peaks in either 0 or 1 of the three cohesin ChIP-seq biological replicates. 137 sex-differential CTCF peaks overlapped sex-differential cohesin peaks, and were thus CAC(∆Coh/∆CTCF); 50 of these 137 CAC peaks were autosomal (Additional file 3: Table S1C). Sex-independent cohesin peaks, and sex-independent CTCF peaks, were, respectively, defined by ranking each peak based on the following ratio: (RIPM-normalized ChIP signal for the merged male sample)/(RIPM-normalized ChIP signal for the merged female sample), performed separately for CTCF and for cohesin. The 1000 peaks whose ratios were closest to 1 were defined as the set of sex-independent CTCF, and cohesin, peaks.

Discovery of intra-TAD loops

CTCF motif discovery was performed using the FIMO option from MEME Suite (v4.10.0), and presence of a motif was defined as a motif score > 10. Intra-TAD loops for female mouse liver were identified using the computational method described previously for male liver [31]. The analysis pipeline was run for all CAC sites and with an initial loop count of 20,000, using the set of default parameters reported for male liver [31]. This analysis yielded 9724 intra-TAD loops with 10,273 loop anchors in female liver; this compares to 9543 intra-TAD loops and 9052 loop anchors identified in male liver [31]. The redundancy in loop anchors is a reflection of nested CAC-mediated loop structures, as was described in other studies using experimentally measured loop identification compared to the computational approach used here; these studies include ChIA-PET analysis of the cohesin subunit SMC1A in mouse embryonic stem cells [26] and Hi-C analysis in human GM12878 cells, where 9448 loops were associated with 12,903 loop anchors [17, 32]. Reciprocal overlap between loops was analyzed using bedtools (bedtools intersect −wa −u −r −f 0.8), as described [31]. A total of 2527 intra-TAD loops were unique to either male or female liver; however, very few had anchors that overlapped a sex-differential CAC site, suggesting that most are not biologically relevant. Supporting this, the loops that were unique to either male or female liver were weaker than the loops shared between male and female livers, and in many cases the loops narrowly met the significance cutoff in one sex but not the other. This finding is similar to our earlier finding that tissue-specific loops are often weaker than those predicted in multiple tissue types [31]. Intra-TAD loops for male and female mouse liver are listed in Additional file 3: Table S1J, and female intra-TAD loop anchors are listed in Additional file 3: Table S1K; a comparable listing for male liver is available in [31].

Computational analysis of 4C-seq datasets

Biological replicates were demultiplexed by index read barcode. As the fastq files for each biological replicate contained sequence reads from multiple viewpoints, the reads in each file were further split based on matches to the reading primer for each viewpoint (Additional file 2: Table S3A). Then, prior to mapping, we used FASTX-Toolkit (v0.0.14) to remove the first 20 nt of sequence from the 5′ end of each read, as this represents the reading primer. For read length consistency, we trimmed 25 nt from the 3′ end of 150 nt read libraries, making them identical in length to the 125 nt libraries (i.e., both were 105 bases long after 5′ and 3′ trimming). Bowtie2 (v2.2.2) was then used to iteratively map the reads to a reduced mouse genome, which comprised all genomic sequences 105 bp upstream and 105 bp downstream from each DpnII cut site in the genome (recognition sequence: GATC), as in [45, 49]). To implement this step, we scanned the genome for all occurrences of the DpnII cut site GATC using the UCSCutils tool oligoMatch (–exactmatch; default parameters). The resulting set of coordinates was then expanded to include sequence 105 bp upstream and 105 bp downstream of each DpnII site (bedtools flank − l 105 − r 105), which was extracted from the mm9 genome using bedtools getfasta. This reduced mouse genome sequence was indexed using bowtie2 prior to mapping. An iterative mapping strategy was implemented because some reads contained multiple ligation junctions, rendering the full-length read unmappable, as described in some Hi-C pipelines [50]. First, non-uniquely mapped reads were trimmed by ~ 10% of their total length (starting from the 3′ end of the read), and mapping was reattempted as above. This process was repeated until a minimum read length of 20 nt was reached, with a step size of 2 nt (for 50-nt-long reads) or a step size of 10 bp (for 125 or 150-nt-long reads) per iteration. These iterative trimming and remapping steps increased the overall percentage of uniquely mapping reads by 2–6% for shorter read lengths libraries, and by 8–22% for longer read length libraries.

After mapping, bedGraph files were smoothed using ucscutils (v. 20130327) with a sliding window of 11 restriction fragments, taking the median value in the window. Smoothed BigWig tracks were normalized by total reads per million to account for differences in sequencing depth. Merged replicates shown in the main figure panels were generated by taking the median signal at a given restriction fragment per viewpoint and sex. The genomic region overlapping the viewpoint fragment was removed prior to BigWig generation for merged replicates. Tracks were visualized in the WashU genome browser with replicates overlaid using the Matplot feature. Additional file 2: Table S3C shows the sequences used for demultiplexing 4C libraries as well as basic statistics for read mapping and quality control. 4C-seq data were also analyzed by using the FourCSeq Bioconductor package [51] to identify interacting regions for each viewpoint. Interacting regions were identified using the ‘addpeak’ command and are reported in Additional file 4: Table S4 (FDR < 0.05 and Z score > 3). Significance of annotated interacting regions shown in individual figures was calculated using the “subjecthits” command (4C-seq signal peak ± 2.5 kb).

Impact of cohesin depletion sex-biased gene expression

RNA-seq data for wild-type and cohesin-depleted (Nipbl-deficient) mouse liver (GSE93431) [21] was analyzed for n = 4 wild-type (control) and n = 4 Nipbl-depleted male mouse liver replicates. Data was RPKM-normalized, and reads were expressed relative to the mean of the wild-type group, which was set = 1 on a per gene basis by dividing the expression value for each individual replicate by the mean of the wild-type group plus a small pseudo count (1e−6) to avoid dividing by zero. Data is presented as mean relative expression ± SD for all plots. To determine the global effects of cohesin depletion on male-biased genes, the set of all expressed, strongly male-biased genes (FPKM > 1, and male/female (M/F) expression ratio > 3) was divided into two groups based on distance from their TSS to the nearest male-biased DHS or male-biased H3K27ac-marked region. Male-biased genes with proximal sex-biased enhancers (n = 29) were defined as having their TSS < 20 kb from the nearest male-biased DHS and from the nearest male-biased H3K27ac peak; and male-biased genes with distal sex-biased enhancers (n = 32) were those with TSS > 20 kb from both such regions. The underlying expression values for all genes are provided in Additional file 5: Table S2A, B.

Analysis of other datasets

Pan-cancer human expression data for Nox4 was obtained for tissues with matched Normal and Tumor expression values and analyzed using the web interface TIMER using default parameters (https://cistrome.shinyapps.io/timer/) [52]. Nox4 expression for male C3H mouse liver and neoplasms was obtained from [53] as processed, normalized RNA-seq counts (Accession # E-MTAB-6972). The following male and female mouse liver datasets were used to generate browser screenshots: DNase-seq [13], male liver CTCF and cohesin ChIP-seq data [31], sex-biased lncRNAs [9], M/F expression ratios [54], and H3K27ac ChIP-seq data [55]. Unless otherwise noted, FPKM and M/F expression ratios for protein coding genes are based on ribosomal RNA-depleted total liver RNA, while lncRNA expression ratios are based on ribosomal RNA-depleted liver nuclear RNA. Bigwig files were RIPM-normalized separately for each experiment type as described above. Genomic regions with significant sex bias are marked by short horizontal bars below each DHS and H3K27ac BigWig signal track. The bars were colored blue to indicate male-biased DHS or H3K27ac regions, and pink to indicate female-biased DHS or H3K27ac regions. A darker shade of color was used to indicate the stringency for the feature: DHS were shaded from dark to light to indicate high, standard, and low stringency for male-biased (blue) and female-biased (pink) DHS, as defined previously [13]. The high and low stringency designations used in that earlier study approximate the stringent and lenient definitions used here. For the sex-biased H3K27ac tracks, a total of four non-overlapping groups were defined based on the magnitude of the fold-change difference in H3K27ac peaks between male and female liver samples (M/F or F/M > 2.5 was defined as strict; and M/F or F/M > 1.5 but < 2.5 were defined as lenient). Both strict and lenient sex-biased H3K27ac peaks were defined using MAnorm with cutoffs of p-adj < 0.05 and read count > 15 for the upregulated peak [48]. These cutoffs resulted in 1583 female-biased H3K27ac peaks (380 strict, 1203 lenient) and 2241 male-biased H3K27ac peaks (604 strict, 1637 lenient).

Statistical analysis

Boxplots, cumulative distribution plots, and statistical analyses were implemented using GraphPad Prism 7. All boxplots are displayed using the Tukey convention, where interquartile range (IQR) was calculated as the difference between the 25th and 75th percentile. Outliers were considered as values falling above the 75th percentile value + 1.5*IQR, or below the 25th percentile value − 1.5*IQR. Whiskers indicate the maximum and minimum values in a set (if no outliers) or from the 75th percentile value + 1.5*IQR to the 25th percentile value − 1.5*IQR (if outliers exist). Boxes indicate the IQR with a horizontal line indicating the median value. Bar graphs show the mean and standard deviation. Unless otherwise indicated, all pairwise comparisons used two-tailed, nonparametric t tests. Comparisons of distributions were performed using a Kolmogorov–Smirnov (KS) test, while comparisons of values used a Mann–Whitney (M-W) test. The results are annotated in individual figures as follows: **** indicates p ≤ 0.0001; *** p ≤ 0.001; ** p ≤ 0.01; * p ≤ 0.05; and not significant (ns) for p > 0.05.

Data availability

CTCF and cohesin (Rad21) ChIP-seq data are available under accession numbers GSE130908 for female liver (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE130908) and GSE102997 for male liver (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE102997). 4C-seq data are available under accession number GSE130911 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE130911).

Results

Sex differences in CTCF and cohesin binding in mouse liver

We used ChIP-seq to identify binding sites for CTCF and the cohesin complex protein Rad21 in both male and female mouse liver. We observed significant sex-difference in factor binding (> 2-fold at FDR < 0.05) at 975 CTCF binding sites and at 1011 cohesin binding sites (Fig. 1a; sites are listed in Additional file 3: Table S1A and B). We applied three criteria to insure the robustness of each sex-differential site: fold-change > 2 and consistency across biological replicates (Additional file 1: Figure S1A–D) at FDR < 0.05 using a negative binomial model implemented in diffReps [47]; differential regions must overlap a genomic region called as a peak in at least one sample; and a minimum of 10 sequence reads in each peak region. Regions showing significant sex differential binding for both CTCF and cohesin represent 137 of all sex-differential sites (Fig. 1b; also see Additional file 1: Figure S1E for alternative filters to define sex-differential binding). An overall trend of sex-biased binding by both factors was seen for an even larger fraction of the sex-differential ChIP-seq peaks, as indicated by aggregate plots and heat maps (Fig. 1c).

Next, we classified the sex-differential CTCF and cohesin-binding sites, as follows: sites where both factors are bound (CAC sites) and show sex-differential binding of either cohesin (∆Coh) or CTCF (∆CTCF), or both factors; cohesin-only binding sites (CNC sites) with sex-differential cohesin binding; and CTCF-only sites (Lone CTCF sites) with sex-differential CTCF binding. CAC sites comprised 45–66% of all sex-differential sites (Additional file 1: Figure S1F) and generally showed stronger factor binding than the sex-differential CNC and Lone CTCF sites (Fig. 2a, Additional file 1: Figure S2A). The strength of factor binding (Fig. 2a), the CTCF motif score (Additional file 1: Figure S2B), and the percentage of sites with a CTCF motif (Additional file 1: Figure S2C) were generally higher for the female-biased sites than the male-biased sites. In contrast, a higher fraction of male-biased than female-biased Lone CTCF sites contained a CTCF motif (66% vs 48%, Additional file 1: Figure S2C), but there was no significant sex difference in normalized ChIP signal or motif score (Fig. 2a, Additional file 1: Figure S2C). The latter sex differences may be driven by additional factors, such as the inhibitory effect of DNA methylation on CTCF binding [56, 57], where the same sequence motif in male and female liver could be preferentially bound in males due to the hypermethylation of DNA seen in female compared to male mouse liver [58]. A majority (51.7%) of sex-biased CTCF peaks with CTCF motifs contain at least one CpG within the core CTCF motif, and therefore could be subject to regulation via sex differences in DNA methylation (Additional file 1: Figure S2D).

Sex-biased binding sites mapping to TAD and intra-TAD boundaries, genes and regulatory elements

We investigated whether the sex-differential binding of CTCF and cohesin is associated with sex-differential segmentation of the genome at the level of TAD and intra-TAD loops (DNA loop domains). TAD positions were derived from low resolution Hi-C experiments performed in mouse liver [34], while intra-TAD loops were based on a validated computational approach that identifies loop positions within TADs [31] of the 137 CAC peaks with significant sex differences in both CTCF and cohesin binding (Fig. 1b), 53 are on autosomes (Additional file 3: Table S1C). 17 of the 137 sex-differential CAC peaks overlap a TAD or intra-TAD loop anchor [31] in either male or female liver, of which 9 are on autosomes (Additional file 3: Table S1C). Ten of the 17 sites are associated with an intra-TAD loop predicted to be present in one sex only, of which 5 contained one of more sex-biased genes. One of these sex-based intra-TAD loop domains contained 6 sex-biased genes from the Cyp2c gene family and 2 sex-biased lncRNA genes (Additional file 1: Figure S2E). Consistent with the low frequency of sex-biased CAC sites at TAD or intra-TAD loop anchors, 88–93% of intra-TAD loops and anchors predicted in male liver were also predicted in female liver (Additional file 1: Figure S2F, Additional file 3: Table S1J). We conclude that CAC-mediated intra-TAD loops are largely conserved between the sexes in mouse liver, and thus do not play a major role in regulating sex-biased gene expression. Rather, TAD and intra-TAD loop domains are likely to indirectly guide enhancer–promoter looping in a sex-biased manner, as described below.

Next, we considered whether the sex-biased binding sites for CTCF and cohesin might help link sex-biased regulatory elements to sex-biased gene promoters. We examined the locations of these sites relative to enhancer DHS, defined as open chromatin regions (i.e., DHS) with a high ratio of H3K4me1 to H3K4me3 ChIP-seq signal [31]. We found that sex-biased CNC sites are much closer to enhancer DHS [median distance of 202 bp (in males) or 119 bp (in females)] than were sex-biased CAC sites (median distance = 8.8 kb (male-biased CAC) and 17.5 kb (female-biased CAC) (Fig. 2b, left). Thus, although sex-biased CNC sites are weaker binding than CAC sites (Fig. 2a), a majority are found at enhancers and may be functional. Of note, female-biased CTCF and cohesin binding sites tended to be closer to TSS than male-biased sites (Fig. 2c), despite equivalent sample quality between the sexes (17% reads in peaks for both male and female samples; Additional file 2: Table S3B). The factors underlying this apparent difference in TSS proximity are unclear.

All classes of sex-biased CTCF and cohesin sites tend to be distal from sex-biased genes: < 20% are within 20 kb of a sex-biased gene, and only 35 to 53% are within the same TAD and could therefore be considered potential cis regulators of sex-biased genes (Fig. 2d). Consistent with the association of cohesin with enhancers [25, 59], 25–29% of sex-biased CNC sites overlap a sex-biased enhancer (either sex-biased DHS or sex-biased H3K27ac peaks), as compared to only 7–11% of sex-biased CACs and Lone CTCF sites (Fig. 2e). Overall, however, a majority of all classes of sex-biased CTCF and cohesin sites are quite distant from sex-biased regulatory elements (Fig. 2e), consistent with their being distal regulators of sex-biased gene expression.

Sex-biased liver CNC sites and Lone CTCF sites showed much more tissue-specific binding of CTCF across mouse tissues [60, 61] than did sex-biased liver CAC sites (Additional file 1: Figure S3A, lower vs upper panels). Significant differences in the tissue-specificities of CTCF binding were also seen at male-biased compared to female-biased CAC sites and Lone CTCF sites (Additional file 1: Figure S3A, SB). In addition, liver-expressed genes that mapped to sex-biased CNC sites showed a more liver-specific expression pattern than genes mapping to sex-biased CAC sites, or liver-expressed genes overall (Additional file 1: Figure S3C). This suggests that sex-biased CNCs participate in tissue-specific transcriptional regulation, as was described for CNC peaks generally [25, 28, 29].

Impact of cohesin depletion on distally regulated male-biased genes

A substantial fraction (35–53%) of sex-biased CTCF and cohesin binding sites are found within the same TAD as at least one sex-biased gene (Fig. 2d) and could play a role in DNA looping between sex-biased enhancers and sex-biased gene promoters. Examples of sex-biased CTCF and/or cohesin binding sites that were either proximal (< 20 kb) or distal to sex-biased genes are shown in Fig. 3a–c. Figure 3a shows a female-biased enhancer (F-biased DHS and F-biased H3K27ac mark) with an overlapping female-biased CAC site (green arrows) located 33 kb upstream of the female-biased gene Slc22a29 (F/M ratio = 8.7), the closest TSS. Figure 3b shows Cml5, a male-specific gene (M/F expression ratio = 20.2) with a male-biased CNC site that overlaps a male-biased DHS ~ 3 kb upstream of its TSS (Fig. 3b, red arrow). The adjacent gene, Nat8 (M/F = 4.2), has a male-biased CAC site that overlaps a male-biased DHS positioned ~ 12 kb upstream of the TSS (Fig. 3b, green arrow). Conceivably, the sex-biased binding of cohesin and CTCF at these sites could contribute to looping of the associated sex-biased enhancers to their correspondingly sex-biased gene targets. Finally, Fig. 3c shows two male-biased complement C8 genes [C8a (M/F = 3.3) and C8b (M/F = 2.8)] that are linearly quite distant (~ 1.5 Mb) from a cluster of strongly male-biased enhancers near the 5′ end of the same TAD. The known TAD structure of this genomic region suggests these enhancers are spatially more proximal to the C8 genes than they are to than the linearly much closer Oma1 gene, located just inside the adjacent TAD (also see Additional file 1: Figure S4A).

Loss of cohesin binding in male mouse liver, achieved by depletion of the cohesin loading factor, Nipbl, leads to a loss of distal enhancer–promoter contacts and an increase in local ectopic contacts, which can activate proximal genes [21]. Using this public RNA-seq data for cohesin-depleted male mouse liver, we compared the effects of cohesin loss on the expression of male-biased genes with proximal sex-biased enhancers versus those that have only distal (> 20 kb) sex-biased enhancers. Figure 3d shows the relative changes in expression in cohesin-depleted compared to wild-type male mouse liver for the three genes shown in Fig. 3c. Expression of C8a and C8b decreased, by 98% and 82%, respectively, upon loss of chromatin-bound cohesin in male liver, while expression of Oma1 increased modestly (+ 22%), perhaps by an enhancer hijacking mechanism [62]. Further work in this mouse model is needed to characterize the effects of cohesin depletion in female mouse liver, and specifically to determine if Oma1 expression becomes sex biased due to de novo enhancer–promoter interactions. In contrast, the male-biased genes with proximal sex-biased enhancers, Cml5 and Nat8 (Fig. 3a), showed no significant change in expression following cohesin loss (Additional file 1: Figure S4B). We verified this requirement of cohesin for expression of distally regulated but not proximally regulated male-biased genes. Thus, male-biased genes with distal (> 20 kb) sex-biased regulatory elements were significantly more sensitive to loss of cohesin than male-biased genes with proximal sex-biased enhancers (median decrease in expression upon cohesin loss: 14.3- vs. 1.4-fold; Fig. 3e). This finding likely results from a requirement for cohesin for distal interactions, via either a direct or an indirect looping mechanism. Conceivably, for sex-biased genes with nearby sex-biased regulatory elements, enhancer–promoter loops required for gene expression can be maintained over short genomic distances by transcription factors such as Mediator [25] or YY1 [63], and without a need for cohesin.

4C-seq analysis of sex-biased chromatin interactions in mouse liver

We performed 4C-seq analysis centered on six viewpoints at sex-biased enhancers in five distinct genomic regions to determine whether a sex-bias in enhancer–promoter loops is associated with sex-biased gene expression in mouse liver. Our findings (Fig. 4) are based on n = 3 individual biological replicates per sex, whose 4C-seq interaction profiles are displayed in Additional file 1: Figure S6. We also indicate locations of computationally predicted CAC-mediated intra-TAD loops (purple arcs) in mouse liver [31], which are similar in size (median 151 kb) and abundance (9543 loops total) to loop domains [32] and Insulated Neighborhoods [26, 64].

A1bg region

We anchored a 4C-seq viewpoint at a female-biased enhancer near the strongly female-biased gene A1bg (F/M = 155) (Fig. 4a, vertical pink bar; Additional file 1: Figure S6A). Nearby are 12 female-biased, mono-exonic nuclear-enriched lncRNAs [10] in three clusters across the genomic region displayed. Robust interactions were observed in female but not male livers between the viewpoint enhancer and three genomic regions (red arrows): a strong female-biased enhancer (right arrow), the promoter of A1bg (middle arrow), and a region downstream of A1bg that contains a cluster of four female-specific lncRNAs (lncRNAs 12590-93; left arrow), where we observed the strongest interactions. The lncRNAs in this cluster are more highly expressed (Additional file 1: Figure S5C) and are more consistently female-biased across various RNA-seq datasets than the other two lncRNA clusters (Additional file 1: Figure S5D). The maximum expression of these 12 lncRNAs ranged from 0.31 to 2.71 FPKM in female liver compared to 0 to 0.02 FPKM in male liver (Additional file 1: Figure S5C).

The precise relationship between the female-biased expression of these lncRNAs and the female-bias in 3D interactions with the distal enhancer is not known. The interaction may be regulatory in nature (e.g., an enhancer–promoter interaction, as with any gene) or it could be facilitated by one or more of the 12 nuclear-enriched, female-biased lncRNAs, as was described for the lncRNAs Xist [65], Firre [66], and Haunt [67]. Alternatively, the female-specific interactions shown may be primarily those of regulatory enhancers driving expression of several female-specific genes—including A1bg and multiple lncRNAs. The female-biased CTCF binding seen at both interacting regions (right and left arrows) lends mechanistic support for the latter proposal, with CTCF mediating enhancer–promoter and enhancer–enhancer interactions. As CTCF is known to interact with lncRNAs in a functional manner, and with high affinity [68], these two mechanisms are not mutually exclusive; one or more of these highly female-specific lncRNAs (Additional file 1: Figure S5D) could function in a cis-acting manner to selectively guide CTCF binding and interactions unique to female liver.

Gm4794 and Sult3a1 region

We used 4C-seq to interrogate an enhancer viewpoint proximal to the highly female-biased gene Gm4794 (F/M = 704; also known as Sult3a2) and also two female-biased mono-exonic lncRNAs, lnc8820 and lnc8821 (F/M = 34 and 90) (Fig. 4b, Additional file 1: Figure S6B). The enhancer viewpoint is distal to two other female-biased protein coding genes, Sult3a1 (F/M = 288) and Rsph4a (F/M = 34). We observed female-biased 4C-seq interactions with the proximal promoter of Gm4794 (left red arrow) and with strong female-biased enhancers overlapping lnc8820 and lnc8821. Unlike the enhancer neighboring A1bg, these interactions are associated with female-biased CNC peaks present at both the viewpoint enhancer and the downstream interacting enhancer. Despite the absence of CAC insulator elements across this genomic region (and consequently, the absence of intra-TAD loops), all focal interactions are local, within ~ 35 kb of the viewpoint. However, the viewpoint enhancer also made weak interactions to a broad, ~ 45 kb region extending from ~ 15 kb upstream of the female-biased promoter of Sult3a1 to ~ 10 kb beyond the gene body (Fig. 4b, red bracket). This 45-kb region contains several robust female-biased DHS and H3K27ac peaks, but lacks cohesin binding, which may account for the lack of strong, focal interactions with the viewpoint enhancer. Gm4794 and Sult3a1 may both interact with the enhancer that overlaps lnc8820 and lnc8821, but this weak (and perhaps indirect) association cannot be captured by proximity ligation under standard 4C-seq conditions. Reciprocal 4C-seq experiments anchored at these ncRNA TSS or alternative 4C methods with increased sensitivity [69] may be needed to validate these weaker interactions.

C9 region

We examined a 4C-seq promoter viewpoint placed at the complement factor C9 gene (M/F = 3.5) to investigate whether intra-TAD loops can indirectly coordinate enhancer–promoter contacts preferential to one sex in a genomic region without sex-biased CTCF or cohesin binding. C9 overlaps an antisense lncRNA that shows a sevenfold greater male-bias in expression (lnc12340; M/F = 26) (Fig. 4c, Additional file 1: Figure S6C). A strongly male-biased enhancer region lies ~ 230 kb upstream of the TSS of C9, and is characterized by male-biased DHS and H3K27ac peaks, whereas the TSS of C9 has only a male-biased DHS. The far upstream enhancer and the TSS of C9 both fall just outside of (< 10 kb from) a nested pair of intra-TAD loops that encompass both TSS of Dab2 (Fig. 4c, bottom track) and at least partially insulate Dab2, whose expression in male liver is 87-fold lower than C9 (FPKM = 0.7 vs 61). The promoter region of C9 interacts with the cluster of far upstream enhancers (Fig. 4c, red arrow), with stronger interactions seen in male liver. Weaker, mostly non-focal interactions were seen between the C9 promoter and several sites within the nested intra-TAD loops. This insulation-by-looping mechanism allows the strong male-biased enhancer to bypass a more proximal gene, Dab2, to drive expression of C9; it also helps explain the 87-fold higher expression of C9 compared to Dab2 in male liver. The shorter isoform of Dab2 shows weak, male-specific interaction, despite a lack of sex-bias in its expression. Given the interactions with the strong, far upstream enhancer region, we hypothesized that the expression of C9 and the antisense lnc12340 would be sensitive to the loss cohesin of binding. Indeed, both genes showed a three to fivefold decrease in expression in cohesin-depleted mouse liver (p < 0.01 for both), while the insulated gene Dab2 showed no significant change in expression (Fig. 4d). While there was an apparent increase in Dab2 expression, that increase did not reach statistical significance (Fig. 4d). Of note, the loss of CAC loop insulation upon cohesin depletion would be expected to result in such an increase, which would allow the upstream enhancer to more freely interact with the promoters of both Dab2 isoforms. Together, these findings support a model whereby an active, male-biased enhancer preferentially interacts with the highly expressed male-biased gene C9, while bypassing the lowly expressed intervening gene Dab2 (also see Discussion and model in Fig. 6).

Nudt7 is a highly expressed, male-biased gene with distal male-biased enhancers within the same intra-TAD loop (M/F = 2.5; FPKM = 132 in male liver). Although there are some weak male-biased DHS within the gene body of Nudt7 (Fig. 4e), there is no apparent sex bias at its shared promoter with the sex-biased lncRNA gene lnc7430 (M/F = 4.0; FPKM = 2.9 in male liver). Approximately 22 kb and 39 kb upstream of the TSS of Nudt7 are two male-biased enhancers with male-biased DHS and H3K27ac marks; the latter also is proximal (~ 2.5 kb upstream) to the TSS of the male-biased lnc7423 (M/F = 7.0; FPKM = 2.8 in male liver). This enhancer cluster is one of 503 super-enhancers found in both male and female liver [31]. We observed male-biased 4C-seq interactions between the enhancer viewpoint and a neighboring male-biased enhancer, and also with the Nudt7 and lnc7430 promoter(s) (Fig. 4e, red arrows; Additional file 1: Figure S6D). In contrast, we did not observe focal interactions in either sex between the enhancer viewpoint and a sex-independent enhancer 15.7 kb upstream of Nudt7. Both Nudt7 and lnc7430 were strongly down regulated in male liver upon cohesin depletion (Fig. 4f), suggesting their expression is dependent on interactions facilitated by the intra-TAD loop encompassing this genomic region. Expression of lnc7423 was not significantly reduced, perhaps due to its closer proximity to strong male-biased enhancers.

Sex-independent, nested intra-TAD loops restrict Nox4 to proximal enhancer–promoter interactions

NADPH oxidase 4 (Nox4) exhibits male-biased expression in mouse liver (M/F = 7.7) and may contribute to a number of liver pathologies whose incidence or severity is male-biased [70, 71]. Nox4 is highly upregulated in tumor compared to healthy liver tissue of mice that spontaneously develop liver tumors (Additional file 1: Figure S7A), and in humans, Nox4 is upregulated in hepatocellular carcinoma and other cancers (Additional file 1: Figure S7B). Mouse Nox4 is located within a pair of nested intra-TAD loops (Fig. 5a, bottom), and has a strong male-biased enhancer 11.5 kb upstream of the TSS and a strong male-biased DHS 125 kb downstream of the TSS. However, only the upstream region has H3K27ac (active enhancer) marks (Fig. 5a). We placed 4C-seq viewpoints at both the upstream region (viewpoint VP1, at − 11.5 kb; red vertical highlight in Fig. 5a) and the downstream region (viewpoint VP2, at + 125 kb; green vertical highlight) to investigate chromatin interactions with each putative regulatory region. Interactions with the downstream DHS at VP2 were limited to the domain defined by the pair of 3′ anchors of the nested intra-TAD loops, consistent with these loops insulating from distal interactions (Fig. 5a, green bracket at bottom). As a result, VP2 did not interact with the promoter of Nox4 or with the -11.5 kb enhancer in either male or female liver. Furthermore, VP2 did not show any consistent male-biased interactions, despite its location at a strong male-biased DHS.

In contrast, the − 11.5 kb enhancer at VP1 showed male-biased interactions with several genomic regions, including the Nox4 promoter and two genomic regions that are in a transcribed-like chromatin state [14] in male liver, but are in an inactive state in female liver (Fig. 5a, regions Y1 and Y2; purple and red in chromatin state maps at top, respectively). Further, both VP1-interacting regions include a short sequence in an enhancer state (Fig. 5a, green band within purple region; see Additional file 1: Figure S7C). We also observed 4C-seq interactions between VP1 and a region just upstream of the shorter intra-TAD loop 3′ anchor 31 kb downstream of the Nox4 TSS in both male and female liver (Fig. 5a, red arrow; Additional file 1: Figure S7D). The high intra-domain interactions that we observed for VP1 and the shorter nested intra-TAD loop, and the low frequency of inter-domain interactions, are indicative of the insulation activity of these loops [31]. Immediately downstream of the 3′ loop anchor, in region Y2, we observed 4C-seq interactions specific to male liver, which may be due to increased movement of cohesin beyond this loop anchor as a result of the increased transcription of Nox4 in male liver. Movement of cohesin has been directly linked to transcription both in vitro and in vivo [18, 72, 73], and may result in a more dynamic loop structure and weaker local insulation at the 3′ loop anchor of the shorter nested loop. The VP1-interacting regions Y1 and Y2 are both in a transcribed chromatin state only in male liver, consistent with the male-biased 4C-seq interactions of VP1 with both regions.

These findings support the predicted model of two nested intra-TAD loops, with the smaller enclosed loop insulated from the larger enclosing loop. Domain predictions for other mouse tissues, based on computational methods and experimentally observed looping in mouse embryonic stem cells [26, 60], support the conclusion that the genomic regions defined by VP1 and VP2 are in separate domains (Fig. 5a, bottom). Accordingly, only the -11.5 kb enhancer at VP1 would be predicted to interact with the Nox4 promoter. Generally, the cis regulatory elements relevant for the regulation of Nox4 appear to be contained within the smaller intra-TAD loop. It is less clear what regulatory function the male-biased DHS at VP2 plays, as it does not interact with Nox4 or with the downstream gene, Tyr, which is not expressed in liver (FPKM < 0.01 in both sexes).

Discussion

We investigated sex differences in autosomal 3D genome organization in the mouse liver model, focusing on sex-based differences in chromatin binding and interactions involving cohesin and CTCF, which mediate long-range DNA looping interactions that segment mammalian genomes into megabase-scale TAD domains and their shorter intra-TAD domains. We identified 1847 binding sites for cohesin and/or CTCF that show significant differential occupancy between male and female mouse liver; however, very few of these sites were associated with sex differences in TAD or intra-TAD loop anchors. Given the low resolution of Hi-C datasets available for mouse liver [34], it is possible that we are underestimating the number of TADs present, which limits our interpretation of TAD loop anchors to only a subset of the true total. Furthermore, high resolution Hi-C would also increase our confidence in determining which specific CTCF sites anchor long-range loops. A majority of the sex-biased binding sites classified as cohesin-non-CTCF (CNC) sites (but only a minority of cohesin-and-CTCF (CAC) and Lone CTCF sites) mapped to distal enhancers, and a major subset of these overlapped sex-biased enhancers. This finding indicates a role for cohesin in sex-biased enhancer activity, and is consistent with the extensive sex bias in chromatin state seen at sex-biased enhancers but not promoters [14]. We also found that 77% (72 of 96) of sex-biased CNCs that overlap sex-biased enhancers are > 20 kb from a TSS of sex-biased gene (median distance 238 kb), consistent with the well-characterized role of cohesin in mediating distal enhancer–promoter interactions [25, 27]. Furthermore, male-biased genes with distal but not proximal sex-biased enhancers were much more sensitive to cohesin depletion than genes with proximal sex-biased enhancers, implicating cohesin in long-range enhancer interactions regulating these sex-biased genes. Finally, by applying circularized chromosome conformation capture with sequencing (4C-seq) to sex-biased enhancer viewpoints in five genomic regions, we established that sex differences in chromatin interactions are a common feature of sex-biased gene expression in the liver, and we elucidated how chromatin interactions link sex-biased genes to distal sex-biased enhancers, guided both directly and indirectly by cohesin and/or CTCF looping.

Although TADs and intra-TADs are largely conserved across tissues, 20-30% of all such CAC-mediated loops are cell type-specific [31]. Nevertheless, when comparing male and female mouse liver, which show extensive hormone-determined differences in epigenetic state [14, 54, 74], we did not find evidence for widespread formation of sex-specific intra-TAD loops. Rather, we found that intra-TAD loops in mouse liver are largely sex-independent and devoid of sex-biased CTCF or cohesin binding at their CAC anchors. These loops do have the ability, however, to indirectly facilitate sex-dependent chromatin interactions. Thus, a sex-independent intra-TAD loop was shown to insulate the super-enhancer-associated male-biased gene Nudt7, and nested intra-TAD loops insulating Nox4 restricted the promoter of this male-biased gene from an intronic enhancer while enabling interactions with a cluster of upstream enhancers. Although the currently available Hi-C data for mouse liver is low resolution, our pair of 4C-seq viewpoints in the Nox4 region support a nested loop model that is supported by computational predictions for liver and by SMC1 ChIA-PET experimental data in embryonic stem cells. Furthermore, our analysis of female-biased gene regions revealed female-biased proximal enhancer–promoter interactions in the Sult3a gene region associated with female-biased cohesin binding, as well as female-biased interactions between the A1bg promoter, a far distal (> 100 kb) enhancer, and distal female-biased CTCF binding sites. Together, these findings support the proposal that CTCF and cohesin contribute in both direct and indirect ways to the formation of sex-biased enhancer–promoter contacts in mouse liver (Fig. 6).

Our analysis of the male-biased complement factor gene C9 provides an interesting example of sex-independent CAC-looping that indirectly facilitates sex-biased enhancer–promoter contacts. C9 interacts strongly with a distal (> 200 kb upstream) male-biased enhancer while bypassing the weakly expressed (and sex-independent) Dab2 gene region, which is insulated by a nested pair of intra-TAD loops. These nested loops, in turn, bring the TSS of C9 into much closer proximity to a cluster of far upstream male-biased enhancers than would be achieved based on linear genomic distance alone (Fig. 6b). Furthermore, we observed more frequent contacts between C9 and the far upstream enhancers in male compared to female liver, despite the absence of any male-biased binding of CTCF or cohesin to help explain sex differences in contact frequency. Conceivably, the male-biased DHS located in the upstream region facilitate direct enhancer–promoter interactions specifically in male liver by working in the context of other known looping factors. Mechanistically, these interactions could be driven by proteins, such as YY1 [63] or Mediator [25], or by non-protein factors such as enhancer-RNAs (eRNAs) [75, 76]. Of note, 3 of the 4 male-biased enhancers far upstream of C9 are actively transcribed to produce bidirectional eRNA transcripts in male mouse liver [77], and all four regions are bound by the protein YY1 [78] (Additional file 1: Figure S6E).

A majority of the 1847 sex-biased binding sites for CTCF and cohesin identified here are intergenic and distal to sex-biased genes and enhancers. Furthermore, many of these sex-biased binding sites are located in TADs without any known sex-biased genes. Although those sites lack any obvious link to sex-biased gene expression in untreated liver, they could have a priming effect and contribute to sex-specific responses reported for hepatic stressors, such as high fat diet [79] and xenobiotic exposure [80], which we recently showed can induce a sex-biased gene response in TADs whose genes do not show a sex-bias in expression in the unstressed state [81]. Just as short-term feeding of a high fat diet can leave a lasting epigenetic memory in the form of epigenetic modifications [82], sex-biased binding of CTCF and/or cohesin may differentially prime each sex for distinct looping patterns that enable the observed sex-biased responses to chemical exposure or dietary stressors.

While CAC sites at TAD and intra-TAD boundaries have a well-established role as anchors that enable loop domain-level nuclear organization [17, 26, 27, 31, 32], non-anchor CAC sites may directly link enhancers to promoters or to other enhancers, and thereby contribute to interactions governing tissue-specific gene expression [35, 40]. A majority of sex-differential CAC binding occurs at non-anchor CAC sites (Additional file 3: Table S1), a subset of which may mediate long-distance interactions involving sex-biased enhancers and gene promoters. Specific examples described here include the enhancer–enhancer and enhancer–promoter contacts that we identified by 4C-seq for A1bg in the context of female-biased CTCF binding. Similarly, more than half of male-biased liver CTCF binding occurred at Lone CTCF sites, which we found are closer than CAC sites to gene TSS, and can also play a non-canonical role in looping between enhancers and promoters [35].

4C-seq interactions between sex-biased enhancer viewpoints and distal sex-biased lncRNAs [9, 10] were found in three of the five sex-biased genomic regions we investigated. In one example, the highly female-biased gene A1bg is nearby several strongly female-biased, nuclear-enriched mono-exonic lncRNAs, several of which are transcribed from genomic loci that show female-specific interactions with the distal female-specific enhancer viewpoint that we interrogated. Enhancer-associated lncRNAs have been defined as intergenic transcripts with enhancer chromatin marks whose expression is tissue-restricted and is associated with increased expression of nearby expressed protein coding genes [83, 84]. Based on our findings, we suggest that functional enhancer-associated lncRNAs might be identified by their looping interactions with enhancer sequences, which can be determined globally using high throughput interaction methods, such as Hi-C [85].

The sex-biased cohesin and CTCF binding sites described here were discovered using livers from adult (8 week) mice, and likely encompass only a subset of all sex-differential cohesin and CTCF binding sites across the lifespan of a mouse, given the dramatic changes in sex-biased gene expression that occur during prenatal and especially postnatal liver development [86, 87]. The sex-biased binding of cohesin and CTCF to liver chromatin in adult liver is expected to be regulated by pituitary growth hormone secretion, which is sex-dependent and produces the sex-dependent plasma growth hormone profiles that regulate the vast majority of sex differences in the adult liver, including differences in gene expression [74, 88], transcription factor binding [89, 90], and chromatin states [14, 54]. Given that CTCF binding to DNA can be inhibited by DNA methylation [56, 57, 91], the hypomethylation of enhancer sequences seen in male compared to female liver [58] could contribute to male-specific CTCF binding at such sites. Such an effect is expected to become more pronounced with age, given the increased male hypomethylation reported in older mice [58, 92].

The sex-specific patterns of pituitary growth hormone secretion regulating sex-biased liver gene expression emerge at puberty, and have been implicated in the dynamic regulation of liver chromatin states in both male and female adult mouse liver [54, 74]. We do not know when during mouse liver development the sex-differential chromatin interactions described here are first established, whether they constitute a relatively fixed (static) 3D framework governing transcription in male and female nuclei, or alternatively, whether they respond dynamically to the temporal changes in plasma growth hormone profiles that regulate sex differences in liver chromatin states. The potential for dynamic, reversible changes in DNA looping was demonstrated in a study where auxin-induced cohesin degradation led to a loss of virtually all DNA loops after 40 min, followed by their reestablishment within one hour of auxin withdrawal and cohesin reintroduction [22]. Extensive changes in chromatin interactions also occur during circadian oscillations in mouse liver [93,94,95] and during macrophage differentiation [96]. Our finding that the expression of distally regulated male-biased genes is highly sensitive to cohesin depletion illustrates the importance of distal enhancer–promoter interactions in maintaining male-biased gene expression. Of note, however, this response is largely not a direct result of a loss of sex differential cohesin and/or CTCF binding at distal male-biased enhancers, as only 5 of 32 distally regulated male-biased genes have male-biased CTCF/cohesin binding at their distal male-biased enhancers. Alternatively, the sex-biased enhancer–promoter and enhancer–enhancer loops that we describe here might be determined by intrinsic sex differences [88], or might be established by early postnatal hormone exposures that program liver gene expression [97]. Studies to detect dynamic changes in DNA looping and quantify changes in chromatin interaction strength, e.g., during the course of a male plasma growth hormone pulse [74], will likely require improvements to the 4C-seq protocol, including the use of unique molecular identifiers for more accurate quantification of interactions [69] and the elimination of any PCR artifacts associated with over-amplification, which are difficult to address using conventional 4C-seq methods [45] and may have decreased the magnitude of the apparent sex differences in chromatin interactions seen in our work.

In conclusion, we employed the mouse liver model with its extensive sex differences in gene expression to study sex differences in nuclear organization and DNA looping interactions in a non-reproductive tissue exposed to sex-unique patterns of hormone stimulation. We determined that male-biased genes with distal but not proximal sex-biased enhancers are particularly sensitive to the loss of cohesin binding. Furthermore, while most sex-biased binding sites for CTCF and cohesin were found to be distal from sex-biased genes, a subset likely contributes to sex-biased looping between regulatory elements in cis, as exemplified by the female-biased DNA looping interactions observed for A1bg. In addition, sex-independent CAC-looping may indirectly provide sex specificity to chromatin interactions by insulating male-biased genes such as Nudt7, or by bringing a sex-biased gene into closer proximity to a cluster of sex-biased enhancers, as demonstrated for C9. Together, these findings illustrate the direct and indirect contributions that cohesin and CTCF can make to sex-biased gene expression in the liver, and may be broadly applicable to other biological systems where distal regulation of gene expression is of interest.

Availability of data and materials

The datasets generated and/or analyzed during the current study are included in this published article and its supplementary information files. Raw and processed sequencing files generated in this study are available at GEO under accession numbers GSE130908, GSE102997, and GSE130911.

References

Mayne BT, Bianco-Miotto T, Buckberry S, Breen J, Clifton V, Shoubridge C, Roberts CT. Large scale gene expression meta-analysis reveals tissue-specific, sex-biased gene expression in humans. Front Genet. 2016;7:183.
Article PubMed PubMed Central CAS Google Scholar
Qu K, Zaba LC, Giresi PG, Li R, Longmire M, Kim YH, Greenleaf WJ, Chang HY. Individuality and variation of personal regulomes in primary human T cells. Cell Syst. 2015;1(1):51–61.
Article CAS PubMed PubMed Central Google Scholar
Rinn JL, Rozowsky JS, Laurenzi IJ, Petersen PH, Zou K, Zhong W, Gerstein M, Snyder M. Major molecular differences between mammalian sexes are involved in drug metabolism and renal function. Dev Cell. 2004;6(6):791–800.
Article CAS PubMed Google Scholar
Zhang Y, Klein K, Sugathan A, Nassery N, Dombkowski A, Zanger UM, Waxman DJ. Transcriptional profiling of human liver identifies sex-biased genes associated with polygenic dyslipidemia and coronary artery disease. PLoS ONE. 2011;6(8):e23506.
Article CAS PubMed PubMed Central Google Scholar
Ruggieri A, Barbati C, Malorni W. Cellular and molecular mechanisms involved in hepatocellular carcinoma gender disparity. Int J Cancer. 2010;127(3):499–504.
Article CAS PubMed Google Scholar
Guy J, Peters MG. Liver disease in women: the influence of gender on epidemiology, natural history, and patient outcomes. Gastroenterol Hepatol (N Y). 2013;9(10):633–9.
Google Scholar
Buzzetti E, Parikh PM, Gerussi A, Tsochatzis E. Gender differences in liver disease and the drug-dose gender gap. Pharmacol Res. 2017;120:97–108.
Article CAS PubMed Google Scholar
Clodfelter KH, Holloway MG, Hodor P, Park SH, Ray WJ, Waxman DJ. Sex-dependent liver gene expression is extensive and largely dependent upon signal transducer and activator of transcription 5b (STAT5b): STAT5b-dependent activation of male genes and repression of female genes revealed by microarray analysis. Mol Endocrinol. 2006;20(6):1333–51.
Article CAS PubMed Google Scholar
Melia T, Hao P, Yilmaz F, Waxman DJ. Hepatic long intergenic noncoding RNAs: high promoter conservation and dynamic, sex-dependent transcriptional regulation by growth hormone. Mol Cell Biol. 2016;36(1):50–69.
Article CAS PubMed Google Scholar
Melia T, Waxman DJ. Sex-biased lncRNAs inversely correlate with sex-opposite gene coexpression networks in diversity outbred mouse liver. Endocrinology. 2019;160(5):989–1007.
Article CAS PubMed PubMed Central Google Scholar
Hao P, Waxman DJ. Functional roles of sex-biased, growth hormone-regulated MicroRNAs miR-1948 and miR-802 in young adult mouse liver. Endocrinology. 2018;159(3):1377–92.
Article CAS PubMed PubMed Central Google Scholar
Waxman DJ, Holloway MG. Sex differences in the expression of hepatic drug metabolizing enzymes. Mol Pharmacol. 2009;76(2):215–28.
Article CAS PubMed PubMed Central Google Scholar
Ling G, Sugathan A, Mazor T, Fraenkel E, Waxman DJ. Unbiased, genome-wide in vivo mapping of transcriptional regulatory elements reveals sex differences in chromatin structure associated with sex-specific liver gene expression. Mol Cell Biol. 2010;30(23):5531–44.
Article CAS PubMed PubMed Central Google Scholar
Sugathan A, Waxman DJ. Genome-wide analysis of chromatin states reveals distinct mechanisms of sex-dependent gene regulation in male and female mouse liver. Mol Cell Biol. 2013;33(18):3594–610.
Article CAS PubMed PubMed Central Google Scholar
Phillips JE, Corces VG. CTCF: master weaver of the genome. Cell. 2009;137(7):1194–211.
Article PubMed PubMed Central Google Scholar
Ghirlando R, Felsenfeld G. CTCF: making the right connections. Genes Dev. 2016;30(8):881–91.
Article CAS PubMed PubMed Central Google Scholar
Sanborn AL, Rao SS, Huang SC, Durand NC, Huntley MH, Jewett AI, Bochkov ID, Chinnappan D, Cutkosky A, Li J, et al. Chromatin extrusion explains key features of loop and domain formation in wild-type and engineered genomes. Proc Natl Acad Sci USA. 2015;112(47):E6456–65.
Article CAS PubMed PubMed Central Google Scholar
Davidson IF, Goetz D, Zaczek MP, Molodtsov MI, Huis In ‘t Veld PJ, Weissmann F, Litos G, Cisneros DA, Ocampo-Hafalla M, Ladurner R, et al. Rapid movement and transcriptional re-localization of human cohesin on DNA. EMBO J. 2016;35(24):2671–85.
Article CAS PubMed PubMed Central Google Scholar
Fedoriw AM, Stein P, Svoboda P, Schultz RM, Bartolomei MS. Transgenic RNAi reveals essential function for CTCF in H19 gene imprinting. Science. 2004;303(5655):238–40.
Article CAS PubMed Google Scholar
Xu H, Balakrishnan K, Malaterre J, Beasley M, Yan Y, Essers J, Appeldoorn E, Tomaszewski JM, Vazquez M, Verschoor S, et al. Rad21-cohesin haploinsufficiency impedes DNA repair and enhances gastrointestinal radiosensitivity in mice. PLoS ONE. 2010;5(8):e12112.
Article PubMed PubMed Central CAS Google Scholar
Schwarzer W, Abdennur N, Goloborodko A, Pekowska A, Fudenberg G, Loe-Mie Y, Fonseca NA, Huber W, Haering CH, Mirny L, et al. Two independent modes of chromatin organization revealed by cohesin removal. Nature. 2017;551(7678):51–6.
Article PubMed PubMed Central Google Scholar
Rao SSP, Huang SC, Glenn St Hilaire B, Engreitz JM, Perez EM, Kieffer-Kwon KR, Sanborn AL, Johnstone SE, Bascom GD, Bochkov ID, et al. Cohesin loss eliminates all loop domains. Cell. 2017;171(2):305–320.e324.
Article CAS PubMed PubMed Central Google Scholar
Nora EP, Goloborodko A, Valton AL, Gibcus JH, Uebersohn A, Abdennur N, Dekker J, Mirny LA, Bruneau BG. Targeted degradation of CTCF decouples local insulation of chromosome domains from genomic compartmentalization. Cell. 2017;169(5):930–944.e922.
Article CAS PubMed PubMed Central Google Scholar
Ciosk R, Shirayama M, Shevchenko A, Tanaka T, Toth A, Shevchenko A, Nasmyth K. Cohesin’s binding to chromosomes depends on a separate complex consisting of Scc2 and Scc4 proteins. Mol Cell. 2000;5(2):243–54.
Article CAS PubMed Google Scholar
Kagey MH, Newman JJ, Bilodeau S, Zhan Y, Orlando DA, van Berkum NL, Ebmeier CC, Goossens J, Rahl PB, Levine SS, et al. Mediator and cohesin connect gene expression and chromatin architecture. Nature. 2010;467(7314):430–5.
Article CAS PubMed PubMed Central Google Scholar
Dowen JM, Fan ZP, Hnisz D, Ren G, Abraham BJ, Zhang LN, Weintraub AS, Schuijers J, Lee TI, Zhao K, et al. Control of cell identity genes occurs in insulated neighborhoods in mammalian chromosomes. Cell. 2014;159(2):374–87.
Article CAS PubMed PubMed Central Google Scholar
Phillips-Cremins JE, Sauria ME, Sanyal A, Gerasimova TI, Lajoie BR, Bell JS, Ong CT, Hookway TA, Guo C, Sun Y, et al. Architectural protein subclasses shape 3D organization of genomes during lineage commitment. Cell. 2013;153(6):1281–95.
Article CAS PubMed PubMed Central Google Scholar
Faure AJ, Schmidt D, Watt S, Schwalie PC, Wilson MD, Xu H, Ramsay RG, Odom DT, Flicek P. Cohesin regulates tissue-specific expression by stabilizing highly occupied cis-regulatory modules. Genome Res. 2012;22(11):2163–75.
Article CAS PubMed PubMed Central Google Scholar
Schmidt D, Schwalie PC, Ross-Innes CS, Hurtado A, Brown GD, Carroll JS, Flicek P, Odom DT. A CTCF-independent role for cohesin in tissue-specific transcription. Genome Res. 2010;20(5):578–88.
Article CAS PubMed PubMed Central Google Scholar
Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, Hu M, Liu JS, Ren B. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485(7398):376–80.
Article CAS PubMed PubMed Central Google Scholar
Matthews BJ, Waxman DJ. Computational prediction of CTCF/cohesin-based intra-TAD loops that insulate chromatin contacts and gene expression in mouse liver. Elife. 2018;7:e34077.
Article PubMed PubMed Central Google Scholar
Rao SS, Huntley MH, Durand NC, Stamenova EK, Bochkov ID, Robinson JT, Sanborn AL, Machol I, Omer AD, Lander ES, et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell. 2014;159(7):1665–80.
Article CAS PubMed PubMed Central Google Scholar
Dixon JR, Jung I, Selvaraj S, Shen Y, Antosiewicz-Bourget JE, Lee AY, Ye Z, Kim A, Rajagopal N, Xie W, et al. Chromatin architecture reorganization during stem cell differentiation. Nature. 2015;518(7539):331–6.
Article CAS PubMed PubMed Central Google Scholar
Vietri Rudan M, Barrington C, Henderson S, Ernst C, Odom DT, Tanay A, Hadjur S. Comparative Hi-C reveals that CTCF underlies evolution of chromosomal domain architecture. Cell Rep. 2015;10(8):1297–309.
Article CAS PubMed PubMed Central Google Scholar
Ren G, Jin W, Cui K, Rodrigez J, Hu G, Zhang Z, Larson DR, Zhao K. CTCF-mediated enhancer–promoter interaction is a critical regulator of cell-to-cell variation of gene expression. Mol Cell. 2017;67(6):1049–1058.e1046.
Article CAS PubMed PubMed Central Google Scholar
Splinter E, Heath H, Kooren J, Palstra RJ, Klous P, Grosveld F, Galjart N, de Laat W. CTCF mediates long-range chromatin looping and local histone modification in the beta-globin locus. Genes Dev. 2006;20(17):2349–54.
Article CAS PubMed PubMed Central Google Scholar
Huang J, Li K, Cai W, Liu X, Zhang Y, Orkin SH, Xu J, Yuan GC. Dissecting super-enhancer hierarchy based on chromatin interactions. Nat Commun. 2018;9(1):943.
Article PubMed PubMed Central Google Scholar
Donohoe ME, Zhang LF, Xu N, Shi Y, Lee JT. Identification of a Ctcf cofactor, Yy1, for the X chromosome binary switch. Mol Cell. 2007;25(1):43–56.
Article CAS PubMed Google Scholar
Park HJ, Li J, Hannah R, Biddie S, Leal-Cervantes AI, Kirschner K, Flores Santa Cruz D, Sexl V, Gottgens B, Green AR. Cytokine-induced megakaryocytic differentiation is regulated by genome-wide loss of a uSTAT transcriptional program. EMBO J. 2016;35(6):580–94.
Article CAS PubMed Google Scholar
Guo Y, Xu Q, Canzio D, Shou J, Li J, Gorkin DU, Jung I, Wu H, Zhai Y, Tang Y, et al. CRISPR inversion of CTCF sites alters genome topology and enhancer/promoter function. Cell. 2015;162(4):900–10.
Article CAS PubMed PubMed Central Google Scholar
Oti M, Falck J, Huynen MA, Zhou H. CTCF-mediated chromatin loops enclose inducible gene regulatory domains. BMC Genomics. 2016;17:252.
Article PubMed PubMed Central CAS Google Scholar
Whalen S, Truty RM, Pollard KS. Enhancer–promoter interactions are encoded by complex genomic signatures on looping chromatin. Nat Genet. 2016;48(5):488.
Article CAS PubMed PubMed Central Google Scholar
Yang Y, Zhang R, Singh S, Ma J. Exploiting sequence-based features for predicting enhancer–promoter interactions. Bioinformatics. 2017;33(14):i252–60.
Article CAS PubMed PubMed Central Google Scholar
Dekker J, Rippe K, Dekker M, Kleckner N. Capturing chromosome conformation. Science. 2002;295(5558):1306–11.
Article CAS PubMed Google Scholar
van de Werken HJ, Landan G, Holwerda SJ, Hoichman M, Klous P, Chachik R, Splinter E, Valdes-Quezada C, Oz Y, Bouwman BA, et al. Robust 4C-seq data analysis to screen for regulatory DNA interactions. Nat Methods. 2012;9(10):969–72.
Article PubMed CAS Google Scholar
Landt SG, Marinov GK, Kundaje A, Kheradpour P, Pauli F, Batzoglou S, Bernstein BE, Bickel P, Brown JB, Cayting P, et al. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res. 2012;22(9):1813–31.
Article CAS PubMed PubMed Central Google Scholar
Shen L, Shao NY, Liu X, Maze I, Feng J, Nestler EJ. diffReps: detecting differential chromatin modification sites from ChIP-seq data with biological replicates. PLoS ONE. 2013;8(6):e65598.
Article CAS PubMed PubMed Central Google Scholar
Shao Z, Zhang Y, Yuan GC, Orkin SH, Waxman DJ. MAnorm: a robust model for quantitative comparison of ChIP-Seq data sets. Genome Biol. 2012;13(3):R16.
Article CAS PubMed PubMed Central Google Scholar
Raviram R, Rocha PP, Muller CL, Miraldi ER, Badri S, Fu Y, Swanzey E, Proudhon C, Snetkova V, Bonneau R, et al. 4C-ker: a method to reproducibly identify genome-wide interactions captured by 4C-Seq experiments. PLoS Comput Biol. 2016;12(3):e1004780.
Article PubMed PubMed Central CAS Google Scholar
Servant N, Varoquaux N, Lajoie BR, Viara E, Chen CJ, Vert JP, Heard E, Dekker J, Barillot E. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol. 2015;16(1):259.
Article PubMed PubMed Central CAS Google Scholar
Klein FA, Pakozdi T, Anders S, Ghavi-Helm Y, Furlong EE, Huber W. FourCSeq: analysis of 4C sequencing data. Bioinformatics. 2015;31(19):3085–91.
Article CAS PubMed PubMed Central Google Scholar
Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108–10.
Article CAS PubMed PubMed Central Google Scholar
Aitken SJ, Ibarra-Soria X, Kentepozidou E, Flicek P, Feig C, Marioni JC, Odom DT. CTCF maintains regulatory homeostasis of cancer pathways. Genome Biol. 2018;19(1):106.
Article PubMed PubMed Central CAS Google Scholar
Lau-Corona D, Suvorov A, Waxman DJ. Feminization of male mouse liver by persistent growth hormone stimulation: activation of sex-biased transcriptional networks and dynamic changes in chromatin states. Mol Cell Biol. 2017;37(19):e00301-17.
Article PubMed PubMed Central Google Scholar
Lau-Corona D, Bae WK, Hennighausen L, Waxman DJ. Sex-biased genetic programs in liver metabolism and liver fibrosis are controlled by EZH1 and EZH2. PLos Genet. 2020;16:e1008796.
Article CAS PubMed PubMed Central Google Scholar
Hark AT, Schoenherr CJ, Katz DJ, Ingram RS, Levorse JM, Tilghman SM. CTCF mediates methylation-sensitive enhancer-blocking activity at the H19/Igf2 locus. Nature. 2000;405(6785):486–9.
Article CAS PubMed Google Scholar
Renda M, Baglivo I, Burgess-Beusse B, Esposito S, Fattorusso R, Felsenfeld G, Pedone PV. Critical DNA binding interactions of the insulator protein CTCF: a small number of zinc fingers mediate strong binding, and a single finger-DNA interaction controls binding at imprinted loci. J Biol Chem. 2007;282(46):33336–45.
Article CAS PubMed Google Scholar
Reizel Y, Spiro A, Sabag O, Skversky Y, Hecht M, Keshet I, Berman BP, Cedar H. Gender-specific postnatal demethylation and establishment of epigenetic memory. Genes Dev. 2015;29(9):923–33.
Article CAS PubMed PubMed Central Google Scholar
Seitan VC, Faure AJ, Zhan Y, McCord RP, Lajoie BR, Ing-Simmons E, Lenhard B, Giorgetti L, Heard E, Fisher AG, et al. Cohesin-based chromatin interactions enable regulated gene expression within preexisting architectural compartments. Genome Res. 2013;23(12):2066–77.
Article CAS PubMed PubMed Central Google Scholar
Shen Y, Yue F, McCleary DF, Ye Z, Edsall L, Kuan S, Wagner U, Dixon J, Lee L, Lobanenkov VV, et al. A map of the cis-regulatory sequences in the mouse genome. Nature. 2012;488(7409):116–20.
Article CAS PubMed PubMed Central Google Scholar
Lin S, Lin Y, Nery JR, Urich MA, Breschi A, Davis CA, Dobin A, Zaleski C, Beer MA, Chapman WC, et al. Comparison of the transcriptional landscapes between human and mouse tissues. Proc Natl Acad Sci USA. 2014;111(48):17224–9.
Article CAS PubMed PubMed Central Google Scholar
Haller F, Bieg M, Will R, Korner C, Weichenhan D, Bott A, Ishaque N, Lutsik P, Moskalev EA, Mueller SK, et al. Enhancer hijacking activates oncogenic transcription factor NR4A3 in acinic cell carcinomas of the salivary glands. Nat Commun. 2019;10(1):368.
Article CAS PubMed PubMed Central Google Scholar
Weintraub AS, Li CH, Zamudio AV, Sigova AA, Hannett NM, Day DS, Abraham BJ, Cohen MA, Nabet B, Buckley DL. YY1 is a structural regulator of enhancer–promoter loops. Cell. 2017;171(7):1573–1588.e1528.
Article CAS PubMed PubMed Central Google Scholar
Hnisz D, Day DS, Young RA. Insulated neighborhoods: structural and functional units of mammalian gene control. Cell. 2016;167(5):1188–200.
Article CAS PubMed PubMed Central Google Scholar
Engreitz JM, Pandya-Jones A, McDonel P, Shishkin A, Sirokman K, Surka C, Kadri S, Xing J, Goren A, Lander ES, et al. The Xist lncRNA exploits three-dimensional genome architecture to spread across the X chromosome. Science. 2013;341(6147):1237973.
Article PubMed PubMed Central CAS Google Scholar
Hacisuleyman E, Goff LA, Trapnell C, Williams A, Henao-Mejia J, Sun L, McClanahan P, Hendrickson DG, Sauvageau M, Kelley DR, et al. Topological organization of multichromosomal regions by the long intergenic noncoding RNA Firre. Nat Struct Mol Biol. 2014;21(2):198–206.
Article CAS PubMed PubMed Central Google Scholar
Yin Y, Yan P, Lu J, Song G, Zhu Y, Li Z, Zhao Y, Shen B, Huang X, Zhu H, et al. Opposing roles for the lncRNA haunt and its genomic locus in regulating HOXA gene activation during embryonic stem cell differentiation. Cell Stem Cell. 2015;16(5):504–16.
Article CAS PubMed Google Scholar
Kung JT, Kesner B, An JY, Ahn JY, Cifuentes-Rojas C, Colognori D, Jeon Y, Szanto A, del Rosario BC, Pinter SF, et al. Locus-specific targeting to the X chromosome revealed by the RNA interactome of CTCF. Mol Cell. 2015;57(2):361–75.
Article CAS PubMed PubMed Central Google Scholar
Schwartzman O, Mukamel Z, Oded-Elkayam N, Olivares-Chauvet P, Lubling Y, Landan G, Izraeli S, Tanay A. UMI-4C for quantitative and targeted chromosomal contact profiling. Nat Methods. 2016;13(8):685–91.
Article CAS PubMed Google Scholar
Liang S, Kisseleva T, Brenner DA. The role of NADPH oxidases (NOXs) in liver fibrosis and the activation of myofibroblasts. Frontiers in physiology. 2016;7:17.
Article PubMed PubMed Central Google Scholar
Sun Q, Zhang W, Zhong W, Sun X, Zhou Z. Pharmacological inhibition of NOX4 ameliorates alcohol-induced liver injury in mice through improving oxidative stress and mitochondrial function. Biochimica et biophysica acta General subjects. 2017;1861(1 Pt A):2912–21.
Article CAS PubMed Google Scholar
Borrie MS, Campor JS, Joshi H, Gartenberg MR. Binding, sliding, and function of cohesin during transcriptional activation. Proc Natl Acad Sci USA. 2017;114(7):E1062–71.
Article CAS PubMed PubMed Central Google Scholar
Busslinger GA, Stocsits RR, van der Lelij P, Axelsson E, Tedeschi A, Galjart N, Peters JM. Cohesin is positioned in mammalian genomes by transcription. CTCF and Wapl. Nature. 2017;544(7651):503–7.
Article CAS PubMed PubMed Central Google Scholar
Connerney J, Lau-Corona D, Rampersaud A, Waxman DJ. Activation of male liver chromatin accessibility and STAT5-dependent gene transcription by plasma growth hormone pulses. Endocrinology. 2017;158(5):1386–405.
Article CAS PubMed PubMed Central Google Scholar
Li W, Notani D, Ma Q, Tanasa B, Nunez E, Chen AY, Merkurjev D, Zhang J, Ohgi K, Song X, et al. Functional roles of enhancer RNAs for oestrogen-dependent transcriptional activation. Nature. 2013;498(7455):516–20.
Article CAS PubMed PubMed Central Google Scholar
Sigova AA, Abraham BJ, Ji X, Molinie B, Hannett NM, Guo YE, Jangi M, Giallourakis CC, Sharp PA, Young RA. Transcription factor trapping by RNA in gene regulatory elements. Science. 2015;350(6263):978–81.
Article CAS PubMed PubMed Central Google Scholar
Fang B, Everett LJ, Jager J, Briggs E, Armour SM, Feng D, Roy A, Gerhart-Hines Z, Sun Z, Lazar MA. Circadian enhancers coordinate multiple phases of rhythmic gene transcription in vivo. Cell. 2014;159(5):1140–52.
Article CAS PubMed PubMed Central Google Scholar
Schwalie PC, Ward MC, Cain CE, Faure AJ, Gilad Y, Odom DT, Flicek P. Co-binding by YY1 identifies the transcriptionally active, highly conserved set of CTCF-bound regions in primate genomes. Genome Biol. 2013;14(12):R148.
Article PubMed PubMed Central CAS Google Scholar
Keleher MR, Zaidi R, Hicks L, Shah S, Xing X, Li D, Wang T, Cheverud JM. A high-fat diet alters genome-wide DNA methylation and gene expression in SM/J mice. BMC Genomics. 2018;19(1):888.
Article CAS PubMed PubMed Central Google Scholar
Lodato NJ, Melia T, Rampersaud A, Waxman DJ. Sex-differential responses of tumor promotion-associated genes and dysregulation of novel long noncoding RNAs in constitutive androstane receptor-activated mouse liver. Toxicol Sci. 2017;159(1):25–41.
Article CAS PubMed PubMed Central Google Scholar
Lodato NJ, Rampersaud A, Waxman DJ. Impact of CAR agonist ligand TCPOBOP on mouse liver chromatin accessibility. Toxicol Sci. 2018;164(1):115–28.
Article CAS PubMed PubMed Central Google Scholar
Leung A, Trac C, Du J, Natarajan R, Schones DE. Persistent chromatin modifications induced by high fat diet. J Biol Chem. 2016;291(20):10446–55.
Article CAS PubMed PubMed Central Google Scholar
Marques AC, Hughes J, Graham B, Kowalczyk MS, Higgs DR, Ponting CP. Chromatin signatures at transcriptional start sites separate two equally populated yet distinct classes of intergenic long noncoding RNAs. Genome Biol. 2013;14(11):R131.
Article PubMed PubMed Central Google Scholar
Gil N, Ulitsky I. Production of spliced long noncoding RNAs specifies regions with increased enhancer activity. Cell Syst. 2018;7(5):537–547.e533.
Article CAS PubMed PubMed Central Google Scholar
Lieberman-Aiden E, van Berkum NL, Williams L, Imakaev M, Ragoczy T, Telling A, Amit I, Lajoie BR, Sabo PJ, Dorschner MO, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science. 2009;326(5950):289–93.
Article CAS PubMed PubMed Central Google Scholar
Conforto TL, Waxman DJ. Sex-specific mouse liver gene expression: genome-wide analysis of developmental changes from pre-pubertal period to young adulthood. Biol Sex Diff. 2012;3:9.
Article CAS Google Scholar
Lowe R, Gemma C, Rakyan VK, Holland ML. Sexually dimorphic gene expression emerges with embryonic genome activation and is dynamic throughout development. BMC Genomics. 2015;16:295.
Article PubMed PubMed Central CAS Google Scholar
Wauthier V, Sugathan A, Meyer RD, Dombkowski AA, Waxman DJ. Intrinsic sex differences in the early growth hormone responsiveness of sex-specific genes in mouse liver. Mol Endocrinol. 2010;24(3):667–78.
Article CAS PubMed PubMed Central Google Scholar
Zhang Y, Laz EV, Waxman DJ. Dynamic, sex-differential STAT5 and BCL6 binding to sex-biased, growth hormone-regulated genes in adult mouse liver. Mol Cell Biol. 2012;32(4):880–96.
Article CAS PubMed PubMed Central Google Scholar
Conforto TL, Steinhardt GFT, Waxman DJ. Cross talk between gh-regulated transcription factors HNF6 and CUX2 in adult mouse liver. Mol Endocrinol. 2015;29(9):1286–302.
Article CAS PubMed PubMed Central Google Scholar
Wang H, Maurano MT, Qu H, Varley KE, Gertz J, Pauli F, Lee K, Canfield T, Weaver M, Sandstrom R, et al. Widespread plasticity in CTCF occupancy linked to DNA methylation. Genome Res. 2012;22(9):1680–8.
Article CAS PubMed PubMed Central Google Scholar
Takasugi M, Hayakawa K, Arai D, Shiota K. Age- and sex-dependent DNA hypomethylation controlled by growth hormone in mouse liver. Mech Ageing Dev. 2013;134(7–8):331–7.
Article CAS PubMed Google Scholar
Xu Y, Guo W, Li P, Zhang Y, Zhao M, Fan Z, Zhao Z, Yan J. Long-range chromosome interactions mediated by cohesin shape circadian gene expression. PLoS Genet. 2016;12(5):e1005992.
Article PubMed PubMed Central CAS Google Scholar
Yeung J, Mermet J, Jouffe C, Marquis J, Charpagne A, Gachon F, Naef F. Transcription factor activity rhythms and tissue-specific chromatin interactions explain circadian gene expression across organs. Genome Res. 2018;28(2):182–91.
Article CAS PubMed PubMed Central Google Scholar
Mermet J, Yeung J, Hurni C, Mauvoisin D, Gustafson K, Jouffe C, Nicolas D, Emmenegger Y, Gobet C, Franken P, et al. Clock-dependent chromatin topology modulates circadian transcription and behavior. Genes Dev. 2018;32(5–6):347–58.
Article CAS PubMed PubMed Central Google Scholar
Phanstiel DH, Van Bortle K, Spacek D, Hess GT, Shamim MS, Machol I, Love MI, Aiden EL, Bassik MC, Snyder MP. Static and dynamic DNA loops form AP-1-bound activation hubs during macrophage development. Mol Cell. 2017;67(6):1037–1048.e1036.
Article CAS PubMed PubMed Central Google Scholar
Banerjee S, Das RK, Shapiro BH. Feminization imprinted by developmental growth hormone. Mol Cell Endocrinol. 2019;479:27–38.
Article CAS PubMed Google Scholar
Lowe R, Gemma C, Rakyan VK, Holland ML. Sexually dimorphic gene expression emerges with embryonic genome activation and is dynamic throughout development. BMC Genomics. 2015;16(1):295.
Article PubMed PubMed Central CAS Google Scholar

Download references

Acknowledgements

We thank Dr. Andy Rampersaud of the Waxman Laboratory for developing the ChIP-seq analysis pipeline used in this study.

Funding

Supported in part by NIH grants DK121998 and DK33765 from the NIH (to DJW) and a National Sciene Foundation predoctoral fellowship (DGE-1247312, to BJM). The funding sponsor played no role in the design of the study and collection, analysis, and interpretation of data or in writing the manuscript.

Author information

Authors and Affiliations

Department of Biology and Bioinformatics Program, Boston University, 5 Cummington Mall, Boston, MA, 02215, USA
Bryan J. Matthews & David J. Waxman

Authors

Bryan J. Matthews
View author publications
You can also search for this author in PubMed Google Scholar
David J. Waxman
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

BJM and DJW conceived and designed the study. BJM carried out all of the laboratory experiments and data analyses. BJM and DJW jointly wrote and edited the manuscript for publication. Both authors read and approved the final manuscript.

Corresponding author

Correspondence to David J. Waxman.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Figure S1

. Consistency of biological replicates and additional features of sex-biased differential CTCF/cohesin peaks. A-D. Consistency of biological replicate ChIP-seq samples, with merged tracks and annotations as in Figs. 2 and 3. Here, we additionally show all cohesin (n = 3 per sex) and CTCF ChIP-seq replicates (n = 4 per sex). Blue, samples from male mouse liver; red, samples from female mouse liver are shown in red. A. Biological replicates for a male-biased “Lone CTCF” site nearby Mir2136. B. Biological replicates for two female-biased CAC(ΔCoh) peaks nearby 5530601H04Rik. Higher cohesin signal is observed across replicates for both the weaker upstream peak and a strong downstream peak. CTCF signal is not significantly different at either peak. C. Biological replicates for a female-biased CAC peak (ΔCoh AND ΔCTCF) upstream of strong sex-independent peak in an intron of Moh9. D. Biological replicates for several sex-independent CAC peaks of various peak strengths. Replicates are consistent at all sites regardless of peak strength. E. The limited overlap of sex-biased CTCF and sex-biased cohesin binding, seen in Fig. 1B, is not an artifact of the thresholds used for peak filtering. In the data shown here, similar to Fig. 1B, the overlap of sex-biased CTCF and sex-biased cohesin binding sites is limited. Shown are Venn diagrams with the number of overlapping sex-biased cohesin (Left; blue) and CTCF sites (Right; purple) for the following groups (top to bottom): (1) All sex-differential sites output by diffReps without filtering for overlap of a MACS2 ChIP-seq peak; (2) All sex-differential diffReps sites (rather than peaks, which may contain multiple sites) from diffReps that overlap a MACS2 peak for a given factor; (3) All sex-differential hotspots identified by diffReps, which is an alternate method in the diffReps software package to identify differential sites. Specifically, this last approach looks for clusters of differentially co-regulated sites that might be missed by simple overlap analysis. Overlap for all Venn diagrams is defined as ≥ 1 bp overlapping using bedtools. In some cases, two Rad21 peaks overlap a CTCF peak, or vice versa, and therefore, the number of overlapping cohesin (Rad21) sites does not necessarily equal the number of CTCF sites (hence, two numbers in the Venn overlap). F. Pie charts showing the class distribution of each sex-biased CTCF/cohesin peak set (from top to bottom): male-biased ΔCohesin peaks, female-biased ΔCohesin peaks, male-biased ΔCTCF peaks, and female-biased ΔCTCF peaks. For each of these four groups, the fraction of peaks at CAC sites is shown in purple while the fraction of peaks at either CNC (for ΔCohesin) or Lone CTCF (for ΔCTCF) is shown in blue. The total number of differential peaks in each group is indicated below each chart. Overall, female-biased sites comprised a higher percentage of CAC sites than male-biased sites. Consequently, a larger percentage of male-biased peaks are CNC peaks (for ΔCohesin peaks) and Lone CTCF peak (for ΔCTCF peaks). Peak numbers here differ slightly from Fig. 1B for cohesin differential peaks, but not CTCF, because of our approach to categorizing peaks as CNC or CAC for cohesin peaks (see “Materials and methods”). For CTCF we defined CAC peaks as genomic regions bound by CTCF that were also bound by cohesin in a majority of individual cohesin replicates (2 or 3 out of a total n = 3 per sex). Using the same approach for cohesin, we defined CAC peaks as genomic regions bound by cohesin that were also bound by CTCF in a majority of individual CTCF replicates (3 or 4 out of a total n = 4 per sex). If a peak was bound by none or only a minority of replicates for the opposite factor then it was considered Lone CTCF (in the case of CTCF; 0 or 1 cohesin replicates overlapping) or CNC (in the case of cohesin; 0 or 1 CTCF replicates overlapping). As CTCF has n = 4 replicates, if a cohesin peak is bound by exactly 2 individual CTCF replicates (of the same sex) then it is not classified and is excluded from downstream analyses. 54 male-biased cohesin peaks overlap 2 male CTCF replicates and 36 female-biased cohesin peaks overlap 2 female CTCF replicates (value of 2 in column H of Additional file 3: Table S1B). All overlaps were performed using bedtools with a minimum overlap of 1 bp, and all comparisons were made separately for males and females. Figure S2. Comparison of sex-biased CTCF/cohesin peaks. A. Female-biased CTCF and cohesin peaks tend to be stronger than male-biased peaks. Shown here are box plots for ChIP-seq signal for CTCF and cohesin for both ΔCohesin and ΔCTCF peaks. These plots differ from those presented in Fig. 2A, which present normalized ChIP signal for the factor with differential signal (i.e., male and female cohesin ChIP-seq signal for ΔCohesin peaks). In aggregate, CAC peaks with significant sex-biased cohesin binding show the same directionality of sex-bias for CTCF (and vice versa), albeit at a reduced magnitude (see also Fig. 1C). The y-axis shows normalized ChIP-seq signal for the groups indicated along the x-axis. Peaks with male-biased and female-biased cohesin binding (Left) and CTCF binding (Right) are presented separately. Each group of 4 box plots represents the male and female ChIP-seq signal for cohesin, followed by the corresponding ChIP-seq signals for CTCF for the same set of peaks. Each plot represents all differential peaks for a given sex (male or female) and factor (CTCF or cohesin). These four datasets are further divided by peak type (CAC or CNC for ΔCohesin peaks, and CAC or Lone CTCF for ΔCTCF peaks), as indicated below the x-axis. Peak scores are calculated by average intra-peak ChIP signal, normalized by the total sequence reads per million in peak (RIPM; see “Materials and methods”). B. Female-biased CAC peaks contain higher quality CTCF motifs than male-biased CAC peaks [p = 0.0212 for CAC(ΔCoh) and p = 0.0023 for CAC(ΔCTCF) peaks; M-W t test], as reflected by the FIMO motif score. This log-likelihood ratio score is a reflection of how close the best intra-peak motif matches the canonical core CTCF motif MA0139.1. There is no significant difference between motif scores for male-biased and female-biased Lone CTCF, or for male-biased and female-biased CNC peaks (p = 0.7671 and p = 0.1329; M-W t-test). The dashed line at FIMO score = 10 reflects the cutoff used to define the presence or absence of a motif in Additional file 1: Figure S2C. C. CTCF Motif frequency, based on presence of CTCF motif (MA0139.1) as identified by FIMO, with a minimum score of 10. The y-axis shows the percent of peaks in each group (separately for male-biased, female-biased, and sex-independent) found to have a CTCF motif within the peak region. A larger fraction of female-biased than male-biased CAC peaks was found to contain a CTCF binding motif. In contrast, a larger fraction of male-biased Lone CTCF peaks contain a CTCF motif, despite no significant difference in peak strength between male-biased and female-biased Lone CTCF peaks. A larger fraction of female-biased CNC peaks contain a CTCF motif, however, the vast majority do not contain CTCF motifs, as expected (< 20% for all groups). In all cases, the percent for each group is comparable to a matched set of sex-independent peaks. D. Proportion of male-biased and female-biased CTCF peaks that have: no CTCF motif (gray), a CTCF motif lacking a CpG dinucleotide (orange), or a CTCF motif containing a CpG dinucleotide (orange). “All” indicates any male-biased or female-biased CTCF peak. E. Female-biased intra-TAD on Chr19 that contains 12 sex-biased genes, shown in blue and red boxes, some of which are lncRNA genes (ncRNA gene designations, in green). Inset at bottom left of the Figure shows CTCF and cohesin (Rad21) ChIP-seq tracks for male and female mouse liver surrounding the 3′ boundary of the female-biased intra-TAD. F. Intra-TAD loops and loop anchors are mostly shared between male and female mouse liver. Using a computational intra-TAD loop prediction algorithm [31], we used the cohesin and CTCF ChIP-seq datasets for male and female mouse liver to identify 9543 intra-TAD loops in male liver [31] and 9724 loops in female liver, respectively. 87.9% of the intra-TAD loops in male liver were also identified in female liver (left), and 93.4% of the male intra-TAD loop anchors are also predicted to be loop anchors in female liver. This finding is consistent with there being a limited number of autosomal CAC peaks with sex differences in CTCF and cohesin binding (53 total) (Additional file 3: Table S1C) To account for nested loop structures, shared loops were defined as loops with a reciprocal overlap of 80% or greater between the loops, as implemented in prior studies of CAC-mediated insulating loops [26, 31]. Figure S3. Tissue conservation of liver sex-differential CTCF and cohesin peaks in ENCODE mouse consortium datasets. A. The x-axis indicates the number of male mouse tissues other than liver where CTCF is bound, out of 15 tissues examined by the ENCODE Consortium. A value of 15 indicates tissue-ubiquitous CTCF binding, and a value of 0 indicates liver-specific CTCF binding. The y-axis shows the proportion of male-biased peaks (blue) or female-biased peaks (red) that fall into a given bin. P-values comparing the distribution of tissue-specificity values for CTCF binding between males and females are indicated in the upper left corner of each plot (KS t-test). Results show that CAC sites (upper panels) are much less tissue-specific than Lone CTCF and CNC sites (lower panels). Further, female-biased CAC peaks are less tissue-specific than male-biased CAC peaks, while male-biased Lone CTCF peaks are less liver-specific than female-biased peaks of the same class. The greater tissue ubiquity of CTCF binding for female-biased CAC peaks could be due to the fact that female-biased CTCF peaks are stronger and contain higher quality CTCF motifs, insofar as stronger peaks show greater conservation for both sex-differential and all CTCF peaks (see panel B, below); however, male-biased Lone CTCF peaks are not significantly stronger, nor do they contain higher quality motifs than the female-biased Lone CTCF peaks. The apparent difference could be due to the fact that the CTCF ChIP-seq data from the non-liver, non-reproductive tissues examined here was obtained by the ENCODE consortium from male mice [60]. Very few male-biased and female-biased CNC peaks were bound by CTCF in other any other mouse tissues (< 20% of the total sex-biased liver CNC sites). This finding provides additional evidence that CNC sites are found at liver-specific cis regulatory elements, and that these sites do not act as insulators in other non-liver tissues (i.e., CTCF binding is lost in liver or gained in some other tissue). B. There is a strong association between CTCF peak strength and tissue conservation of CTCF binding, which likely explains the modestly higher tissue conservation of female-biased CTCF and cohesin peaks seen in panel A. Shown on the y-axis are reads-in-peaks normalized ChIP-seq data for all CTCF peaks, male-biased CTCF peaks, and female-biased CTCF peaks. These are grouped according to the number of non-liver ENCODE tissues with a CTCF peak, where 0 indicates a peak is liver-specific and 15 means all male mouse tissues with ENCODE datasets have CTCF bound at that position, as in panel A. C. The tissue specificity of neighboring genes varies significantly with the class of CTCF/cohesin binding site. Shown are Tau values, where a value close to 1 indicates the pattern of expression across mouse ENCODE RNA-seq datasets is highly tissue-specific, and where Tau values less than ~ 0.3 indicate housekeeping genes. Nearest genes (within 20 kb) were defined based on distance to the TSS, and only liver-expressed genes were considered (FPKM > 1). Tau values were calculated based on the average of two replicates from all tissues except testis, using expression data generated by the ENCODE consortium. Both female-biased and male-biased CNC sites are near genes that generally are more tissue-specific than liver-expressed genes. In addition, genes near female-biased CNC sites are significantly different than genes near similarly sex-biased CACs (p = 0.007; M-W test). This difference is not a reflection of the male-biased or female-biased CAC group used in the comparison, insofar as genes near female-biased CNCs are significantly more liver-specific than genes near male-biased CAC sites (p = 0.0171; M-W test), while the opposite comparison for male-biased CNCs vs female-biased CACs is still not significant (p = 0.07; M-W test). For these analyses, liver-expressed genes are defined by a liver expression value of FPKM > 1 (8810 genes in total) and mapping was based on the closest TSS within 20 kb of a peak. Figure S4. Screenshot of TAD containing C8a/C8b, and cohesin depletion effects. A. Shown is a screenshot with proposed model linking the distal male-biased enhancers and component complement genes C8a and C8b within a single TAD on mouse chromosome 4. This screenshot spans chr4:102960671-104603975. The tracks, normalization, and annotations are as described in Fig. 3. B. For the proximally regulated male-biased genes shown in Fig. 3B (Nat8 and Cml5), depletion of cohesin does not significantly impact gene expression. Figure S5. Quality control of 4C-seq library. A. Agarose gel analysis for quality control of ligated, digested, and re-ligated 4C samples. Lane (i) analyzes a sample after proximity ligation, lane (ii) shows the sample after digestion with the restriction enzyme Csp6i, and lane (iii) shows the sample after self-circularization ligation. Lane (iii) represents the final material used as input for inverse PCR with viewpoint-specific primers (Additional file 2: Table S3A). DNA fragment sizes (in kb) are marked on the left of the gel. B. Agarose gel analysis for quality control of final 4C-seq libraries after the inverse PCR step. A diverse mixture of PCR products is present, as indicated by a smear on the gel, with sizes primarily below ~ 1 kb, which indicate a high-quality library and which allows for efficient sequencing. DNA fragment sizes (in kb) are marked on the right of the gel. C. Shown is RNA-seq expression data for A1bg and 12 mono-exonic lncRNAs (see Fig. 4A), obtained in six separate RNA-seq datasets from CD-1 mouse liver, and one RNA-seq dataset from C57/Bl6 mouse liver. The first two columns indicate the maximum expression (in FPKM) for male and female liver across these datasets. Following this from left to right, the columns indicate the mean expression level of each gene in female liver (FPKM values) for: Total PolyA + unstranded RNA-seq [sequencing series G83; [54]], Total PolyA + unstranded RNA-seq [sequencing series G85; [74]], Total PolyA + stranded RNA-seq [sequencing series G118; [54]], Nuclear PolyA + stranded RNA-seq [sequencing series G119; [54]], Total Ribosomal RNA-depleted stranded RNA-seq [sequencing series G118; [54]], Nuclear Ribosomal RNA-depleted stranded RNA-seq [sequencing series G119; [54]], and Total Ribosomal RNA-depleted stranded RNA-seq [from C57Bl6/J, all others CD-1; [98]]. The final columns indicate the nuclear enrichment for PolyA + RNA datasets and for Ribosomal RNA-depleted RNA-seq datasets (linear scale). Specifically, for PolyA + datasets, this is the FPKM value in data column 6 (“Nuclear Poly A Strnd”) divided by the FPKM value in data column 5 (“Total Poly A Strnd”). Similarly, the final column is calculated for Ribosomal RNA-Depleted Nuclear versus total RNA-seq expression (FPKM in data column 8 “Nuclear RiboM Strnd” divided by the FPKM value in data column 7 “Total RiboM Strnd”). D. Shown are the log2 male/female expression ratios for A1bg and the 12 mono-exonic lncRNAs (Fig. 4A) for the seven RNA-seq datasets described in panel C. Fold change is calculated by EdgeR for all datasets and the order of datasets is the same as in panel C. The final column indicates the number of datasets in which the sex differences indicated are significant (FDR < 0 0.05; EdgeR). Figure S6. 4C-seq data tracks for individual male and female mouse livers. Shown are 4C-seq data for the same genomic regions presented in Fig. 4, but showing the individual 4C-seq data for each of 3 male liver and 3 female liver biological replicates per viewpoint. The gene tracks and sex-biased sites are as described in Fig. 4. These bed file tracks are as follows (from top to bottom): sex-biased H3K27ac peaks, sex-biased DHS, sex-biased cohesin peaks, and sex-biased CTCF peaks. Protein coding genes, sex-biased lncRNA genes, and intra-TAD loops are also shown, where present. A. All six biological replicates for the A1bg enhancer viewpoint (chr15:60733512-60954051). B. All six biological replicates for the Gm4794 enhancer viewpoint (chr10:33418446-33680888). C. All six biological replicates for the C9 promoter viewpoint (chr15:6147917-6461799). D. All six biological replicates for the Nudt7 enhancer viewpoint (chr8:116592444-116707613). E. Distal enhancer regions from panel C with bidirectional eRNA loci and YY1 binding sites indicated in relation to male-biased enhancers (chr15:6,164,877–6,184,670). Figure S7. High expression of Nox4 in hepatocellular carcinoma and 4C-seq biological replicates. A. Nox4 is highly upregulated in tumor relative to normal healthy tissue of mice that spontaneously develop tumors (mouse strain C3H; [53]) (p = 0.0006, M-W t-test). B. In human patient samples, Nox4 is consistently upregulated in tumor tissue relative to normal controls, including in hepatocellular carcinoma (marked by thick black line). Only data for cancer types with matched primary tumor and normal tissue controls are shown, with Nox4 showing significant up regulation in tumors for 14 of 18 matched tissue pairs. Expression datasets are from The Cancer Genome Atlas (TCGA) and were analyzed by the Tumor IMmune Estimation Resource (TIMER) webtool (https://cistrome.shinyapps.io/timer/) with default parameters [52]. The significance of comparisons between Normal and Tumor tissue was calculated by Wilcoxon test and is indicated on the chart as: 0 ≤ *** < 0.001 ≤ ** < 0.01 ≤ * < 0.05. For example, the difference between liver tumor versus normal tissue (p = 1.35E − 25) is indicated as ***. C. Chromatin state key for top browser track in Fig. 5A, which is based on chromatin states in both male and female mouse liver, which were determined for the entire genome based on a 14-state model of chromatin states developed in [14]. D. Male and female 4C-seq biological replicates for the Nox4 enhancer-1 viewpoint (chr7:94248242-94726358). Tracks are as described in Additional file 1: Figure S6. “SB Peaks” indicate the male-biased (blue) and female-biased (pink) for the indicated factor or assay and a darker color indicates a more stringently defined sex-biased region (cutoffs differ for each factor; see “Materials and methods”). These tracks, from top to bottom, are as follows: H3K27ac ChIP-seq, DNase-seq, Rad21 ChIP-seq, and CTCF ChIP-seq and are as described in Fig. 3 and Fig. 4. A complete listing of all such peaks is provided in Additional file 5: Tables S2D and E (K27ac ChIP-seq and DNase-seq) and in Additional file 3: Tables S1D, and 1E (CTCF and cohesin ChIP-seq peaks). E. Male and female 4C-seq biological replicates for the Nox4 enhancer-2 viewpoint (chr7:94248242-94726358). Tracks are as described in Additional file 1: Figure S6 and in panel D of this Figure. F. Three separate models of regulatory domain prediction support the model that Nox4 enhancer-1 viewpoint and enhancer-2 viewpoint (VP1, VP2) do not interact. The position of enhancer-1 viewpoint is indicated by a red arrow, and the position of enhancer-2 viewpoint is indicated by a green arrow. From top to bottom: Enhancer Promoter Units (EPUs) are domains based on the pan-tissue correlation of ChIP-seq signal at enhancer and promoters [60]; intra-TAD loops are computationally predicted based on CTCF and cohesin ChIP-seq data for male mouse liver [31]; cohesin ChIA-PET loops show cohesin-anchored interactions from mouse embryonic stem cells [26].

Additional file 2: Table S3.

A. 4C-seq Primer Sequences. B. ChIP-seq Stats. C. 4C-seq Stats.

Additional file 3: Table S1.

A. CTCF Sex Differential Peaks. B. Cohesin Sex Differential Peaks. C. CAC Sex Differential Peaks. D. CTCF Trackline. E. Cohesin Trackline. F. Male CTCF Peaks. G. Female CTCF Peaks. H. Male Coh peaks. I. Female Coh peaks. J. M/F iTADs. K. F iTAD Anchors.

Additional file 4: Table S4.

4CSeq.

Additional file 5: Table S2.

A. PCG Coh Depl Expr. B. lncRNA Coh Depl Expr. C. M-biased Genes CohDepl. D. H3K27ac trackline. E. DHS trackline.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Cite this article

Matthews, B.J., Waxman, D.J. Impact of 3D genome organization, guided by cohesin and CTCF looping, on sex-biased chromatin interactions and gene expression in mouse liver. Epigenetics & Chromatin 13, 30 (2020). https://doi.org/10.1186/s13072-020-00350-y

Download citation

Received: 19 February 2020
Accepted: 03 July 2020
Published: 17 July 2020
DOI: https://doi.org/10.1186/s13072-020-00350-y

Impact of 3D genome organization, guided by cohesin and CTCF looping, on sex-biased chromatin interactions and gene expression in mouse liver

Abstract

Introduction

Materials and methods

Mouse protocols, extraction of liver nuclei, and chromatin preparation

Chromatin immunoprecipitation and sequencing (ChIP-seq)

4C-seq methods

Computational analysis of ChIP-seq datasets

Discovery of intra-TAD loops

Computational analysis of 4C-seq datasets

Impact of cohesin depletion sex-biased gene expression

Analysis of other datasets

Statistical analysis

Data availability

Results

Sex differences in CTCF and cohesin binding in mouse liver

Sex-biased binding sites mapping to TAD and intra-TAD boundaries, genes and regulatory elements

Impact of cohesin depletion on distally regulated male-biased genes

4C-seq analysis of sex-biased chromatin interactions in mouse liver

A1bg region

Gm4794 and Sult3a1 region

C9 region

Sex-independent, nested intra-TAD loops restrict Nox4 to proximal enhancer–promoter interactions

Discussion

Availability of data and materials

References

Acknowledgements

Funding

Author information

Authors and Affiliations

Contributions

Corresponding author

Ethics declarations

Ethics approval and consent to participate

Consent for publication

Competing interests

Additional information

Publisher's Note

Supplementary information

Additional file 1: Figure S1

Additional file 2: Table S3.

Additional file 3: Table S1.

Additional file 4: Table S4.

Additional file 5: Table S2.

Rights and permissions

About this article

Cite this article

Share this article

Keywords

Epigenetics & Chromatin

Contact us