|Home | About | Journals | Submit | Contact Us | Français|
Chromatin state maps were developed to elucidate sex differences in chromatin structure and their impact on sex-differential chromatin accessibility and sex-biased gene expression in mouse liver. Genes in active, inactive, and poised chromatin states exhibited differential responsiveness to ligand-activated nuclear receptors and distinct enrichments for functional gene categories. Sex-biased genes were clustered by chromatin environments and mapped to DNase-hypersensitive sites (DHS) classified by sex bias in chromatin accessibility and enhancer modifications. Results were integrated with genome-wide binding data for five transcription factors implicated in growth hormone-regulated, sex-biased liver gene expression, leading to the following findings. (i) Sex-biased DHS, but not sex-biased genes, are frequently characterized by sex-differential chromatin states, indicating distal regulation. (ii) Trimethylation of histone H3 at K27 (H3K27me3) is a major sex-biased repressive mark at highly female-biased but not at highly male-biased genes. (iii) FOXA factors are associated with sex-dependent chromatin opening at male-biased but not female-biased regulatory sites. (iv) Sex-biased STAT5 binding is enriched at sex-biased DHS marked as active enhancers and preferentially targets sex-biased genes with sex-differences in local chromatin marks. (v) The male-biased repressor BCL6 preferentially targets female-biased genes and regulatory sites in a sex-independent chromatin state. (vi) CUX2, a female-specific repressor of male-biased genes, also activates strongly female-biased genes, in association with loss of H3K27me3 marks. Chromatin states are thus a major determinant of sex-biased chromatin accessibility and gene expression, with FOXA pioneer factors proposed to confer sex-dependent chromatin opening and STAT5, but not BCL6, regulating sex-biased genes by binding to sites in a sex-biased chromatin state.
The epigenetic environment is controlled by a variety of factors, most notably chromatin modifications, which regulate DNA methylation, chromatin accessibility, and transcription factor (TF) binding. Genomic regions active in a given cell type or regulatory state can be identified by the distinct and characteristic patterns of histone methylations, acetylations, and other chromatin modifications associated with transcribed genes, transcription start sites (TSS), and enhancers. Expressed genes can thus be distinguished from silenced genes and from genes that are poised for expression based on their unique epigenetic environment (1, 2). Enhancers exhibit cell-type-specific patterns of local chromatin modifications, which confer cell-type-specific activity and enable cell-type-specific regulation of distal target genes (3). These characteristic patterns of chromatin modifications can, in turn, be used to predict gene expression in multiple cell types (4).
Hypersensitivity to DNase I cleavage is commonly used to assay chromatin accessibility, which is a key feature of cell-type-specific regulation (5, 6). Regulation of chromatin opening is an important regulatory event in the binding of cell-type-specific TFs, which are invariably (>90%) bound at DNase-hypersensitive sites (DHS) (5). Factors regulating cell-type-specific chromatin opening and TF binding include differences in nucleosome occupancy and positioning, chromatin modifications, and other features of the chromatin environment (7–10). Thus, cell-type-specific gene expression involves a complex interplay between TF activators and repressors, chromatin modifiers, and the chromatin environment. Some TFs can bind to closed chromatin and facilitate nucleosome repositioning and/or chromatin opening (e.g., FOXA “pioneer factors”) (11, 12), while other factors preferentially bind open chromatin marked by cell-type- or condition-specific modifications (13–15). Comprehensive approaches, involving the integration of genome-wide chromatin state maps with binding site maps for TFs and other regulatory factors, are therefore needed to unravel the complexity of mammalian transcriptional networks (13, 14, 16).
Sex differences in gene expression are common and result from differences in genotypic sex as well as sex differences in circulating and local hormonal regulators (17). Sex differences in gene expression are best characterized in the liver, where ~1,000 genes are expressed in a sex-biased manner, as seen in mice, rats, and humans (18–20). These sex differences affect diverse physiological processes, including hepatic steroid, lipid, and drug metabolism (21–23), and contribute to sex differences in cardiovascular disease risk, fatty liver disease, and the development of hepatocellular carcinoma (20, 24–26). Sex differences in liver gene expression are regulated by the temporal pattern of pituitary growth hormone (GH) secretion, which is sex dependent and programmed by androgen and estrogen exposure during the neonatal period (21). In males, GH is released from the pituitary gland in pulses at regular intervals, with little or no plasma GH detectable between GH secretory events, whereas the more frequent pituitary GH release in females gives a more continuous pattern of circulating hormone. Consequently, GH signaling to downstream transcriptional mediators is pulsatile in male liver and is persistent in female liver (27, 28).
Several TFs contribute to the sex-dependent actions of GH in the liver. STAT5, a major effector of the transcriptional actions of GH, is activated by GH-induced tyrosine phosphorylation, which leads to STAT5 dimerization and translocation into the nucleus (29). Liver nuclear STAT5 DNA-binding activity is pulsatile in males but persistent in females, mirroring the temporal plasma GH profiles of each sex (27, 28, 30). STAT5 is essential for sex-biased gene expression in mouse liver, where 75 to 82% of sex-differential expression is lost when Stat5b is inactivated (31). BCL6, a male-biased, GH-regulated transcriptional repressor, competes with STAT5 binding and antagonizes GH-stimulated STAT5 transcriptional activity (30, 32). CUX2 is a female-specific, continuous GH-inducible repressor (33) that downregulates many male-biased genes and activates a subset of female-biased genes in female liver (34). Other liver-expressed TFs, including HNF4A and HNF6, are also important for hepatic sex differences, as indicated by knockout studies and reporter gene analyses and by the enrichment of their motifs at GH-dependent regulatory sites (35–39).
Sex differences in gene expression have been linked to sex differences in GH-regulated chromatin accessibility, as revealed by genome-wide DHS mapping (39). Sex-biased STAT5 binding in mouse liver correlates with sex differences in DHS accessibility and the presence of activating chromatin marks and inversely correlates with the repressive mark trimethylation of histone H3 at K27 (H3K27me3) (30). GH-induced STAT5 DNA binding requires chromatin opening and activating histone modifications, as seen at several early GH response genes in rat liver (40). However, little is known about the chromatin states of sex-biased genes and their DHS regulatory sites and how the overall epigenetic environment might impact binding and regulation by STAT5 and other TFs required for sex-biased liver gene expression.
Here, we characterize chromatin environments on a global scale in both male and female mouse liver based on genome-wide data for four activating chromatin marks, two repressive chromatin marks, and global DHS maps. We identify genes in active, inactive, and poised chromatin states, and we establish the propensity of genes in each state for induction or repression by ligand-activated nuclear receptors, which mediate metabolic responses to myriad steroids, drugs, and environmental chemicals (41–43). Further, we elucidate sex-differential chromatin states at sex-biased genes and sex-biased DHS, and we integrate these chromatin state data with binding site data from chromatin immunoprecipitation with high-throughput sequencing (ChIP-seq) for the GH-regulated TFs STAT5 (30), BCL6 (30), and CUX2 (34) and for the pioneer factors FOXA1 and FOXA2 (11, 12), which impart male bias to hepatocellular carcinoma by facilitating androgen-mediated tumor promotion in males and estrogen-dependent resistance to tumorigenesis in females (44). Our findings reveal the importance of sex-dependent chromatin states in the control of liver gene expression. Strikingly, we find evidence for distinct sets of regulatory mechanisms in each sex. Many female-biased genes are associated with sex-independent chromatin marks and are preferentially repressed by the male-biased repressor BCL6, while a smaller subset, comprised of highly female-biased genes, shows sex differences in proximal chromatin marks and preferential activation by female-enriched STAT5 binding. Furthermore, we identify H3K27me3-based repression as an important mechanism regulating highly female-biased but not male-biased genes, and we identify sex differences in H3K4me1 (methylation of histone H3 at K4) profiles at male-biased but not female-biased DHS that suggest sex-dependent nucleosome repositioning facilitated by FOXA factors, in particular FOXA1. These findings are integrated into a model for GH-regulated transcriptional and epigenetic control of sex-biased gene expression.
Chromatin isolation from formaldehyde cross-linked liver nuclei prepared from individual 7- to 8-week-old male and female CD-1 mouse livers (Charles River Laboratories, Kingston, NY) obtained using protocols approved by the Boston University Institutional Animal Care and Use Committee, followed by sonication and chromatin immunoprecipitation (ChIP), was carried out by E. Laz of this laboratory as reported previously (30). The following ChIP-validated antibodies were used: H3K27me3 (ab6002; Abcam) (2 μg per 10 μg of sonicated DNA), H3K9me3 (ab8898; Abcam) (0.8 μg per 10 μg of sonicated DNA), H3K4me3 (ab8580; Abcam) (0.8 μg per 10 μg of sonicated DNA), H3K4me1 (ab8895; Abcam) (0.8 μg per 10 μg of sonicated DNA), histone H3 acetylated at K27 (H3K27ac) (ab4729), H3K36me3 (ab9050), and normal rabbit IgG (sc-2027; Santa Cruz) (7.5 to 8 μg per 100 μg of sonicated DNA). ChIP DNA was analyzed by quantitative PCR (qPCR) (30) to interrogate genomic regions selected as positive controls and negative controls for each antibody based on initial ChIP-seq results.
Liver genomic DNA isolated by ChIP was prepared for sequencing using a SPRI-TE Nucleic Acid Extractor for size selection (Beckman Coulter Genomics, Danvers, MA) followed by PCR enrichment with bar coding (30). Sample preparation and 35- to 41-nucleotide (nt) single-end read sequencing on an Illumina GAII or Illumina HiSeq2000 instrument, followed by mapping to the mouse genome build mm9 (NCBI 37) using Eland extended (Illumina, Inc.) or Bowtie (45), were carried out at the BioMicro Center at the Massachusetts Institute of Technology (MIT; Cambridge, MA). Raw sequencing reads were obtained for individual livers (n = 3 to 4 biological replicates/sex/mark, except for K9me3 and K36me3 at n = 2 replicates). All reads were trimmed to 35 nt to maintain a constant sequence length between samples. Replicates were evaluated by the percentage of reads in straight peaks (more than five identical reads that do not overlap any other reads and are therefore likely artifactual), by the degree of overlap between peaks detected in each replicate, and by correlation in read count between replicates (see Table S7 and “Analysis of data quality” in the supplemental material). Samples with >3% of the reads in straight peaks generally showed low overlap between biological replicates and were excluded. The final data sets for each chromatin mark ranged from 10.7 million (K9me3) to 52.5 million (K4me1) total reads (see Table S3A in the supplemental material).
To identify peaks and regions of chromatin mark enrichment, each data set, after combining data for biological replicates, was analyzed separately for male and female liver, as follows. K4me1, K27ac, and K4me3 form focal (localized) peaks and were analyzed using MACS, version 1.4.1 (46), with default parameters. K27me3, K9me3, and K36me3 were analyzed using SICER, version 1.1 (47), which identifies broad genomic regions (i.e., broad domains or islands). For K27me3, a window size of 400 bp was used, and a gap size of 2,400 was chosen as most appropriate (see Fig. S12 in the supplemental material). The same parameters were used for K9me3. For K36me3, whose peak regions were not as broad (see Table S3A), a window size of 200 bp and a gap size of 800 bp were used. Sex-enriched peaks and regions were identified for each chromatin mark, as follows. For each MACS peak and SICER region, the number of male sequence reads was compared to the number of female sequence reads after normalization by the total number of reads in common peaks—peaks (or regions) that were detected by MACS (or by SICER) in both male and female liver—in order to account for any bias between male and female livers in the percentage of reads in peaks/regions. The log2 male/female ratio (M value) was then calculated, along with a P value based on a Bayesian model (48) as implemented by Shao et al. (49). Each chromatin mark (focal peak or broad genomic region) was classified as male enriched (M > 1 and P < 0.001), female enriched (M < −1 and P < 0.001), or sex independent (−1 < M < 1 or P > 0.001). Table S3 in the supplemental material lists the numbers of peaks and domains identified for each chromatin mark along with their full data sets and genomic coordinates.
The six chromatin mark data sets along with liver DHS data reported previously (39) were analyzed using ChromHMM (1) to learn a hidden Markov model and to assign chromatin states across the mouse genome. ChromHMM was run using an IgG control and with default parameters. A single joint model was learned using data from both male and female liver. A 15-state model was previously found to be biologically meaningful and consistently recovered across cell types for the human genome (1). We therefore started with 20 states and used the ChromHMM CompareModels module to compare decreasing numbers of states to the 20-state model. We then calculated, for each of the 20 states, the similarity—i.e., the correlation between emission parameters—to its closest state in smaller models (see Fig. S3B in the supplemental material). A 14-state model was chosen as the point after which any further decrease in the number of states in the model resulted in states from the 20-state model being recovered with decreasing similarity (see Fig. S3B and C and “Chromatin states in mouse liver” in the supplemental material).
Analysis using high-throughput sequencing of RNA transcripts of total mouse liver (RNA-seq) was carried out using two independent pools of male and two independent pools of female mouse liver RNA (n = 6 livers/pool). Poly(A)-enriched RNA was used to prepare RNA strand-specific libraries (50) for Illumina sequencing at the BioMicro Center at MIT (Cambridge, MA). Forty-nucleotide sequence tags were mapped to the mm9 genome using Bowtie (45) and mapped to Mouse Genome Informatics (MGI) genes (www.informatics.jax.org/) using SeqMonk (Babraham Bioinformatics, Cambridge, United Kingdom [www.bioinformatics.bbsrc.ac.uk/projects/seqmonk/]). The set of all genes represented in this data set was clustered by their read densities in the 2 kb surrounding the TSS and 2 kb surrounding the transcription end sites (TES), with clustering carried out separately in male and female liver. A data set of 15,533 genes was obtained after the exclusion of short genes (<5 kb in length) to eliminate genes containing TSS-associated sequences that overlap the TES-associated region; 7,864 of the 15,533 genes were defined as liver expressed based on log2 reads per kilobase per million mapped (RPKM) was ≥1 in either male or female liver. At this RPKM but not at lower RPKM values, ~80% of genes contained a K4me3 peak and ~90% of genes overlapped with a K27ac island. For each gene, the region of the TSS ± 1 kb and of the TES ± 1 kb was subdivided into a total of 20 200-bp windows. Read counts were obtained for these 20 windows for seven features—six chromatin marks and DNase hypersensitivity—yielding 140 columns of data. Read counts were normalized by total sequencing library size for each ChIP-seq data set, and log transformed. The 15,533 genes were then clustered by their read densities by k-means clustering using Cluster (51), with a k value of 6. The six clusters of genes were visualized as a heat map (see Fig. S5A in the supplemental material) using Java Treeview (52). Each cluster was characterized as active or inactive based on the relative density of activating chromatin marks (K4me1, K27ac, K4me3, and K36me3) and repressive chromatin marks (K27me3 and K9me3). The frequencies with which each 200-bp bin in the TSS and TES regions, expanded to ±2 kb, were assigned 1 of the 14 chromatin states were calculated for genes in each cluster. Three of the six clusters were identified as active based on high read density of K27ac, K4me3, K36me3, and DNase hypersensitivity. Two clusters were designated poised based on the presence of K27me3 together with K4me1, with the poised chromatin state 12 at the TSS. The sixth cluster was inactive, characterized by K27me3 and low levels of activating marks.
In other analyses, sex-biased genes were clustered by their chromatin mark densities. Sex-biased genes were identified from the RNA-seq data set using edgeR (53) for differential expression analysis, as implemented by P.-Y. Hao and T. Melia of this laboratory. A total of 423 male-biased and 477 female-biased genes (listed in Table S1B and C in the supplemental material) were selected based on the following criteria: adjusted P value of <0.01 for sex difference, corresponding to a minimum sex differences in expression of 1.4-fold; liver expressed at a level of log2 RPKM of ≥1 in either sex; genes located on autosomes; and gene length of ≥2 kb. Each set of sex-biased genes was clustered by read densities in male liver and separately in female liver for the 2-kb regions surrounding the TSS and the TES, as described above, by k-means clustering using a k value of 3. Each of the clusters was defined by its chromatin activity—high, medium, and low—as indicated by the relative abundance of activating and repressive chromatin marks (see Fig. S6A in the supplemental material). Genes were then classified into six classes according to the cluster they fell into in male liver and in female liver (see Table S4A in the supplemental material). Each of the six gene classes was then characterized by mapping each gene to chromatin modifications identified by MACS (46) or SICER (47), which were subjected to a direct male-to-female comparison as described above (“Analysis of ChIP-seq data for chromatin modifications”). K4me3 marks were mapped to the promoter (TSS ± 500 bp), and K27ac marks were mapped to the promoter (TSS ± 500 bp) and separately to the gene body (from TSS up to TES); the other four marks were mapped to the gene body alone. Genes that contained both a male-enriched mark and a female-enriched mark for a particular modification were assigned to the sex-independent mark class. Box plots were generated using the box plot function of R with default parameters.
Gene sets that respond to activators of the nuclear receptors constitutive androstane receptor (CAR) (54, 55), pregnane X receptor (PXR) (54), and peroxisome proliferator-activated receptor alpha (PPARα) and PPARβ/δ (56) and to activators of AhR (57) were obtained from the indicated references. For CAR, the union of CAR-upregulated genes from Tojima et al. (54) and Tian et al. (55) was considered, with a corresponding set used for CAR downregulated genes. Genes up- or downregulated by each of the receptors were tested for enrichment for one of the six clusters (see Fig. S5 in the supplemental material), compared to the background set of all 15,533 genes (as defined above), using Fisher's exact test. For enrichment analysis of sex-biased genes shown in Fig. 5, genes were mapped to ChIP-seq-determined binding sites for STAT5, BCL6 (30), CUX2 (34), FOXA1, or FOXA2 (44). For each set of TFBS, the nearest gene body within 10 kb of a TFBS were defined as the target of the corresponding TF. Gene targets of these TFs were then tested using Fisher's exact test for enrichment for belonging to each of the classes of sex-biased genes (Fig. 5A and andB;B; see also Table S4B and C in the supplemental material) compared to a background set comprised of 8,764 liver-expressed genes (log2 PRKM of >1 in either sex) that were >2 kb in length. For each TF the fold enrichment was calculated as (number of genes in each sex-biased gene class that are a TF target/total number of sex-biased genes in the gene class)/(number of all other liver-expressed genes that are a TF target/number of all other liver-expressed genes). Functional term enrichment analysis was performed using DAVID (58, 59). FOXA1 and FOXA2 ChIP-seq data from Li et al. (44) were analyzed by G. Bonilla of this laboratory to identify sex-biased FOXA binding sites using methods described elsewhere (30).
The set of 72,862 merged liver DHS identified previously (39) was ranked by the male-female ratio after normalization by reads in male-female liver common peaks. The DHS genomic regions were ranked separately by their male-female ratio in DNase hypersensitivity and in K4me1, K27ac, and K27me3 read density over the entire peak region. Overlaps between DHS and TF binding sites (TFBSs) (male-enriched STAT5, female-enriched STAT5, male-enriched FOXA1 and FOXA2, and female-enriched FOXA2 binding sites) were also computed. Female-enriched FOXA1 binding sites did not follow the patterns exhibited by the other sex-biased TF binding sites examined (see Fig. S8F and S9B in the supplemental material) and were excluded from subsequent analyses. A TF binding site was considered to overlap a DHS if the ChIP-seq-identified peak region for the TF overlapped the DHS peak region by at least one base pair. Gene set enrichment analysis (60) was used to calculate a running enrichment score for each of the five types of TF binding sites for each of the four ranked lists: sites ranked by sex ratio in DHS, K4me1, K27ac, and K27me3.
Male-biased, female-biased, and sex-independent DHS (39) were characterized by the presence or absence of the two enhancer-associated chromatin marks, K4me1 and K27ac, and by whether the marks were classified as male enriched or female enriched (|M| > 1 and P < 0.001) (see above) or sex independent. For each DHS, the DHS peak region was defined as a 1-kb window centered at the DHS peak summit. Each DHS was considered to contain a K4me1 or K27ac peak if at least 200 bp of the 1-kb DHS peak region overlapped a K4me1 or a K27ac peak. The requirement for a 200-bp overlap was chosen based on the K4me1 and K27ac read profiles at DHS, which peak at ~300 bp on either side of the DHS summit (see Fig. 2E).
Since K27ac marks identify active enhancers and since K4me1 marks without K27ac marks identify enhancers that may not be active (61), DHS were assigned enhancer categories in a hierarchical manner, in the following order (see Table S6A in the supplemental material): (i) sex-biased K27ac, (ii) sex-biased K4me1, (iii) sex-independent K27ac, (iv) sex-independent K4me1, and (v) no K27ac or K4me1 mark. Any DHS that overlapped both a male-biased mark and a female-biased mark for a particular modification was considered to overlap a sex-independent mark. Each of these categories of DHS was compared to the same background set of sex-independent DHS (58,087 DHS in male liver; 52,187 DHS in female liver). Background sex-independent DHS were those whose nearest gene TSS was not sex biased, and additionally, they were required to be >500 kb from the nearest sex-biased TSS. DHS were considered to contain a STAT5, BCL6, CUX2, FOXA1, or FOXA2 binding site (identified by ChIP-seq) if the DHS overlapped the TF's binding site by at least 1 bp. For each category of DHS by enhancer status, the enrichment for a particular TF was calculated as follows: (number of DHS in a category that overlap a TFBS/total number of DHS in that category)/(number of background sex-independent DHS that overlap a TFBS/total number of sex-independent DHS). The Fisher exact test was used to calculate a P value for each enrichment or depletion. Enrichments having a P value of <0.001 are shown in Fig. 6B.
Raw sequence reads for the six sets of chromatin marks in both male and female mouse liver are available at the Gene Expression Omnibus (GEO) website as series GSE44571 (samples GSM1087069 to GSM1087105 and GSM1181495), and sequence reads for RNA-seq data for gene expression in male and female mouse liver are available as series GSE48109 (samples GSM1168535 to GSM1168538). Mouse liver DNase-seq data sets (39) are available as series GSE21777, and ChIP-seq data sets are available for STAT5 (GSE31578; male and female liver, STAT5 high-activity status), BCL6 (GSE31578; male liver) (30), CUX2 (GSE35985; female liver) (34), and FOXA1 and FOXA2 (44) (untreated male and female liver; E-MTAB-805) in ArrayExpress (www.ebi.ac.uk/arrayexpress/).
Six chromatin marks were mapped genome-wide in male and female mouse liver. All four activating marks (histone H3 K4me1, K4me3, K27ac, and K36me3) showed the anticipated positive correlations with each other and with expression of associated genes, and both repressive marks (histone H3 K9me3 and K27me3) were negatively correlated with gene expression (Fig. 1A; see also Fig. S1A in the supplemental material). Genomic localizations and distributions across gene bodies (see Fig. S2 in the supplemental material) are consistent with previously reported associations of each mark: K4me3 with active promoters (62–64), K36me3 with actively transcribed regions (62, 63, 65), K4me1 and K27ac with enhancers (61, 66), K27me3 with polycomb-mediated gene silencing, and K9me3 with constitutive heterochromatin and nongenic regions (67). K27me3 and K9me3 were enriched in female liver on chromosome X (ChrX) (Fig. 1B; see also Fig. S1B), consistent with their role in X inactivation (68). Conversely, K4me1 and K27ac are enriched in male liver on ChrX, consistent with depletion of K4 methylation from the inactive X chromosome in female cells (69). ChrX was excluded from further analysis to eliminate bias in genome-wide male-female comparisons.
ChromHMM (1) was used to learn 14 chromatin states in male and female liver based on DHS data (39) combined with the above six chromatin marks (Fig. 2A, and andB;B; see also Fig. S3A in the supplemental material). These include inactive states (marked by K27me3 [state 1] or K9me3 [state 3]), a bivalent promoter-associated state (K27me3 plus K4me1 [state 12]), a transcribed state (K36me3 [state 14]), promoter states, which contain K4me3 (states 7 and 8), and several enhancer states, which contain K27ac or K4me1 but lack K4me3 (states 5, 6, 9, 10, and 11). Liver-expressed genes show active promoter states at the transcription start site (TSS) (Fig. 2B, left), with an asymmetric pattern reflecting the asymmetric distribution of DHS and K4me1, K4me3, and K27ac marks at promoters (Fig. 2C, left) (5). Active chromatin states occur more frequently in one sex than the other at GH-regulated DHS whose accessibility differs significantly between male and female liver (sex-biased, GH-responsive DHS ) (Fig. 2D; see also Fig. S4B in the supplemental material) but not at sex-independent DHS (see Fig. S4A). K27me3 states are elevated in male compared to female liver at female-biased DHS (cf. male prominence of K27me3-containing states 1 and 12), whereas we did not find a corresponding female bias in K27me3 states at male-biased DHS (Fig. 2D and andE).E). This indicates a role for K27me3 in sex-biased gene silencing in male but not female liver. Overall, sex differences in chromatin environment are much less pronounced at sex-biased TSS than at sex-biased DHS (Fig. 2F versus versusE;E; see also Fig. S4C), suggesting that many sex-biased genes are controlled by distal regulatory sites.
Genes were clustered by chromatin mark and DHS read densities surrounding the TSS and transcript end site (TES) (see Fig. S5A and Table S1 [for the full data set for 15,533 genes] in the supplemental material). Chromatin state distributions for the resultant six gene clusters are shown in Fig. 3A. Clusters 1 to 3 were characterized by high levels of DHS and the presence of all four activating marks and by low levels of repressive marks and are primarily in chromatin state 7 at the TSS (Fig. 3A). Genes in these active chromatin state clusters show high expression (≥80% at log2 RPKM ≥ 1) (Fig. 3B; see also Fig. S5B). Cluster 1 differs from clusters 2 and 3 in its active enhancer states around the TES (Fig. 3A). Cluster 1 genes also have the highest coverage by broad domains of DNase hypersensitivity, with 70% of the genes being at least half covered by an extended DHS island (39), compared to 16 to 23% for the other two active clusters (see Fig. S5C and D). Genes in clusters 4 to 6 have high K27me3 read densities and, correspondingly, show low expression (Fig. 3B; see also Fig. S5B). Genes in clusters 4 and 5 are characterized by K4me1 in combination with K27me3 (chromatin state 12 at the TSS) (Fig. 3A), i.e., the signature of a poised gene (70), while cluster 6 genes have K27me3 without K4me1 or other activating marks (state 1) (Fig. 3A). The poised-gene clusters show significantly higher mean expression levels than cluster 6 genes (Fig. 3B).
To test the hypothesis that the basal chromatin state informs responsiveness to transcriptional regulators, genes in all six chromatin clusters were compared to the sets of genes that respond to four ligand-activated TFs of the nuclear receptor superfamily: CAR, PXR, PPARα, and PPARβ/δ (71). Genes regulated by these TFs were enriched for one or more active chromatin clusters (Fig. 3C and andD;D; see also Table S1 in the supplemental material). Further, genes downregulated by PXR and PPARβ/δ were enriched for poised chromatin clusters, suggesting that ligands that activate these two TFs induce removal of K4me1 while maintaining K27me3 at target genes. In contrast, genes upregulated by another ligand-activated TF, AhR, were enriched for genes in a poised chromatin state, indicating loss of K27me3 upon ligand activation. Consistent with the above preference for activation of genes already in an active chromatin state, the upregulated gene targets of all four nuclear receptors were significantly depleted from the inactive chromatin state (cluster 6) (Fig. 3C). Thus, the basal chromatin state is apparently a determinant of responsiveness to these ligand-activated TFs.
Functional term enrichment analysis (see Table S2 in the supplemental material) indicated enrichment for liver-specific terms such as drug, lipid, and steroid metabolism in cluster 1, the most active gene cluster. The other two active clusters were enriched for a broad range of non-liver-specific terms, e.g., mitochondrion, RNA processing, protein transport, and cell cycle. The two poised clusters were enriched for terms such as cell adhesion, cell junction, mesenchymal cell differentiation, and extracellular matrix, while the inactive cluster 6 was enriched for genes not expressed in liver (e.g., neurotransmission, muscle, vision, and memory).
Genes that show sex-biased expression (423 male-biased genes and 477 female-biased genes) were clustered by their chromatin mark and DHS densities around the TSS and TES in male liver and, separately, in female liver (see “Classification of sex-biased genes” and also Fig. S6A to C in the supplemental material). Six female-biased and six male-biased gene classes (F1 to F6 and M1 to M6, respectively) were identified based on their chromatin activity classifications in each sex (see Table S4A and Fig. S6A in the supplemental material). Next, sex differences in the chromatin environment of each gene class were characterized by comparing normalized densities of each of the six marks on a genome-wide basis in male versus female liver. Genomic regions that showed significant male enrichment or female enrichment were thus identified for each chromatin mark (see Table S3 in the supplemental material). The frequency of sex-biased and sex-independent chromatin marks was then determined for the genes in each class. Figure 4A presents these data for the two largest classes of sex-biased genes, which, as shown, generally lack sex differences in proximal chromatin marks (classes F1 and M1), as do a majority of all sex-biased genes. In contrast, genes in sex-biased gene classes F3, M3, and M4 show sex differences in proximal chromatin marks at comparatively high frequencies (Fig. 4A; see the full data set in Fig. S6E and Table S1B and C in the supplemental material). Classes F3 and M3 also showed the largest sex differences in expression (Fig. 4B), while classes F3 and M4 showed the highest level of expression (see Fig. S6G). Heat maps of chromatin mark sex ratios are shown for classes F3, M3, and M4 in Fig. 4C. These three classes (F3, M3, and M4) represent only 5.6% of all sex-biased genes, consistent with the conclusion that many sex-biased genes are not in a sex-biased chromatin environment (cf. Fig. 2F). Indeed, taking into account distal chromatin marks, which may represent distal regulatory sites, less than half of sex-biased genes have sex-biased chromatin marks within 10 kb, and only 66 to 69% have them within 100 kb (see Fig. S6F).
Importantly, the highly female-biased class F3 genes are characterized by female-biased activating marks and male-biased repressive marks (note the high frequency of male-enriched K27me3), whereas class M3 and M4 genes frequently show male-biased activating marks but not female-biased repressive marks (Fig. 4A). Our finding of male-enriched K27me3 marks at highly female-biased genes but not the corresponding female-enriched K27me3 marks at highly male-biased genes is supported by a separate analysis in which sex-biased genes were divided into four equal groups based on their ranking in sex ratio in expression (see Fig. S7A in the supplemental material). As shown in Fig. S7A, 30 out of the 118 most highly female-biased genes (25%) have male-enriched K27me3 marks, compared to ≤4% for each of the other quarters. In contrast, none of the male-biased genes have female-enriched K27me3 marks.
Next, we investigated whether sex-biased genes in different chromatin states utilize the same regulatory factors and mechanisms. To identify putative targets of such regulatory factors, the sex-biased genes in each class were mapped to nearby (within 10 kb) binding sites for five TFs implicated in GH-regulated, sex-dependent liver gene expression. Enrichment scores were then calculated for each TF compared to a background set of all liver-expressed genes. Figure 5A and andBB present these data for three GH-regulated TFs: STAT5 (male-enriched, female-enriched, and sex-independent STAT5 binding sites) (30), BCL6, a male-biased repressor (30), and CUX2, a highly female-specific repressor (34). Data are also shown for FOXA1 and FOXA2, which contribute to the sex-biased incidence of hepatocellular carcinoma (44). Several of the sex-biased gene classes were enriched for being targets of sex-biased STAT5 binding, with the highest enrichments seen for genes in classes F3 and M4, whose proximal chromatin marks show sex bias (Fig. 4A and andC).C). In contrast, BCL6 targets were enriched among female-biased genes that do not show sex differences in chromatin marks (classes F1, F4, and F6), suggesting that the repression of female-biased genes in male liver by BCL6 (30) does not introduce (or require) a sex-biased chromatin state. BCL6 targets were also enriched at a subset of male-biased genes, reflecting the STAT5 competition mechanism described earlier (30).
CUX2 targets were highly enriched among male-biased genes in class M4 and, to a lesser extent, in class M1 (Fig. 5A), consistent with the proposed role for CUX2 as a repressor of many male-biased genes in female liver (34). CUX2 targets were also enriched among class F3 genes, consistent with CUX2 directly activating a subset of female-biased genes (34). Female-biased CUX2 targets are more likely than other female-biased genes to have a male-enriched K27me3 mark, especially the target genes in class F3 (Fig. 5C), which show the greatest enrichment of CUX2 binding (Fig. 5A) and the highest frequency of male-biased K27me3 marks (Fig. 4A). Since class F3 genes are also highly enriched for being targets of STAT5 and FOXA2 (Fig. 5A), these three factors may cooperatively regulate the most highly female-biased genes. Overall, male-enriched K27me3 domains at the gene body of female-biased genes coassociate with CUX2 binding (within 10 kb) more than with any other TF examined, and, correspondingly, CUX2 binding is more likely to target female-biased genes with male-enriched K27me3 domains than any other sex-biased mark (see Table S4D and E in the supplemental material). FOXA2 shows female-biased binding enriched at female-biased genes with male-enriched K27me3 domains, and this enrichment is almost as high as the corresponding enrichment for CUX2 (see Table S4D and E), suggesting that CUX2 and FOXA2 cooperate to derepress these genes in female liver.
Among female-biased genes, the preference for STAT5 and FOXA2 to target genes that have proximal sex-biased marks and for BCL6 to target genes that lack sex-biased marks is further supported by enrichment analysis of female-biased genes ranked by gene expression sex ratio (see Fig. S7B in the supplemental material). The enrichment for being a target of female-enriched STAT5 or FOXA2 binding increases with an increasing sex ratio in gene expression. In contrast, BCL6 targets are enriched only among the female-biased genes with the lowest sex ratio in expression. Thus, female-enriched STAT5 and FOXA2 binding sites preferentially map to the most highly female-biased genes (Fig. 5A; see also Fig. S7B), which in turn exhibit sex differences in proximal chromatin marks (Fig. 4A and andB),B), while BCL6 targets weakly female-biased genes, which lack sex differences in proximal chromatin marks.
Female-biased DHS (39) are highly enriched (18-fold), compared to sex-independent DHS, for the presence of female-biased K27ac or K4me1 marks; correspondingly, a 16-fold enrichment of a male-enriched K27ac or K4me1 mark is seen at male-biased DHS (see Table S5A in the supplemental material). Thus, many sex-biased DHS have the marks of sex-dependent enhancers. Moreover, sex-biased gene classes F3 and M3, which comprise the most highly sex-biased genes, are enriched for association (within 250 kb) with sex-biased DHS that have sex-biased K27ac (see Table S5B and C), the mark of an active enhancer (61, 66).
Next, we investigated the relationship between sex bias in chromatin marks and DNase hypersensitivity and sex-biased binding of the five TFs discussed above. Gene set enrichment analysis showed positive correlations between the sex-biased binding of STAT5, FOXA1, and FOXA2 and the ranking by sex ratio of each enhancer mark (K27ac and K4me1) and of DNase hypersensitivity and inverse correlations in the case of K27me3 (Fig. 6A). These findings are supported by our analysis of TF binding site frequency in DHS ranked by sex ratio (see Fig. S8A to E in the supplemental material) and by the correlation between sex ratio in TF binding and sex ratio in DHS and in activating chromatin marks (see Fig. S9A to C in the supplemental material). BCL6 is most frequently bound at sex-independent DHS (see Fig. S8G), and BCL6 binding intensity does not correlate with sex ratio in chromatin marks (see Fig. S9D), consistent with its preferred sex-biased target genes (Fig. 5A) having sex-independent chromatin marks (Fig. 4A). CUX2 binding is enriched at male-biased DHS, consistent with the study of Conforto et al. (34), but is not associated with or correlated with sex differences in chromatin marks at the sites where it is bound (see Fig. S8H and S9E).
To examine more closely the relationship between chromatin modifications and regulation by liver TFs, we integrated TF binding data with the enhancer status of DHS and the chromatin class of their respective gene targets, considering genes within 250 kb to allow for distal regulation. For each TF, we identified significant enrichments (P < 0.001) of its binding sites at subsets of DHS categorized by their enhancer modification patterns (Fig. 6B). Sex-biased DHS that have sex-biased enhancer marks are most highly enriched for being GH responsive, linking sexually dimorphic plasma GH profiles and sex differences in enhancer modifications. Consistent with the results above, female-biased binding of STAT5 and of FOXA2 is most highly enriched at female-biased DHS that have female-biased K27ac marks. Similarly, male-biased STAT5, FOXA1, and FOXA2 binding is most highly enriched at male-biased DHS with male-biased K27ac. This preference for sites with sex-biased K27ac is also seen at sex-independent DHS but with lower enrichments (Fig. 6B). Enrichments of TF binding sites were also determined for subsets of DHS that target (within 250 kb) genes in individual sex-biased gene classes (see Table S6B to D in the supplemental material). Female-biased STAT5 binding was most strongly enriched (22-fold) at DHS mapping to class F3 genes, consistent with the results shown in Fig. 5A. In contrast, the strongest enrichment for male-biased STAT5 binding was at DHS that map to class M1 genes (47-fold), which lack proximal sex-biased chromatin marks (Fig. 4A). Thus, distal regulation by male-enriched STAT5 binding is associated with male-biased genes that lack nearby sex differences at the chromatin level.
BCL6 and CUX2 binding are enriched at both male-biased and sex-independent DHS, in particular for sites with a sex-independent K27ac mark (Fig. 6B), in agreement with the lack of correlation between binding sites for these factors and sex differences in chromatin marks described above. Among DHS with sex-independent K27ac, BCL6 binding is enriched in association with female-biased genes that lack nearby sex-biased chromatin marks (classes F1, F2, F4, and F6), while CUX2 binding is most highly enriched in association with male-biased class M4 genes (7.8-fold) (see Table S6C and D in the supplemental material), in agreement with the results shown in Fig. 5A and andB.B. Bcl6 itself is a target of male-biased STAT5 binding (see Table S1C in the supplemental material) and is an example of a male-biased TF that binds to sex-independent regulatory sites, which enables it to regulate secondary target genes that do not exhibit sex differences in chromatin marks. CUX2 also shows enrichment for binding at sex-independent DHS with female-biased K27ac marks (Fig. 6B). The highest CUX2 enrichment was associated with female class F2 genes (which includes Cux2 itself) and class F3 genes (see Table S6B and D). CUX2 is thus a female-specific TF that regulates female-specific genes via sex-independent DHS regulatory elements.
FOXA1 acts as a liver-enriched pioneer factor (11, 12) that opens chromatin at cell-type-specific enhancers by binding compacted chromatin via a core histone binding motif, resulting in the generation of open chromatin regions (DHS) (72, 73). Given the strong enrichment of sex-biased binding of FOXA1 and of the related FOXA2 at DHS showing a sex bias in DNase hypersensitivity and K27ac and K4me1 marks (Fig. 6), we hypothesized that FOXA1 and/or FOXA2 mediates sex-dependent nucleosome repositioning and sex-differential chromatin opening at sex-biased DHS. Supporting this hypothesis, male-biased DHS show a sex-dependent distribution of K4me1 marks surrounding the DHS summit, with a bimodal peak in male liver and a monomodal peak in female liver (Fig. 7A). This pattern is consistent with sex-dependent nucleosome positioning (13) and is most pronounced at DHS that bind FOXA1 or FOXA2 in a male-enriched manner (cf. the greater dip in K4me1 reads at the summit of male-biased DHS bound by FOXA1 or FOXA2 compared to male-biased DHS not bound by these factors, as shown in Fig. 7B). Moreover, in contrast to the monomodal K4me1 peak seen in female liver at DHS with male-enriched FOXA1 or FOXA2 binding, a trough in K4me1 reads is discernible in female liver at DHS with sex-independent FOXA1 or FOXA2 binding. Further support for sex-dependent nucleosome repositioning at male-biased DHS by FOXA1/FOXA2 is provided by the deeper K4me1 trough at male-biased DHS in male liver when FOXA1 or FOXA2 binding occurs in combination with male-biased STAT5 binding than at sites bound by STAT5 alone (Fig. 7C; see also Fig. S10A and B in the supplemental material).
In contrast to male-biased DHS, DHS that are female biased or sex independent have bimodal K4me1 peaks in both male and female liver (Fig. 7A). Female-biased DHS show no sex difference in K4me1 profiles at sites of female-enriched FOXA2 binding (Fig. 7D), despite the strong enrichment of female-biased FOXA2 binding at female-biased DHS (Fig. 6) and at female-biased genes (Fig. 5A; see also Fig. S8D and S9C and Table S6B in the supplemental material). Thus, while FOXA2 is likely involved in female-biased gene regulation, it apparently does not confer sex-dependent nucleosome repositioning at female-biased DHS. The finding that the FOXA factor(s) is apparently active in nucleosome repositioning at male-biased but not female-biased DHS can be explained if FOXA1, rather than FOXA2, is the key factor in sex-dependent chromatin opening, as suggested by the deficiency of female-enriched compared to male-enriched FOXA1 binding sites (577 and 3,280 sites, respectively) and by the more similar numbers of female-enriched and male-enriched FOXA2 binding sites (780 and 976 sites, respectively). Supporting this proposal, 71% of the male-enriched FOXA2 binding sites found at male-biased DHS cooccur with a male-enriched FOXA1 binding site, whereas only 6% of female-enriched FOXA2 binding sites cooccur with a female-enriched FOXA1 binding site. Moreover, the female-enriched FOXA1 binding sites that are present are unusual insofar as they do not preferentially overlap, nor do they correlate with the magnitude of the sex bias in female-biased DHS and chromatin modifications (see Fig. S8F and S9B). To summarize, the pronounced sex-dependent trough of K4me1 reads at male-biased DHS is indicative of nucleosome repositioning occurring in male but not female liver. This sex-differential pattern of marks is associated with, and may be regulated by, male-enriched FOXA1 binding, although it often coincides with male-enriched FOXA2 binding. The relative deficiency of female-enriched FOXA1 binding sites may account for the absence of a corresponding strong, female-specific trough of K4me1 reads at female-biased DHS.
We have developed genome-wide chromatin state maps for male and female mouse liver based on DHS profiles and patterns of six chromatin marks and then used these maps to elucidate chromatin environments at DHS and across gene bodies, in particular, those associated with sex differences in gene expression, a characteristic of ~1,000 genes in mouse liver. By integrating chromatin maps with TF binding data, we show that the mechanisms of sex-biased gene regulation differ between male and female liver, with male-biased DHS but not female-biased DHS being characterized by a sexually dimorphic K4me1 profile that is most closely associated with male-enriched binding of FOXA pioneer factors (Fig. 8). These FOXA factors, in turn, are proposed to facilitate male-enriched binding of STAT5, which activates male-biased genes by mechanisms that include distal regulation of the male-biased genes that lack nearby sex differences in chromatin marks. Genes that show a weak female bias in expression and are associated with sex-independent chromatin marks are preferentially repressed in male liver by the male-biased repressor BCL6, while other, more highly female-biased genes show sex differences in proximal chromatin marks and are preferentially activated by female-enriched STAT5 binding. Many of the most strongly female-biased genes are regulated by K27me3-based repression in male liver, which may be countered by the stimulatory actions of CUX2 in female liver; however, strongly male-biased genes are not repressed by a corresponding K27me3-based mechanism in female liver (Fig. 8). These insights evidence the power of functional genomic studies and epigenetic analysis for elucidating transcriptional regulatory networks in a complex mammalian tissue.
Clustering of genes by local chromatin mark intensities and DNase hypersensitivity identified genes in three distinct chromatin environments: active (clusters 1 to 3), poised (K27me3 plus K4me1 state at the TSS; clusters 4 and 5), and inactive (low activating marks, high K27me3 levels; cluster 6) (Fig. 3). Cluster 1 genes are commonly involved in liver-specific tasks, such as drug, steroid, and lipid metabolism; genes in clusters 2 and 3 perform many non-liver-specific functions, and poised cluster 4 and 5 genes encompass biological processes activated in response to specific stimuli. Genes responsive to nuclear receptors/TFs that regulate many liver metabolic functions (CAR, PXR, PPARα, and PPARβ/δ) (41–43) and genes responsive to AhR were depleted of inactive cluster genes, indicating that these TFs do not target genes in an inactive chromatin state. A poised chromatin state is preferentially associated with genes downregulated by PXR and PPARβ/δ and with genes upregulated by AhR, indicating that poised chromatin can resolve into either an inactive or an active state, depending on the TF and depending on the target gene. In contrast, many genes upregulated by PXR and PPARβ/δ, as well as genes responsive to CAR and PPARα, are already in an active chromatin state in unstimulated liver.
Sex differences in chromatin environment were commonly seen at male-biased and female-biased DHS but not at sex-biased genes, a majority of which lack sex differences in local chromatin marks. Consistently, genes differentially expressed between cell types show cell-type-specific chromatin states at regulatory sites but common chromatin states at the promoter (3, 6, 74). However, some sex-biased genes showing large differences in expression between sexes do show notable sex differences in local chromatin states (e.g., Cyp2b9 and Hsd3b5) (see Fig. S11 in the supplemental material) while others do not (Sult2a6 and Cabyr) (see Fig. S11). Genes with sex differences in local chromatin environment show the strongest enrichment for local sex-biased STAT5 binding (Fig. 5A and andB).B). Distal regulation is most likely for class M1 genes (the largest class of male-biased genes), which generally lack sex differences in local chromatin marks but show 47-fold enrichment for distal (within 250 kb) male-biased STAT5 binding sites at male-biased DHS with male-biased K27ac marks (see Table S6C in the supplemental material). Globally, sex-biased STAT5 binding at sex-biased DHS is 3.6-fold more common in male than female liver (male-biased STAT5 binding at 624 of 2,714 male-biased DHS versus female-biased STAT5 binding at 173 of 1,333 female-biased DHS), which may enable STAT5 to distally regulate a larger number of sex-biased genes in male liver. Finally, sex-biased STAT5 binding is most highly enriched in association with sex-biased K27ac, the mark of an active enhancer, both at sex-biased and at sex-independent DHS (Fig. 6B). This finding suggests that many of the sex-dependent actions of STAT5 action proceed via sex-biased enhancer sequences and highlights the importance of sex differences in chromatin state for sex-dependent gene expression. A failure to recapitulate a sex-dependent chromatin environment may explain the inability of some sex-biased DHS sequences to confer sex-differential reporter gene expression when they are assayed for enhancer activity in liver in vivo (39).
Our analysis of sex-differentiated chromatin states in combination with binding data for TFs implicated in sex-biased liver gene regulation supports the model depicted in Fig. 8, which shows proposed regulatory events and relationships for the largest classes of sex-biased genes with (F3 and M4) or without (F1 and M1) sex differences in proximal chromatin marks. Plasma GH profiles, which are sexually dimorphic, impart sex-differential patterns of STAT5 activation (75, 76) and sex-differential expression of Bcl6 (32) and Cux2 (33) and regulate DHS opening and closing in a sex-biased manner (39). This sex-dependent and plasma GH pattern-dependent remodeling of regulatory sites involves the enhancer-associated modifications K4me1 and K27ac, as evidenced by the strong enrichment for GH responsiveness among DHS that show sex bias in these two activating marks (Fig. 6B). The pioneer factor FOXA1, which initiates chromatin opening that activates enhancers during liver development (77), and the related FOXA2 show sex-dependent binding that correlates closely with sex bias in DHS and in K4me1 and K27ac marks (Fig. 6; see also Fig. S8 and S9 in the supplemental material) and are proposed to induce the observed sex-differential opening at male-biased DHS in mouse liver. Supporting this proposal, male-biased DHS show a sexually dimorphic K4me1 distribution around the DHS summits that is FOXA dependent (Fig. 7). This pattern is reminiscent of unstimulated LNCaP cells, where androgen receptor binding sites lacking FOXA1 binding have a broad K4me2 peak compared to a bimodal peak surrounding the summit of androgen receptor binding sites associated with FOXA1 binding (8). Furthermore, cell-type-specific binding by FOXA1 is enriched at sites with cell-type-specific K4me1 and K4me2 marks, which are required for FOXA1 binding at enhancers (72). Thus, FOXA1 is recruited to sites marked by cell-type-specific K4 methylation, which it utilizes as an epigenetic signature for chromatin opening. By extension, our finding that FOXA1 and FOXA2 are recruited in a male-biased manner to sites with male-enriched K4me1 marks in a bimodal distribution suggests that FOXA induces repositioning of nucleosomes to flank the site, thereby increasing DNase hypersensitivity in male liver at the peak center. Since 71% of the male-enriched FOXA2 binding sites at male-biased DHS overlap a male-enriched FOXA1 binding site, FOXA1 could be the key factor for sex-dependent chromatin opening. Indeed, at female-biased DHS, where FOXA2, but not FOXA1, commonly exhibits female-enriched binding, sex differences in the pattern of K4me1 read density are not apparent (Fig. 7D).
Both bimodal and monomodal patterns of K4me1 marks are seen at FOXA2 and HNF4A binding sites in liver, but the bimodal loci are associated with more highly expressed genes and are more likely to confer liver specificity (13). Further, K4me1 peaks that are bimodal in liver and monomodal in islets are associated with liver-specific genes, and a corresponding relationship holds for bimodal peaks in islets and islet-specific genes (13). The bimodal K4me1 peak is associated with nucleosome depletion at the peak center, as confirmed by nucleosome positioning analysis at a few representative loci (13). In liver, bimodal K4me1 loci occupied by FOXA2 are most likely to be cooccupied by HNF4A (13), and STAT5 binding sites in female liver show significant cooccupation by HNF4A, FOXA1, and FOXA2 (78). HNF4A is required for sex-biased gene expression and may act cooperatively with STAT5 to influence gene expression (35, 36). Thus, FOXA1, and perhaps FOXA2, when recruited to male-biased DHS marked by male-biased K4me1, may act in cooperation with HNF4A and/or STAT5 to regulate sex-biased genes in male liver. Of note, about half of all FOXA2 binding sites in mouse liver occur at nucleosome-bound sites, including some sites that exhibit bimodal K4me1 distribution (79). Consequently, the bimodal distribution of K4me1 cannot be taken as a definitive indication of nucleosome repositioning at every such site (79). Further studies are needed to identify sites of sex-dependent, GH-regulated nucleosome remodeling and to elucidate the individual roles of FOXA1 and FOXA2 in these remodeling events.
Sex-dependent nucleosome positioning at male-biased DHS but not female-biased DHS may explain our previous finding that male-enriched STAT5 binding sites are characterized by open chromatin in male liver versus relatively closed chromatin in female liver (30), which suggests that there is a specific, sex-dependent mechanism to open these sites in male liver. In contrast, female-enriched liver STAT5 binding sites are relatively open in both males and females (30) and therefore may not require sex-dependent chromatin opening by FOXA or other factors. Other mechanisms, such as the persistence of STAT5 binding activity in female liver (30), may be sufficient to confer female-biased gene expression for the target genes of female-enriched STAT5 binding sites.
BCL6, which is expressed in a male-biased manner subject to regulation by GH and STAT5, is a potent transcriptional repressor that competes with STAT5 for binding to chromatin and is associated with repression of many female-biased genes in mouse liver (30, 32). BCL6 preferentially binds sex-independent DHS with sex-independent chromatin marks (see Fig. S8G in the supplemental material), and the female-biased targets of BCL6 are enriched for genes that lack local sex-biased chromatin marks (Fig. 5A and and8B8B and andC;C; see also Table S6C and D in the supplemental material) and show weak sex bias (see Fig. S7B in the supplemental material). These findings support our earlier proposal that sex-biased genes deficient in nearby sex-biased DHS are in part regulated by TFs whose expression or activity is sex dependent (e.g., BCL6), without a need for sex bias in the chromatin environment at the TF binding site (39). The enrichment that BCL6 shows for weakly but not strongly female-biased targets may reflect the modest sex difference that characterizes BCL6 expression in mouse liver (32).
The high female specificity of CUX2 expression (~100-fold) (33) effectively restricts its binding to female liver chromatin, where CUX2 binding is enriched at male-biased DHS and at sites with male-enriched STAT5 binding, and is associated with repression of ~30% of male-biased genes (34). Unlike the female-biased targets of BCL6, male-biased genes that show the highest enrichment for being a proximal target of CUX2 (class M4 genes) are characterized by sex differences in local chromatin marks and also show the highest enrichment of male-biased binding of STAT5 and FOXA1 (Fig. 5B and and8A).8A). The male-gene-repressive actions of CUX2 in female liver are also associated with CUX2 binding at male-biased DHS and at sites with male-enriched STAT5 binding (34). This suggests that CUX2 may keep male-biased DHS closed in female liver and thereby contribute to the male bias in chromatin accessibility and STAT5 binding seen at those sites, which may be opened by FOXA1 in male liver. CUX2 can also activate female-biased genes in female mouse liver (34), and, consistently, CUX2 binding was enriched at both female-biased and sex-independent DHS that map to several classes of female-biased genes (see Table S6B and D in the supplemental material). Furthermore, female-biased targets of CUX2 are frequently marked by enrichment for K27me3 at the gene body in male liver (Fig. 5C). Thus, the activating role of CUX2 may involve chromatin remodeling at target genes—in this case, derepression of female-biased genes by removal of K27me3 repressive marks across the gene body (Fig. 8B). Presumably, the male bias in K27me3 marks at these genes is introduced at puberty, when the female-specific expression of CUX2 and a majority of its targets first emerge (80).
Comparison of the epigenetic environment and chromatin states that characterize male and female mouse liver and integration of these data with gene expression and TF binding data revealed distinct mechanisms of sex-biased gene regulation in male and female liver. Sex-dependent K27me3-mediated repression is shown to be an important mechanism of repression of highly female-biased but not of male-biased genes, and a sexually dimorphic K4me1 profile, which likely indicates sex-dependent nucleosome repositioning by FOXA pioneer factors, is seen at male-biased but not female-biased regulatory sites. TFs implicated in sex-biased gene expression were found to exhibit distinct relationships with epigenetic marks: STAT5-mediated activation was associated with sex-biased chromatin modifications, while BCL6-mediated repression was primarily associated with sex-independent chromatin modifications, both at BCL6 binding sites and at its target genes. Finally, the female-specific CUX2 may repress male-biased genes by interfering with STAT5 binding, while CUX2 activation of highly female-biased genes may involve removal of K27me3 repressive marks.
We thank the following lab members for their contributions to this study: Ekaterina Laz for ChIP sample preparation, Andy Rampersaud for confirmatory qPCR analysis of ChIP samples, Peng-Ying Hao and Tisha Melia for providing RNA-seq data, and Gracia Bonilla for initial analysis of sex-biased FOXA1 and FOXA2 binding sites.
Author contributions were as follows: A.S. and D.J.W. conceived and designed the study; A.S. carried out all of the computational studies, including analysis of ChIP-seq and DNase-seq data sets, and chromatin state analysis and prepared figures and tables for publication; A.S. and D.J.W. jointly wrote and revised the manuscript.
This work was supported in part by NIH grant DK33765 (to D.J.W.).
Published ahead of print 8 July 2013
Supplemental material for this article may be found at http://dx.doi.org/10.1128/MCB.00280-13.