Introduction

Growth hormone (GH) regulates hepatic expression of enzymes and transporters that play critical roles in lipid metabolism and in the detoxification of many drugs and other lipophilic foreign chemicals [1, 2]. Dysregulation of hepatic GH signaling can lead to liver metabolic disorders, including the development of fatty liver disease and non-alcoholic steatohepatitis, with males more susceptible than females, as seen in both mice and humans [3–5]. Correspondingly, GH regulates many liver-expressed genes in a sex-dependent manner, enabling each sex to meet its specific metabolic and hormonal requirements [6]. This program of sex-dependent gene expression is controlled by the sex-dependent temporal patterns of pituitary GH secretion [7, 8], which emerge at puberty but are programmed earlier in life by neonatal exposure to androgen [9–11]. Pituitary GH secretion, and consequently plasma GH profiles, are intermittent (pulsatile) in pubertal and adult males, whereas they are near-continuous (persistent) in pubertal and adult females, as seen in rats, mice and humans [12]. These plasma GH profiles, in turn, regulate sex-specific gene transcription in the liver through both positive regulatory mechanisms (class I sex-biased genes) and negative regulatory mechanisms (class II sex-biased genes) [6].

The sex-dependent hepatic actions of GH require the transcription factor STAT5 [13], which is activated by phosphorylation on a single tyrosine residue catalyzed by the GH receptor-associated tyrosine kinase JAK2 [14]. Each successive male plasma GH pulse induces a cycle of STAT5 tyrosine phosphorylation, dimerization and nuclear translocation, followed by tyrosine dephosphorylation and recycling of STAT5 back to the cytosol in time to reset the overall signaling pathway for the subsequent plasma GH pulse [1]. In male mouse liver, GH pulse-activated STAT5 binds to the DNA motif TTCNNNGAA at genomic sites strongly enriched for proximity to male-biased genes, but the mechanistic relationship between STAT5 binding and sex differences in chromatin accessibility and epigenetic marks is unknown. In female liver, persistent activation of STAT5 by the near-continuous presence of circulating GH is associated with a significant enrichment of STAT5 binding nearby female-biased genes [15], however, the underlying mechanisms, including chromatin features that distinguish these STAT5 binding sites from male-biased STAT5 binding sites are poorly understood.

Liver STAT5 is activated by male plasma GH pulses within minutes, enabling STAT5 to rapidly induce the transcription of several male-biased genes [16]; however, a majority of male-biased genes respond slowly to changes in GH secretory patterns [17], likely reflecting the time required for secondary changes, including changes in histone modifications and the underlying chromatin state [18]. Continuous infusion of GH in male mice overrides the endogenous plasma GH pulses and substantially feminizes liver gene expression over a period of days, with female-biased genes already in an active chromatin state in male liver often responding earlier than genes in an inactive chromatin state [17]. Changes in chromatin accessibility occur at genomic regions discovered as DNase-I hypersensitive sites (DHS), a hallmark of epigenetic regulation. These DHS can be assayed by DNase-seq, which identified ∼70,000 DHS in male and female mouse liver [19], including thousands of enhancers, promoters, and insulator regions [20]. More than 4,000 of the 70,000 liver DHS show sex differences in chromatin accessibility, as well as sex-biased binding of STAT5 and other GH-regulated transcription factors. Together, these factors reinforce sex differences in liver gene expression [18, 19] and determine downstream sex differences in disease susceptibility [21–23]. Importantly, male-biased transcription factor binding is strongly enriched at male-biased DHS located nearby male-biased genes, and female-biased transcription factor binding is strongly enriched at female-biased DHS found nearby female-biased genes. Sex differences in chromatin structure and accessibility are thus key features of sex-differential liver gene expression [17, 18]. However, little is understood about the mechanisms linking the sex-dependent temporal GH secretion patterns to the robust sex differences in chromatin accessibility and transcription factor binding that regulate liver gene expression.

Here, we elucidate the relationship between the sex-dependent patterns of plasma GH stimulation of hepatocytes and sex differences in liver chromatin accessibility. We identify more than 800 male-biased enhancer DHS regions, where the direct binding of plasma GH pulse-activated STAT5 induces a dynamic cycle of male liver chromatin opening and closing at sites that comprise 31% of all male-biased DHS. Thus, the pulsatility of plasma GH stimulation per se confers significant male bias in chromatin accessibility and STAT5 binding at a substantial fraction of the genomic sites linked to sex-biased liver gene expression and liver disease.

Furthermore, we establish that a single physiological replacement dose of GH given as a pulse to hypophysectomized (hypox) mice recapitulates, within 30 min, the pulsatile re-opening of chromatin seen in pituitary-intact male mouse liver. We elucidate key epigenetic features distinguishing this dynamic, STAT5-driven mechanism of male-biased chromatin opening from that operative at a second, distinct class comprised of static male-biased DHS, which are constitutively open in male liver but closed in female liver. Finally, our analysis of histone marks enriched at each class of sex-biased DHS indicates that distinct epigenetic mechanisms mediate sex-specific gene repression in each sex. Together, these studies identify pulsatile chromatin opening as a novel mechanism controlling sex differences in chromatin accessibility and transcription factor binding closely linked to sex differences in gene expression and liver disease.

Results

Sex-biased DHS are primarily distal enhancers that target sex-biased genes

Open (accessible) chromatin regions identified in mouse liver by DNase-seq analysis have been classified as enhancers, insulators and promoters based on their chromatin mark patterns and CTCF binding activities (n = 70,211 standard reference DHS set) [20] (Table S2A). Several thousand of these DHS show greater chromatin accessibility in male than female liver (n = 2,729 male-biased DHS) or greater accessibility in female than male liver (n = 1,366 female-biased DHS) [19]. These sex-biased DHS regions were enriched for enhancer marks, with >85% classified as enhancers or weak enhancers in mouse liver, as compared to 66% for all sex-independent DHS (Fig. 1A).

Fig. 1.

Assignment of the sex-biased liver DHS to their putative gene targets (closest RefSeq gene or multi-exonic lncRNA gene transcription start site in the same topologically associating domain (TAD)) [20] revealed that sex-biased DHS were highly enriched for mapping to genes showing the corresponding sex bias in the level transcription, but not for genes whose expression shows the opposite sex bias (Table S2B, Table S3). Cumulative plots of the distance from each DHS to its target gene revealed that a majority of sex-biased DHS classified as enhancers or insulators are distal to their target genes. Sex-biased DHS with H3K4me3 marks, which comprise only 1-2% of all sex-biased DHS, are proximal to their target genes, validating their classification as promoter DHS (Fig. 1B). Thus, a majority of sex-biased DHS have the marks of positively acting distal enhancers, consistent with our finding that sex-dependent DNA looping at an intra-TAD scale is common in mouse liver [24].

No major differences in chromatin state distributions between male-biased, female-biased and sex-independent enhancer DHS were seen based on separate chromatin state maps developed for male and female liver [18]. Thus, all three sets of enhancer DHS showed a high frequency (50-66%) of chromatin state E6, whose emission probabilities (Fig. 1C) indicate a high frequency of DHS in combination with the activating chromatin marks H3K27ac and H3K4me1, and a lower frequency (12-16%) of chromatin state E5 (high frequency of DHS but low frequency of H3K27ac and H3K4me1) (Fig. 1D). Sex-biased and sex-independent insulator DHS also showed similar chromatin state distributions (state E5 ∼ state E6), as did sex-biased and sex-independent promoter DHS, which were primarily in state E7, characterized by a high frequency of H3K4me3 and H3K4me1 marks (Fig. 1D).

The absence of major differences in overall chromatin state distributions between male-biased and female-biased liver DHS prompted us to investigate other factors that may provide insight into underlying mechanisms regulating sex-differences in hepatic chromatin accessibility and their link to sex differences in gene expression.

Impact of endogenous pulses of GH and STAT5 activity in male liver

Sex differences in pituitary GH secretion – pulsatile in males vs. near continuous (persistent) in females – regulate the sex-dependent expression of hundreds of genes in adult mouse liver. This regulation requires the GH-activated transcription factor STAT5 [25, 26]. We hypothesize that the pulsatile activation of STAT5 seen in male liver in direct response to plasma GH stimulation [15, 16, 27] dynamically alters the chromatin accessibility landscape of male mouse liver. More specifically, we propose that the repeated activation of STAT5 in male mouse liver by plasma GH pulses induces dynamic cycles of chromatin opening and closing at a subset of liver DHS, and that this response can be discovered by comparing chromatin accessibility profiles in livers of individual male mice euthanized at a peak vs. at a trough of hepatic STAT5 activity (Fig. 2A, top). To test this hypothesis, we analyzed DNase-seq libraries prepared from liver nuclei purified from 21 individual adult male mice. In parallel, we determined the STAT5 DNA-binding activity of each liver by EMSA (electrophoretic mobility shift analysis) of whole liver extracts.

Fig. 2.

Ten of the 21 livers had high STAT5 EMSA activity (STAT5-high activity livers) and 11 livers had very low or no detectable STAT5 EMSA activity (STAT5-low activity livers) (Fig. 2B, Fig. S1). Next, we performed DNase-seq analysis on genomic DNA fragments released by light DNase-I digestion of nuclei purified from each liver, followed by diffReps analysis [28] comparing the DNase-seq-released DNA fragments from each group of livers. We thus discovered genomic regions showing significant differential chromatin accessibility between STAT5-high and STAT5-low male livers. Principal component analysis using either the top 200 or the top 600 most significant diffReps-identified differentially accessible regions revealed that 18 of the 21 livers gave patterns of DNase-released fragments that correlate with STAT5 activity. Two of the STAT5-high livers and one of the STAT5-low livers were outliers (Fig. 2C). These outliers were also identified by their discordant DNase-seq read count distributions (Fig. 2D) and were excluded from all downstream analyses.

We re-analyzed the remaining 18 STAT5-high and STAT5-low livers using diffReps then filtered the output list to retain regions identified as DHS peak regions by MACS2 analysis of the same 18 DNase-seq datasets. We thus identified n=2,832 genomic sites where chromatin opening is associated with high liver STAT5 activity and n=123 other sites where chromatin opening is associated with low liver STAT5 activity (Table S4A). As each liver represents a time point of peak (STAT5-high livers) or trough (STAT5-low livers) levels of GH pulse-activated liver STAT5 activity (Fig. 2A), the STAT5-high/STAT5-low differential sites identify genomic regions where chromatin dynamically opens or closes in male mouse liver in close association with GH pulse activation of STAT5. The greater chromatin opening in STAT5-high livers as compared to STAT5-low male livers (and as compared to female livers) was visualized in aggregate plots of normalized DNase-seq cuts across the 2,832 genomic regions (Fig. 2E, left). Chromatin accessibility was greater in the STAT5-low livers at the 123 STAT5-low diffReps-identified genomic sites, where STAT5 activation is associated with chromatin closing (Fig. 2E, right). This pattern, where individual male mouse livers largely show either high or low DNase-seq read count distributions at the top differential genomic sites, was also seen in an independent set of nine male liver DNase-seq samples from a second mouse strain generated by the ENCODE consortium: 4 of the 9 livers showed read distributions very similar to the STAT5-high livers, while 5 of the 9 livers were similar to the STAT5-low livers (Fig. 2D, green bars).

Dynamic vs static liver DHS

Comparison of the set of 2,832 STAT5-high genomic sites with our standard reference set of 70,211 liver DHS revealed that 31% (n=834) of the 2,729 male-biased liver DHS described above are STAT5-high sites, i.e., they respond dynamically to endogenous male plasma GH pulses of liver STAT5, whereas only 0.5% of female-biased liver DHS and 2% of sex-independent liver DHS showed this dynamic STAT5 response (Fig. 2F, green). This strong enrichment of GH/STAT5 pulse-induced chromatin opening at male-biased DHS as compared to sex-independent DHS (ES = 12.9; p < 1E-05) identifies intermittent chromatin opening induced by male plasma GH pulses as a mechanism that can explain the male bias in chromatin accessibility for a significant subset (31%) of male-biased DHS.

To determine whether STAT5 binding per se is an important feature of the observed pulsatile changes in chromatin accessibility, we examined normalized DNase-I cut site aggregate plots for the subset comprised of n=1,307 male-biased DHS that bind STAT5 in ChIP-seq analyses. STAT5-high livers showed the highest level of chromatin opening, followed by up to a 2.8-fold lower level of chromatin opening at those same genomic regions in STAT5-low male livers and in female livers (Fig. 3A, plot 1). In contrast, male-biased DHS that did not bind STAT5 (n=1,422 DHS) showed nearly equal chromatin accessibility in STAT5-high vs. STAT5-low male livers but ∼2-fold lower accessibility in female livers (Fig. 3A, plot 2). Thus, the male bias in accessibility at the STAT5-bound but not at the non-STAT5 bound male-biased DHS can be explained by STAT5-induced pulsatile chromatin opening in male liver. It is also apparent that the male bias of the non-STAT5 bound DHS set is associated with chromatin closing in female liver (Fig. 3A, plot 2, black). The conclusion that chromatin is relatively closed at these 1,422 DHS in female liver is also supported by comparing their mean normalized DNase cutting frequency to that of the full genome-wide set of 53,404 STAT5-unbound, sex-independent DHS (Fig. 3A, plot 2, black versus plot 4). Chromatin accessibility was much higher at the subset of sex-independent DHS that bound STAT5 (n=12,712 DHS, plot 3). Importantly, the high chromatin accessibility at these sex-independent DHS was seen in livers of both sexes and was largely invariant between STAT5-high and STAT5-low male livers. We conclude that STAT5 binding is associated with increased chromatin opening, and that dynamic chromatin opening and closing associated with STAT5 binding is a defining characteristic of a distinct subset of male-biased DHS but occurs infrequently at female-biased and sex-independent DHS.

Fig. 3.

The above findings allow us to identify two distinct subsets of male-biased DHS, which differ in their responsiveness to liver STAT5 activation: 1) dynamic male-biased DHS are characterized by a male bias in chromatin accessibility linked to an increase in chromatin opening following the binding of STAT5 when activated by a plasma GH pulse in male liver every 3-4 h [1, 16] (n=834 dynamic male-biased DHS, 31%); and 2) static male-biased DHS are open in male liver constitutively, i.e., throughout the pulsatile, on/off cycles of GH-induced STAT5 activity, and are comparatively closed in female liver (n=1,895 static male-biased DHS, 69%). Supporting this conclusion, mean chromatin opening at the set of 834 dynamic male-biased DHS was 3.8-fold higher in STAT5-high male livers than in STAT5-low male livers or in female livers (Fig. 3A, plot 5). In contrast, the higher chromatin accessibility in male than female liver at the set of 1,895 static male-biased DHS was largely independent of the male liver’s STAT5 activity status. Of note, the mean level of male liver chromatin opening at those 1,895 sites was 2.6-fold lower than the peak level of accessibility seen at the dynamic male-biased DHS regions in STAT5-high male livers (Fig. 3A, plot 6 vs. plot 5) (Table 1A).

Table 1:

The mean level of chromatin accessibility at the set of static female-biased DHS (1,359 sites, Fig. 2F) was 2.7 to 2.9-fold higher in female liver than in male liver, where accessibility was independent of male plasma GH/STAT5 pulses (Fig. 3A, plot 7). Overall, the accessibility of the static female-biased DHS in female liver was very similar to that of the genome-wide set of 64,584 static sex-independent DHS (Fig. 3A, plot 8 vs. plot 7; peak normalized DNase-I activity value of 699 vs 683, Table S6).

STAT5 binding is closely associated with dynamic male-biased DHS

The presence of a canonical STAT5 motif, TTCNNNGAA, is a distinguishing feature of dynamic male-biased DHS: 81% of 834 dynamic male-biased DHS contained a STAT5 motif versus only 38% of 1,895 static male-biased DHS (c.f., 25% background motif frequency at 64,584 sex-independent static DHS) (Fig. 3B). In addition, multiple STAT5 motifs are more frequently found at dynamic DHS than at static DHS (Fig. 3C). Consistent with this, STAT5 binding, determined experimentally by ChIP-seq, occurs at a greater fraction of dynamic than static male-biased DHS (85% versus 32%; Fig. S2A), and the level of STAT5 binding (normalized STAT5 ChIP-seq read counts) was significantly higher at the STAT5-bound subsets of the dynamic vs. static DHS sets (Fig. S2B). De novo motif discovery supported these findings, with top scoring motifs matching the STAT5B motif found in dynamic DHS (Fig. S3A, Fig. S3B) but not in static DHS sets (Fig. S3C-S3E). A close association between STAT5 binding and chromatin opening is also indicated by the higher chromatin accessibility at the DHS subsets associated with STAT5 binding (Fig. S2C, top row vs. bottom row).

These findings establish a close link between STAT5 binding, which is pulsatile in male liver [15], and the repeated opening and closing of chromatin at dynamic DHS. Nevertheless, pulsatile chromatin opening was also found at 124 male-biased DHS (15% of all dynamic male-biased DHS) that did not bind STAT5 (Fig. S2C, plot 1B). Moreover, STAT5 binding alone is not sufficient to insure dynamic, male-biased chromatin opening and closing, insofar as 90% of all sex-independent DHS that bind STAT5 do not undergo dynamic chromatin opening and closing (Fig. S2C, plot 5A (11,473 sites, 90.3%) vs plot 4A (1,239 sites, 9.7%)). Furthermore, only 54% of male-biased DHS that bind STAT5 are dynamic DHS (710 sites), the other 46% (597 sites) being static male-biased DHS, even though they bind STAT5 (Fig. S2C, plot 1A vs 2A), albeit at a lower level than the dynamic male-biased DHS and similar to the STAT5-bound static sex-independent DHS (Fig. S2B). Thus, factors other than pulsatile STAT5 binding per se are required for dynamic chromatin opening and closing to occur. Finally, dynamic DHS showed a strong preference for male-biased chromatin accessibility, with the frequency of dynamic chromatin opening and closing being 5.6-fold higher at male-biased DHS with STAT5 bound (54%) than at the sex-independent DHS with STAT5 bound (9.7%).

Enrichment of other GH-regulated liver transcription factors at dynamic vs static DHS

We used published ChIP-seq data to investigate whether dynamic and static male-biased DHS differ with respect to the binding of other, sex-biased transcription factors (Table 1B, Table S7A). We examined two GH/STAT5-regulated transcriptional repressors that reinforce STAT5 regulation of liver sex differences. One factor, BCL6, is a male-biased protein that can compete for STAT5 binding to chromatin and preferentially represses expression of female-biased genes in male liver [15, 18, 22, 29]. The second factor, CUX2, is a female-specific repressor protein whose binding sites in female liver are enriched nearby male-biased genes, which enables CUX2 to repress those genes in female liver [30, 31]. Both repressors showed significant enrichment for binding to genomic regions defined by male-biased DHS as compared to a background set of static sex-independent DHS. Thus, male-biased BCL6 binding was significantly enriched at the set of dynamic male-biased DHS (ES=2.0, p=E-14) but showed minimal enrichment at static male-biased DHS (ES=1.3, p=1.6E-04). BCL6 binding was also enriched at dynamic sex-independent DHS (Table S7B). In contrast, CUX2 was most highly enriched for binding in female liver at genomic sites that correspond to static male-biased DHS (ES=4.6, p=E-40), and de novo motif discovery identified a CUX2 motif as the top enriched motif in this, but not the other DHS sets (Fig. S3C). This binding of CUX2 may contribute to the greater closure of those DHS in female as compared to male liver (c.f., Fig. 3A, plot 6). Finally, male-biased STAT5 binding sites showed 2.5-fold greater enrichment at dynamic than at static male-biased DHS (ES=58, p=E-301 vs ES=24, p=E-190) (Table S7B), consistent with our findings, above, implicating pulsatile STAT5 binding in dynamic chromatin opening at those sites.

We reanalyzed published mouse liver ChIP-seq data for FOXA1 and FOXA2 [32], pioneer factors implicated in chromatin opening [33, 34], to identify sex-dependent binding sites for each factor. Male-biased binding sites discovered for each FOXA factor showed strong enrichment for binding at male-biased DHS, with ∼2-fold higher enrichment seen at the static male-biased DHS (for FOXA1: ES=6.1 (static DHS) vs. ES=3.7 (dynamic DHS); for FOXA2: ES=44 (static DHS) vs. ES=21 (dynamic DHS) (Table 1B, Table S7B, Table S8). Consistent with this, de novo motif discovery identified a Fox family factor, FOXI1, as a close match for one of the top enriched motifs in the set of static but not in the set of dynamic male-biased DHS (Fig. S3C vs Fig. S3A, p-value 9.2E-06 vs. 2.0E-03). Finally, female-biased FOXA2 binding sites, but not female-biased FOXA1 binding sites, showed strong enrichment for static female-biased DHS (ES=18, p=E-75; Table S7B). Taken together, these findings support the proposal that FOXA2, and to a lesser extent FOXA1, contribute to sex-dependent chromatin opening, in particular at static DHS.

Impact of hypophysectomy on liver chromatin accessibility

Surgical removal of the pituitary gland (hypophysectomy) ablates pituitary GH secretion and thereby abolishes liver STAT5 activation [16], leading to widespread loss of sex-specific liver gene expression [6]. We hypothesized that by ablating GH-induced STAT5 activation and DNA binding, hypophysectomy will lead to closure of many of the open chromatin regions that regulate sex-specific gene expression. Furthermore, we proposed that restoration of a pulsatile GH signal, in the form of a single exogenous GH pulse given by i.p. injection, will reopen chromatin at many of those sites (Fig. 4A). To test this hypothesis, we used DNase-seq to compare the chromatin accessibility profiles of liver nuclei from hypox male and hypox female mice to those of pituitary-intact control liver nuclei. Hypophysectomy induced changes in chromatin accessibility at several thousand sites, including large numbers of sex-independent DHS, with many more genomic regions undergoing chromatin closing than chromatin opening in male liver, but not in female liver (Fig. 4B, Table S4). Importantly, male-biased DHS that responded to hypophysectomy were almost exclusively closed in male liver following hypophysectomy (Fig. 4C, first two rows, % opening vs % closing). A much smaller percentage of static female-biased DHS responded to hypophysectomy, and the responses observed generally involved chromatin opening in male liver and chromatin closing in female liver (Fig. 4C, Table S4H). Thus, sex-biased DHS are subject to both positive and negative pituitary hormone regulation and respond differently between the sexes.

Fig. 4.

DHS responses to pituitary ablation link to class I and class II sex-biased gene regulatory responses

Hypophysectomy abolishes GH-regulated liver sex differences via two distinct mechanisms, which identify two classes of sex-biased genes [6]: hypophysectomy leads to down regulation of class I sex-biased genes in the sex where the gene is more highly expressed, and it leads to up regulation of class II sex-biased genes in the sex where the gene shows lower expression (Fig. S4A, model). Building on this classification, we hypothesized that DHS that close following hypophysectomy are enhancers mapping to sex-biased genes subject to either positive regulation by GH (class I sex-biased genes) or negative regulation by GH (class II sex-biased genes). Supporting this proposal, the DHS that close in male liver following hypophysectomy showed strong, specific enrichment (ES=5.7, p=1.3E-69) for mapping to class I male-specific genes, which are down regulated in hypox male liver; whereas the DHS that close in hypox female liver showed strong, specific enrichment (ES=7.8, p=6.6E-46) for mapping to class I female-specific genes, which are down regulated in hypox female liver (Fig. 4D). In contrast, the DHS that open in hypox male liver were specifically enriched for mapping to class II female-specific genes (ES=3.6, p=1.8E-15), which are up regulated (de-repressed) in male liver following hypophysectomy; whereas the DHS that open in hypox female liver were specifically enriched albeit only moderately (ES=1.9, p=4.8E-04) for mapping to class II male-specific genes, which are up regulated in female liver following hypophysectomy (Fig. S4B). The enrichment patterns exhibited by these four sets of hypophysectomy-responsive DHS mirror the corresponding sex-specific gene response patterns (Fig. 4D, bottom). Thus, all four DHS sets respond to hypophysectomy in a manner consistent with positively-acting regulatory elements linked to sex-biased genes. DHS that are linked to class I sex-biased genes, and whose chromatin closes following hypophysectomy, require pituitary hormone to maintain open chromatin; when pituitary hormones are ablated by hypophysectomy, chromatin closes and their class I sex-biased target genes are repressed. In contrast, DHS that are linked to class II sex-biased genes of the opposite sex bias, and whose chromatin opens following hypophysectomy, are kept in a closed chromatin state by pituitary hormone; when pituitary hormones are ablated, chromatin opens locally at those DHS and their class II sex-biased target genes are de-repressed: expression of female-biased class II genes increases in hypox male liver and expression of male-biased class II genes increases in hypox female liver.

A single exogenous GH pulse rapidly reopens chromatin at dynamic DHS in hypox male liver

Next, we investigated whether GH is the pituitary factor whose loss accounts for the widespread chromatin closing seen in hypox male liver. Fig. 4B shows that exogenous GH treatment induces chromatin opening within 30 min at more than 3,500 DHS, including 71% of the 1,487 dynamic DHS that closed in male liver following hypophysectomy (Fig. 4E). The fraction of reopening chromatin regions increased to 83% when considering the set of male-biased dynamic DHS that closed following hypophysectomy (Table S2A, column AK). Chromatin reopening was sustained for at least 90 min and then decreased substantially back toward baseline by 240 min (Fig. 4B; Fig. 5A, plot 1), at which time liver STAT5 signaling has terminated [16]. Importantly, the mean level of chromatin re-opening induced by an exogenous GH pulse at dynamic male-biased DHS was very similar to that of the same DHS set in STAT5-high male liver (Fig. 5A, plot 1 vs. Fig. 3A, plot 5; peak DNase-I activity value 1886 vs. 1848). Thus, the rapid chromatin opening induced by a single exogenous pulse of GH recapitulates the effects of an endogenous GH pulse in opening dynamic male-biased DHS.

Fig. 5.

Hypophysectomy led to closure of many fewer static male-biased DHS than dynamic male-biased DHS (Fig. 4C). Exogenous GH treatment partially reversed chromatin closing at static male-biased DHS within 30 min, with the effect persisting, even after 240 min (Fig. 5A, plot 2). This persistence is consistent with static male-biased DHS remaining constitutively open in male liver when liver STAT5 activity dissipates between endogenous plasma GH pulses. Static female-biased DHS also showed a decrease in mean chromatin opening following hypophysectomy, with further decreases in accessibility seen following GH pulse treatment (Fig. 5A, plot 3). Finally, in an important control, the mean chromatin accessibility at the set of 64,584 static sex-independent DHS was unchanged following hypophysectomy or GH pulse stimulation (Fig. 5A, plot 4).

Stratification of the full set of 2,729 male-biased DHS by the presence or absence of STAT5 binding, as determined by ChIP-seq, highlighted the STAT5 dependence of the rapid increases in chromatin accessibility stimulated by GH treatment (Fig. 5B, plot 1 vs plot 2). GH stimulated chromatin reopening at the static male-biased DHS was also primarily associated with the STAT5-bound DHS subset (Fig. S5, plot 2A vs 2B). In contrast, GH induced a much smaller increase in chromatin opening at the STAT5-bound subset of sex-independent DHS (Fig. 5B, plot 3 vs plot 1), where chromatin is already open and largely independent of plasma GH pulses (Fig. 3A, plot 3). GH decreased chromatin accessibility within 30 min at the 123 genomic sites with lower chromatin accessibility in STAT5-high liver as compared to STAT5-low liver (c.f., Fig. 2E), followed by a return to baseline by 240 min (Fig. 5C); this response pattern is consistent with the repression of chromatin accessibility at these sites by endogenous GH/STAT5 pulses (Fig. 2E). Finally, sex-independent DHS that do not bind STAT5 were unresponsive to both hypophysectomy and GH pulse treatment (Fig. 5B, plot 4).

A subset comprised of 354 dynamic, STAT5-high responsive DHS that close following hypophysectomy did not reopen even 240 min after GH pulse treatment (Fig. 4E, Fig. S6). These 354 DHS showed significantly lower differential chromatin accessibility between STAT5-high and STAT5-low livers than did the dynamic STAT5-high DHS subset that responded to an exogenous GH pulse (Fig. S6C, set 4 vs set 5). Conceivably, these 354 dynamic DHS may require multiple GH pulses to reopen or may be co-dependent on other pituitary-regulated hormones ablated by hypophysectomy.

Distinct sex-biased histone mark patterns at static and dynamic sex-biased DHS

Our initial analyses revealed no major differences between dynamic and static male-biased DHS regarding the distribution of enhancer vs insulator vs promoter classifications (Fig. S7A) or their overall chromatin state distributions (Fig. S7B). We therefore examined both classes of male-biased DHS for differences in sex-biased histone marks (Table 1C). Male-biased enhancer marks (H3K27ac and H3K4me1) showed strong enrichment at both dynamic and static male-biased DHS when compared to a background set of 64,584 static, sex-independent DHS (Fig. 6A, Table S7B). In contrast, male-biased H3K36me3 marks, which are characteristic of transcribed regions but have also been shown to inhibit the spread of PRC2-catalyzed H3K27me3 repressive marks [35], were enriched at static but not at dynamic male-biased DHS. In addition, static but not dynamic male-biased DHS were significantly depleted of female-biased enhancer marks (H3K27ac and H3K4me1) (Fig. 6B). This result is consistent with the above finding that static male-biased sites are in a comparatively closed state in female liver (Fig. 3A, plot 2; model in Fig. 6C). Static female-biased DHS showed strong enrichments for female-biased enhancer marks (H3K27ac and H3K4me1), for female-biased H3K4me3 promoter marks, and for female-biased H3K36me3 transcribed region marks (Fig. 6B), a pattern that was very similar to the enrichments seen at static male-biased DHS (Fig. 6A).

Fig. 6.

Examination of two repressive histone marks, H3K27me3 and H3K9me3, revealed their use in a unique way in each sex to enforce sex differences in chromatin states at sex-biased DHS (Fig. 6C). Male-biased H3K27me3 marks were specifically associated with static female-biased DHS (Fig. 6A), consistent with our prior work showing deposition of these marks by the Ezh1/Ezh2 enzymatic component of PRC2 as a specific mechanism to repress many female-biased genes in male liver [36]. In contrast, female-biased H3K9me3 repressive marks were significantly enriched at both dynamic and static male-biased DHS (Fig. 6B); however, we did not find any corresponding enrichment of male-biased H3K9me3 marks at female-biased DHS. This novel finding suggests that H3K9me3 marks are specifically used to repress male-biased genes in female liver, contrasting with the specific use of H3K27me3 marks to repress female-biased genes in male liver.

Chromatin state analysis uncovers bivalent-like states at closed sex-biased DHS

We sought to delineate chromatin state differences at sex-biased DHS in each sex by computing the enrichment (or depletion) of each of 14 distinct chromatin states for each DHS set. These 14 chromatin states are defined by a panel of activating and repressive histone-H3 marks and by the presence of open chromatin (Fig. 1C) and were determined separately for male and female mouse liver [18]. Enrichments were calculated using a genome-wide background set comprised of all 64,584 static sex-independent DHS (Fig. 7). In male liver, enhancer state E6 (characterized by the presence of DHS, H3K27ac and H3K4me1; Fig. 1C) was enriched at both dynamic and static male-biased DHS, and at dynamic sex-independent DHS. Promoter states E7 and E8 were significantly depleted, consistent with the paucity of promoter states in the full set of male-biased DHS (Fig. 1A). Dynamic but not static male-biased DHS showed strong depletion in male liver of the inactive chromatin states E1 (H3K27me3 marks) and E2 (major marks not identified) and of state E14 (H3K36me3 marks) (Fig. 7A). Finally, we observed strong enrichment of the inactive state E2 at both dynamic and static male-biased DHS in female liver (Fig. 7B).

Fig. 7.

Dynamic male-biased DHS were enriched for the enhancer state E9 in female liver (ES=2.42, p=1.6E-69; Fig. 7B) but not male liver (Fig. 7A). The emission parameters of state E9 are very similar to those of state E6, namely, it shows a high frequency of the activating chromatin marks H3K27ac and H3K4me1 but lacks the open chromatin (DHS) feature characteristic of state E6 (Fig. 1C). Thus, in female liver, dynamic male-biased DHS contain histone marks that typify an active enhancer but are inactive due to the comparatively closed state of their chromatin. Indeed, the extent of chromatin opening at these DHS in female liver is equivalent to the level of chromatin opening seen at the same set of sites in male liver between GH/STAT5 activity pulses (Fig. 3A, plot 5; model in Fig. 6C). Dynamic male-biased DHS were also enriched for chromatin state E10 in female liver (ES=2.89, p=1.9E-18) but not male liver. State E10 shows a high frequency of H3K27ac marks, but not H3K4me1 marks, and lacks DHS. This pattern of chromatin state E9 and E10 enrichment at dynamic male-biased DHS in female liver, combined with the strong enrichment of H3K9me3 repressive marks (Fig. 6B), indicates that the genomic regions encompassing dynamic male-biased DHS are in a bivalent-like chromatin state in female liver (Fig. 6C). Such a bivalent state could facilitate chromatin opening under certain pathological conditions, e.g., chemical exposure or biological stress.

Static female-biased DHS showed little or no enrichment of specific chromatin states in female liver when compared to the genome-wide background set of static sex-independent DHS (Fig. 7B). This supports a model whereby their female bias in accessibility is largely due to active suppression in male liver (Fig. 6C). Indeed, the genomic regions encompassing static female-biased DHS were strongly enriched in male liver for inactive chromatin states E1 (repressive mark H3K27me3), E2 and E4 (Fig. 7A). Static female-biased DHS showed weak enrichment for enhancer states E10 and E11 in male liver, which are respectively characterized by the active enhancer marks H3K27ac and H3K4me1 but devoid of DHS. This is analogous to the finding, above, that male-biased DHS are enriched for chromatin states replete with active histone marks but deficient in DHS in female liver. Static female-biased DHS also showed modest enrichment in male liver for bivalent state E12, which is characterized by a mixture of activating histone marks and the presence of repressive H3K27me3 marks, consistent with the strong enrichment of the latter marks seen in Fig. 6A, above, and the overall conclusion that at least a subset of static female-biased DHS is in a bivalent-like state in male liver. This pattern is analogous to the bivalent state adopted by dynamic male-biased DHS in female liver, except for the use of distinct marks to effect sex-specific repression in each sex, namely: H3K9me3 in female liver and H3K27me3 in male liver (Fig. 6C). Finally, promoter-like states E7 and E8 were significantly depleted from static female-biased DHS (Fig. 7), consistent with female-biased DHS largely being gene distal sex-biased enhancers.

Distinct gene targets and enriched biological processes of dynamic and static sex-biased DHS

Dynamic and static male-biased DHS both showed strong enrichment for mapping to male-biased gene targets (8-fold and 11.7-fold enrichments, respectively), as did static female-biased DHS for female-biased gene targets (12.5-fold enrichment) when DHS were mapped to the single nearest transcription start site in the same TAD (Table S2B). Further, mapping DHS to putative target genes using GREAT, which typically maps each DHS to two genes, revealed examples where many male-biased genes may be regulated by multiple male-biased DHS. In some cases, all of the associated male-biased DHS were static male-biased DHS (e.g., Cyp4a12a, with 7 static male-biased DHS), while in other cases they were a mixture of dynamic and static male-biased DHS (e.g., Cyp7b1, with 7 dynamic and 8 static male-biased DHS) (Table S9A). Top female-biased genes enriched for nearby static female-biased DHS included Cux2 and three Cyp3a genes, with 10-15 female-biased DHS each (Table S10C). Sex-independent genes that are well-established direct targets of STAT5 [37] include Igf1, a target of 12 dynamic sex-independent DHS and 3 dynamic male-biased DHS, and Socs2, a target of 8 dynamic sex-independent DHS (Table S10D, Table S10A). Other analyses revealed that all three sets of sex-biased DHS were significantly enriched for mapping to hypophysectomy class I-responsive sex-biased genes as compared to hypophysectomy class II-responsive sex-biased genes (Table S2B). This finding indicates that the sex-biased DHS sets are enriched for positively-acting enhancers that close following hypophysectomy, which leads to decreased expression of their class I (i.e., hypophysectomy repressed) sex-biased gene targets.

The gene targets of each sex-biased DHS set were enriched for distinct but partially overlapping Gene Ontology (GO) Biological Processes, as determined by GREAT analysis. Top enriched GO terms common to both dynamic and static male-biased DHS included lipid metabolic process and steroid metabolic process, consistent with the major role that GH plays in the male-prevalence of fatty liver development and liver metabolic disease [3–5]. Unique enriched terms for dynamic male-biased DHS included cell adhesion, gland development and hepaticobiliary development, whereas gene targets of static male-biased DHS were uniquely enriched for cellular response to glucocorticoid stimulus and EGF receptor signaling, among others (Fig. S8, Table S10, bottom).

Static female-biased DHS gene targets were uniquely enriched for long chain fatty acid metabolic pathway and related terms. Finally, the set of dynamic sex-independent DHS mapped to gene targets enriched for terms related to glucose transport and metabolism, IGF receptor signaling and fibroblast proliferation (Fig. S8). This finding is consistent with the widespread metabolic effects that GH has in both sexes, including complex effects on glucose uptake and glucose oxidation, and on hepatic gluconeogenesis and glycogenolysis [38].

Discussion

Sex differences in chromatin accessibility are a central epigenetic feature that enables regulatory proteins, such as the GH-activated transcription factor STAT5, to bind chromatin and regulate hepatocyte gene transcription in a sex-specific manner. However, the underlying mechanisms controlling sex differences in liver chromatin accessibility, which occur at more than 4,000 distinct sites across the mouse genome, are poorly understood. GH, working through its sex-dependent temporal patterns of secretion by the pituitary gland, is the major hormonal factor controlling STAT5-dependent sex differences in liver gene transcription, but little is known about the potential of GH to directly regulate sex differences in the epigenetic landscape required for sex-dependent transcriptional outputs. To address this question, we elucidated the impact of male plasma GH pulses on global patterns of liver chromatin accessibility by analyzing livers from a population of 18 individual male mice euthanized at either a peak or a trough of plasma GH pulse-stimulated hepatic STAT5 DNA-binding activity.

Our findings establish that the naturally occurring, endogenous plasma GH pulses characteristic of males induce dynamic cycles of chromatin opening and closing at several thousand DNase-I hypersensitive sites (DHS) in male mouse liver and that these events comprise one of two major mechanisms regulating the male bias in liver chromatin accessibility. Analysis of sex-dependent transcription factor binding patterns, histone marks and chromatin states elucidated key features distinguishing this dynamic mechanism of male-biased enhancer activation from that of static, GH pulse-unresponsive male-biased DHS and from female-biased DHS (see model, Fig. 6C). GH thus acts at three distinct steps to regulate sex-dependent hepatocyte gene expression, all three involving the GH-activated transcription factor STAT5 (Fig. 8): GH pulse-induced chromatin opening at dynamic male-biased DHS driven by pulsatile GH activation of STAT5, as discussed further, below; direct transcriptional activation of sex-biased genes by GH-activated STAT5; and GH/STAT5-dependent transcriptional regulation of downstream repressors, such as the female-specific CUX2, which binds to male-biased enhancers in female liver to reinforce sex differences in gene transcription.

Fig. 8.

Endogenous GH/STAT5 pulses activate dynamic male-biased DHS

GH pulse-induced chromatin opening in male mouse liver is shown to be directly controlled by the endogenous male rhythm of plasma GH stimulation of hepatocytes, the major liver cell type contributing to sex-biased liver gene expression [39]. The GH pulse-regulated chromatin regions identified here responded to changes in plasma GH levels in a rapid and dynamic manner, with extensive GH-induced chromatin opening occurring within 30 min, as seen when hypophysectomized male mice were given a physiological replacement dose of GH pulse by intraperitoneal injection (Fig. 4A). This time course is consistent with the rapid activation of GH receptor signaling to STAT5, which occurs within 5 min in cell culture [40] and within 15 min in vivo in a hypophysectomized rat model [41]. Importantly, 85% of the dynamic male-biased DHS identified here (710 of 834 dynamic male-biased DHS) were bound by STAT5, which appears to be a key driver of these GH pulse-induced chromatin opening events.

Chromatin closing followed the termination of nuclear STAT5 signaling, which is complete within 4 h [16] and occurs in sufficient time to reset the GH receptor-JAK2 signaling complex and resensitize hepatocytes before the next plasma GH pulse [40] (Fig. 2A). STAT5 binding and chromatin opening at dynamic male-biased DHS were both significantly lower in livers of female mice, where plasma GH levels and liver STAT5 activity are persistent (near-continuous) yet ineffective at maintaining an open chromatin state. Indeed, the extent of chromatin opening at these male-biased DHS in female liver is very similar to the level in livers from male mice euthanized when liver STAT5 activity is low, a time inferred to be between plasma GH pulses (STAT5-low male livers; Fig. 3A, panel 5). Pulsatile chromatin opening stimulated by endogenous plasma GH pulses is thus a unique mechanism for establishing and maintaining male-biased chromatin accessibility and transcription factor binding at 834 male-biased DHS, corresponding to 31% of all male-biased DHS in the liver.

GH pulse-responsive DHS were discovered by comparing global chromatin accessibility patterns in a set of 8 livers collected from mice euthanized at a peak of GH-activated STAT5 DNA-binding activity (STAT5-high livers) to those of 10 other livers from mice euthanized between pulses of STAT5 DNA-binding activity, i.e., when liver STAT5 activity is very low or undetectable (STAT5-low livers). Three other male livers gave DNase-seq profiles inconsistent with their EMSA-determined STAT5 DNA-binding activity. These outlier livers may have come from male mice euthanized just after the onset of a STAT5 pulse [27], when more time is needed to induce chromatin opening (2 outliers), or shortly after STAT5 is deactivated by tyrosine dephosphorylation [42] but prior to the reversal of chromatin opening (1 outlier). Importantly, the top 200 genomic regions showing differential chromatin opening between STAT5-high and STAT5-low livers also separated a group of nine C57Bl/6 male mouse liver DNase-seq samples from the ENCODE consortium into two distinct classes, corresponding to the accessibility patterns of STAT5-high activity livers (n=4) and STAT5-low activity livers (n=5), respectively. Thus, the overall GH-induced signaling pathways and downstream epigenetic events leading to dynamic chromatin opening and closing at these specific genomic regions is robust across studies and mouse strains.

STAT5-high male livers showed the highest levels of chromatin accessibility among the male-biased DHS that bind STAT5. Similarly, a higher level of chromatin opening was seen at sex-independent DHS that bind STAT5 compared to those that do not bind STAT5 (Fig. 3A, panel 3 vs 4). There is thus a close linkage between STAT5 binding and the extent of chromatin opening. This finding is consistent with the proposal that STAT5 acts as a pioneer factor to enable chromatin opening, as was reported for STAT5 action at the Il9 gene locus in Th9 T-cells [43]. It should be noted, however, that pulsatile chromatin opening also occurred at the small subset (15%) of dynamic male-biased DHS that do not bind STAT5, suggesting other GH receptor signaling pathways [44] may play a role in male-biased chromatin opening at those sites. Furthermore, pulsatile STAT5 activation and pulsatile STAT5 DNA binding alone are not sufficient to insure pulsatile chromatin opening, insofar as many sex-independent DHS do bind STAT5, yet are constitutively open in both male and female liver (Fig. S2C, plot 5A). Finally, plasma GH pulse stimulation decreased chromatin accessibility at a small number of liver DHS, most of which were sex-independent. The mechanism for this GH-stimulated decrease in accessibility and its physiological significance are unknown.

Sex-dependent regulation of static male-biased DHS

We identified a second, but less well-defined mechanism that controls the male bias in chromatin accessibility at a distinct DHS set, comprised of 1,895 static male-biased DHS. These DHS represent 69% of all male-biased DHS and mapped to a set of target genes and enriched biological processes distinct from but overlapping with those of the dynamic male-biased DHS. Static male-biased DHS are constitutively open in male liver, with chromatin accessibility largely unchanged across the peaks and valleys of GH-induced liver STAT5 activity. Nevertheless, GH regulates the sex bias of these static DHS, as evidenced by their extensive closure in livers of male mice given a continuous infusion of GH (Table S2E), which mimics the female plasma GH pattern and substantially feminizes liver gene expression [17, 45].

Thus, in contrast to dynamic male-biased DHS, the sex biased chromatin accessibility of static male-biased DHS is primarily due to their closed chromatin state in the persistent GH signaling environment of female liver. Static male-biased DHS are relatively deficient in STAT5 binding when compared to dynamic male-biased DHS, but showed significant enrichment for binding the female-specific repressor protein CUX2 [30], which we infer contributes to their closure in both female liver and in continuous GH-infused male liver (where CUX2 is induced to female-like levels [17, 46]) through its established transcriptional repressor activity [47]. Finally, static male-biased DHS showed strong enrichment of male-biased binding sites for the pioneer factors FOXA1 and FOXA2, which may help maintain their constitutively open chromatin state in male liver.

Novel insights into underlying mechanisms from chromatin state analysis

Both dynamic and static male-biased DHS are largely (∼95%) promoter-distal enhancer DHS, as indicated by their flanking histone modifications. This supports our prior finding that sex-dependent intra-TAD DNA looping mechanisms are common in mouse liver [24] and indicates such looping mechanisms likely bring both sets of male-biased DHS into closer proximity with their sex-biased promoter targets. Dynamic male-biased DHS were enriched for the active enhancer state E6 in male liver but were enriched for enhancer states E9 and E10 in female liver.

Enhancer states E9 and E10 are characterized by a high frequency of same activating chromatin marks as chromatin state E6, namely H3K27ac and H3K4me1 (E9) or H3K27ac alone (E10), but unlike E6 they are both deficient in open chromatin (DHS) (Fig. 1C). Thus, in female liver, dynamic male-biased DHS regions contain active histone marks but are in a comparatively closed chromatin state. This may in part be due to their enrichment for female-biased H3K9me3 (repressive) histone marks, which gives them bivalent character (Fig. 6C). A bivalent chromatin state may protect dynamic male-biased DHS from irreversible silencing [48] and give them the potential for chromatin opening and activation of their latent enhancer activity in female liver, e.g., in response to chemical exposure or biological stress [49]. In contrast, static male-biased DHS are depleted of female-biased enhancer marks (H3K27ac and H3K4me1) and are in an inactive chromatin state in female liver. Notably, the bivalent female chromatin state of dynamic male-biased DHS is distinct from the poised state these genomic regions apparently adopt in male liver between plasma GH pulses (i.e., in STAT5-low male liver) (Fig. 6C) and in hypophysectomized male liver, where a single GH pulse is all that is required to induce rapid chromatin opening, even after several weeks of pituitary hormone ablation.

The set of female-biased DHS identified here is almost entirely (99%) static, i.e., GH pulse-unresponsive, and showed strong enrichment for three female-biased activating chromatin marks, H3K27ac, H3K4me1 and H3K4me3. These DHS also showed strong enrichment for female-biased FOXA2 binding, which likely contributes to their high chromatin accessibility in female compared to male liver (Fig. 6C). The static female-biased DHS also showed strong enrichment for the repressive histone mark H3K27me3 in male liver, which is expected to contribute directly to their closed, inactive chromatin state. Moderate enrichments for several chromatin states characterized by the presence of activating chromatin marks but lacking in DHS (states E10, E11 and E12) was also observed. This indicates heterogeneity within the set of static female-biased DHS, some of which have bivalent chromatin character in male liver. Evidence for heterogeneity of these female-biased DHS also comes from their varied time-dependent loss of H3K27me3 marks and from their differential time courses for chromatin opening when male mice are given GH as a continuous infusion [17]. Heterogeneity was also seen in the extent of female-biased target gene de-repression when H3K27me3 marks are lost from livers of male mice deficient in the H3K27-trimethylase enzyme complex PRC2 [36]. Further study is needed to elucidate the mechanistic basis for these time-dependent epigenetic responses and their associated sex-biased gene expression changes [17].

Novel insight into the fundamental underlying epigenetic mechanisms of sex-biased gene regulation comes from our discovery that distinct repressive histone marks, H3K27me3 and H3K9me3, are used in a unique way in each sex to enforce sex differences in chromatin states at sex-biased DHS. Whereas male-biased H3K27me3 repressive marks are highly enriched at, and are specifically associated with, static female-biased DHS, in agreement with our prior work [18, 36], we show here that female-biased H3K9me3 repressive marks are specifically enriched at both dynamic and static male-biased DHS, without a corresponding enrichment of male-biased H3K9me3 marks at static female-biased DHS. This unexpected finding supports the proposal that sex-specific gene repression is mediated by two distinct mechanisms, with H3K9me3 marks specifically used to repress certain male-biased regulatory sites in female liver, and H3K27me3 marks specifically used to repress female-biased regulatory sites in male liver.

Finally, we discovered that sex-biased H3K36me3 marks are a unique distinguishing feature of static sex-biased DHS, with male-biased H3K36me3 marks being highly enriched at static male-biased DHS but not at dynamic male-biased DHS, and female-biased H3K36me3 marks highly enriched at static female-biased DHS (Fig. 6C). H3K36me3 marks are classically associated with the demarcation of actively transcribed genes [50] but are also used to maintain cell type identity by inhibiting the spread of H3K27me3 repressive marks at cell type-specific enhancers [35, 51]. The enrichment of H3K36me3 marks at static male-biased DHS described here could thus be an important mechanism to maintain sex-dependent hepatocyte identity by keeping static male-biased enhancers constitutively open and free of H3K27me3 repressive marks in male liver, and similarly for H3K36me3 marks enriched at static female-biased DHS in female liver. Further study is needed to elucidate the underlying mechanisms whereby these and the other sex-specific histone marks discussed above are deposited on chromatin in a sex-dependent and site-specific manner and the roles that GH plays in regulating these epigenetic events.

Discovery of regulatory elements associated with class I and class II sex-biased genes

DNase-seq profiling identified more than 5,000 liver DHS that open or close following hypophysectomy. Strikingly, a single physiological replacement dose of GH restored, within 30 min, liver STAT5 activity and chromatin accessibility at 83% of dynamic male-biased DHS that closed following hypophysectomy. These chromatin sites are maintained in a primed (poised) state in hypox male liver over a period of several weeks, despite the prolonged deficiency of GH and other pituitary hormones, including gonadal steroids, and the dysregulation of several thousand liver-expressed genes [6, 16]. DHS that close in male liver following hypophysectomy were enriched for proximity to class I male-specific genes, which are down regulated in hypox male liver, whereas DHS that close in hypox female liver showed the strongest enrichment for class I female-specific gene targets, which are down regulated in hypox female liver (Fig. S4, model). Furthermore, DHS that open in male liver following hypophysectomy were enriched for proximity to class II female-specific genes and DHS that open in female liver were enriched for proximity to class II male-specific genes (Fig. 4D). These enrichments are consistent with the definition of class II sex-biased genes, i.e., genes that are up regulated by hypophysectomy in the sex where they show lower expression in intact mice. Together, these findings lend strong support for the functional role of each class of sex-biased and pituitary hormone-dependent DHS as regulatory elements with intrinsic, positively acting enhancer potential, keeping gene expression on in intact liver by maintaining open chromatin, e.g., at male-biased DHS repressed in hypox male liver in the case of class I male-biased genes, or keeping gene expression off by maintaining a closed chromatin state, e.g., at female-biased DHS induced in hypox male liver in the case of class II female-biased genes. Further study will be required to elucidate the underlying regulatory factors and molecular mechanisms governing these gene regulatory circuits, in particular those controlling the de-repression of class II sex-biased genes following pituitary hormone ablation.

Pituitary GH secretory patterns vs. gonadal steroids as regulators of sex-biased liver chromatin accessibility and gene expression

While testosterone has a well-established role in programming hypothalamic control of pituitary GH secretory patterns [9–11], it is also possible that androgens and estrogens could regulate sex differences in hepatocytes directly at the epigenetic or transcriptional level. However, our findings support the proposal that plasma GH patterns, and not gonadal steroids, dominate epigenetic control of liver sex differences. First, the ability of a single exogenous plasma GH pulse to rapidly reopen dynamic male-biased DHS closed by hypophysectomy – in the face of ongoing ablation of pituitary stimulated gonadal steroid production and secretion – implicates GH signaling per se in the direct regulation of chromatin accessibility for this class of male-biased DHS. Second, GH regulates the sex bias of static male-biased DHS as well, as evidenced by their widespread closure in male liver following continuous GH infusion (Table S2E). It is important to note, however, that hepatocyte-specific knockout of androgen receptor (AR) does, in fact, dysregulate ∼15% of sex-biased genes, albeit with a much lower effect size than global AR knockout [52] due to the systemic disruption of the somatotropic axis and circulating GH secretory profiles [53, 54]. Conceivably, AR could regulate these genes by a direct binding mechanism, acting either alone or in concert with GH-activated STAT5 to keep chromatin open constitutively at a subset of static male-biased DHS, of which 32% undergo at least partial closure in male liver following hypophysectomy (Fig. 4C). Estrogen receptor (ERα) likely plays only a minor role in regulating sex-biased liver DHS enhancers, given the lack of effect of hepatocyte-specific ERα knockout on sex-biased liver gene expression [22] and our finding that only 12% of static female-biased DHS close in female liver following hypophysectomy, which decreases circulating estradiol levels [55].

Methods

Animal treatments and DNase-seq analysis

All mouse work was carried out in accordance with ARRIVE Essential guidelines 2.0 [56] for study design, sample size, randomization, experimental animals and procedures, and statistical methods, and with approval of the Boston University Institutional Animal Care and Use Committee. Male and female CD1 mice (ICR strain, strain code # 022), 7-weeks old and purchased from Charles River Laboratories (Wilmington, MA), were kept on a 12-hour light/dark cycle with food and water without restriction. Livers were collected from individual untreated mice between 8 and 10 weeks of age. Where indicated, male and female mice were hypophysectomized (hypox) by the supplier at ∼7-8 weeks of age.

Completeness of hypophysectomy was verified by the absence of weight gain over a 2-3 week period and by the lack of Mup protein in urine samples (SDS-PAGE analysis) [16]. Hypox male mice were treated with recombinant rat GH (purchased from Prof. Arieh Gertler, Protein Laboratories Rehovot Ltd, Rehovot, Israel), given as a single intraperitoneal injection at 125 ng of GH per gram body weight [6], or vehicle (control), and were euthanized by cervical dislocation 30, 90 or 240 min later [16]. This dose is 12 times lower than the supraphysiological dose of GH widely used to study reactivation of hepatic Igf1 following hypophysectomy {Alzhanov, 2015 #77}. A protein extract prepared from a small piece of each liver to be used for DNase-seq analysis (see below) was assayed for liver STAT5 DNA-binding activity by EMSA, as described [16]. Results of this assay allowed us to classify each individual male mouse as STAT5-high activity (n=10) or STAT5-low activity (n=11) at the time of euthanasia and liver collection. EMSA was also used to verify the ablation of liver STAT5 activity following hypophysectomy and the effectiveness of exogenous GH administration at restoring endogenous STAT5-high liver levels of STAT5 EMSA activity within 30 min.

DNase-seq libraries, Illumina sequencing, and data processing

Livers used for these analyses were obtained from intact male and female mice, from hypox male and hypox female mice, and from hypox male mice given a single injection of GH and euthanized 30, 90 or 240 min later (n=6-12 livers per group). Nuclei were purified from individual fresh livers as described [19] and stored at −80 °C. Nuclei were subsequently treated with DNase I to release DNA fragments from hypersensitive genomic regions to identify liver DNase hypersensitive sites (DHS) [16]. The DNA fragments released from each liver were combined (3 livers per pool) to give n=2 to n=4 independent pooled samples, which were used to prepare DNase-seq libraries for each mouse group (Table S1). In the case of STAT5-high and STAT5-low liver DHS analysis, however, DNase-seq libraries were prepared directly from DNase-I digested fragments released from each individual liver without pooling, as described below. DNase-seq libraries were prepared using the NEBNext Ultra Library Prep Kit for Illumina (New England Biolabs). Illumina sequencing was performed on an Illumina HiSeq instrument to a depth of 12 to 33 million mapped 40-50 bp single-end sequence reads per sample. Detailed sequencing statistics are shown in Table S1. Raw and processed data are available at www.ncbi.nlm.nih.gov/gds under accession numbers GSE131848 and GSE131852 (SuperSeries GSE131853).

Sequencing data was analyzed using a custom in-house DNase-seq pipeline [49]. Briefly, the pipeline processes raw FASTQ files and outputs various control metrics, including FASTQC reports, confirmation of read length, verification of the absence of read strand bias, and identification of contaminating adaptor sequences using Trim Galore (https://github.com/FelixKrueger/TrimGalore). Reads were mapped to mouse genome mm9 using Bowtie2 (v2.2.6) [57]. Regions of DNase hypersensitivity (DHS) were discovered as peaks identified by MACS2 (v2.1.0.20150731) [58] using the option (-nomodel –shift −100 –extsize 200) to inhibit read shifting, and the option (-keep-dup) to retain all reads that contribute to the peak signal. Peaks were discovered for each individual DNase-seq library then filtered to remove peaks that overlap ENCODE blacklisted regions [59] as well as peaks comprised of 5 or more identical reads that do not overlap any other read (“straight peaks”). Nine additional individual male mouse liver DNase-seq datasets generated by the ENCODE consortium (FASTQ files downloaded at http://genome.ucsc.edu/cgi-bin/hgFileUi?db=mm9&g=wgEncodeUwDnase; GEO accession: GSM1014195, liver replicates #5 to #13) were analyzed using the same DNase-seq pipeline (cell line, liver; strain, C57BL/6; age, adult 8 weeks; sex, male). These DNase-seq datasets were single end sequencing reads, 36 bp in length, with a sequencing depth of 19 to 37 million total mapped reads per sample.

Liver DHS classification and DHS gene target expression data

A set of 72,862 mouse liver DHS regions previously identified by DNase-seq analysis of male and female mouse liver [19] was filtered to remove 515 DHS regions that overlapped ENCODE blacklisted regions. 2,136 DHS regions that could not be classified according to a 5-class DHS model defined previously [20] were also removed to obtain a standard reference set comprised of 70,211 mouse liver DHS. Each DHS was designated as a promoter, weak promoter, enhancer, weak enhancer, or insulator DHS based on the 5-class model, which is primarily based on ChIP-seq signals for the H3 histone marks H3K4me1and H3K4me3, combined with CTCF ChIP-seq binding in adult male mouse liver [20]. Weak enhancers were identified by their low levels of H3K27ac ChIP-seq signal compared to DHS classified as enhancers, combined with a distance from RefSeq gene transcription start sites inconsistent with a weak promoter designation [20]. Further, individual DHS were designated male-biased (n=2,729), female-biased (n=1,366), or sex-independent (n= 66,116) based on significant differences in chromatin accessibility between male and female mouse liver [19]. Each DHS was assigned a single putative gene target (RefSeq gene or multi-exonic lncRNA gene {Melia, 2019 #21}) corresponding to the closest transcription start site within the same TAD [20], except as noted. Table S2A lists DHS target genes, their overlap with liver STAT5 binding sites determined by ChIP-seq [15], their DNase-seq activity (chromatin accessibility) in male and female mouse liver, and their responses to hypophysectomy and to hypophysectomy + GH treatment. For some analyses, including Gene Ontology (GO) enrichments, GREAT (v.4.0.4) (http://great.stanford.edu/public/html/index.php) [60, 61] was used to identify target genes using default settings for the Basal + extension gene association rule, which typically assigns two genes to each DHS. GREAT output is provided in Table S9 and Table S10.

Differential DHS between STAT5-high and STAT5-low activity livers

Nuclei purified from each individual STAT5-high activity (n=10) and STAT5 low-activity mouse liver (n=11) (Table S1), classified based on EMSA as described above, were analyzed by DNase-seq to discover genomic regions (i.e., DHS) that responded dynamically to endogenous plasma GH pulses and the associated changes in liver STAT5 DNA-binding activity, as follows. diffReps analysis [28] was used to discover genomic sites that were more open or were more closed (DNase-seq normalized intensity |fold-change| > 2 and FDR < 0.05 [Benjamini-Hochberg adjusted p-value]) for the set of EMSA-identified STAT5-high activity livers as compared to the set of EMSA-identified STAT5-low activity livers. The diffReps nucleosome option (200 bp window size) was used and the option (-frag) was set to zero for all comparisons. The differential sites that diffReps identified as significant were further analyzed by two methods, principal component analysis and boxplot analysis, to determine the distributions of normalized sequence read counts (reads per kilobase per million mapped reads) for each individual DNase-seq library. Three outlier DNase-seq samples were thus identified: two from livers designated STAT5-high activity based on EMSA (samples G74A_M1 and G74A_M2), which gave DNase-seq read count distributions (top 200 and top 600 diffReps-identified differential sites, ordered by decreasing fold-change in chromatin opening) more similar to STAT5-low activity livers; and one from a STAT5-low EMSA activity liver (sample G92_M5), which gave a DNase-seq read count distribution more similar to STAT5-high activity livers. These 3 outlier liver DNase-seq libraries were excluded from all downstream analysis. Next, we implemented diffReps analysis using the same cutoffs described above to compare the DNase-seq activity profiles of the remaining EMSA-identified STAT5-high activity (n=8) and STAT5-low activity (n=10) DNase-seq samples. The resultant set of diffReps differential sites was overlapped with a MACS2 DHS peak union list, which was obtained by merging the MACS2 DHS peak calls from all 18 male mouse liver DNase-seq samples using the BEDTools Merge command. This overlap analysis yielded a final set of 70,767 MACS2 DHS peak union sites (Table S4B), of which 2,832 MACS2-identified DHS that were more open (more accessible chromatin state) and 123 DHS were in a more closed chromatin state based on their overlap with the diffReps-identified STAT5-high vs. STAT5-low livers differential regions (Fig. 2E). The other 67,812 DHS regions were designated static DHS, as they did not overlap a diffReps-identified STAT5-high vs. STAT5-low differential region.

Comparison of these 70,767 MACS2-identified DHS peak union sites with our standard reference set of 70,211 mouse liver DHS identified n=2,373 reference set DHS that overlapped a STAT5-high differential site (Fig. 2F), n=98 reference set DHS that overlapped a STAT5-low differential site, and n=45,754 liver DHS that overlapped a static DHS site (Table S2A, columns H and I). The remaining 21,986 standard reference set DHS did not overlap the 70,767 DHS peak union list and were also labeled static DHS. We applied the designation of sex bias determined previously [19], namely, male-biased, female-biased, or sex-independent, to each DHS that overlapped the reference set 70,211 liver DHS (Table S2A, column F). The standard reference set liver DHS were further designated dynamic DHS if they overlapped a STAT5-high > STAT5-low differential DHS (i.e., a GH/STAT5 pulse-opened DHS). Liver DHS not identified as dynamic (including the 98 STAT5-low > STAT5-high DHS peak union sites that overlap a standard reference set DHS) were designated static with respect to chromatin accessibility changes induced by endogenous GH-induced STAT5 pulses. Sex-specific DHS were thus designated dynamic male-biased DHS, static male-biased DHS, dynamic female-biased DHS, or static female-biased DHS with respect to GH/STAT5-induced chromatin opening (Table S2A, column I). Sex-independent DHS were similarly designated dynamic or static sex-independent DHS (Fig. 2F).

DNase-I cut site aggregate plots

Normalized DNase-I cut site aggregate plots were generated using an input DNase-seq dataset and the set of input genomic regions (DHS sequences) to be analyzed by sequence read counting. First, FASTQ files from DNase-seq biological replicates were concatenated to obtain a single combined replicates file for each treatment group. For example, for male liver, we generated a single combined STAT5-high DNase-seq FASTQ file by merging the n=8 biological replicate STAT5-high liver FASTQ files, and separately, we generated a single combined STAT5-low sample by merging the n=10 biological replicate STAT5-low DNase-seq FASTQ files (Table S1). Combined replicate FASTQ files were generated for intact female, hypox female, hypox male, and each of the hypox male + GH treatment time point DNase-seq samples in the same manner. Second, for each combined replicate file, a BED file comprised of positive and negative strand reads was processed to determine the DNase-I cut site, which corresponds to the 5’-end of each sequence read. Third, the set of input genomic regions was processed to generate a list of 2 kb regions centered at the midpoint of each DHS. The BEDTools Coverage command using the (-d) option was used to count the number of DNase-I cut sites at each nt position of the 2 kb midpoint-centered regions, thus producing a read count matrix composed of 2,000 read counts for each input genomic region. Fourth, a custom R script was used to load the read count matrix, calculate the sum of read counts at each nt position, normalize the raw read counts by the number of reads in the subset of 70,211 standard reference DHS to be analyzed, normalize by the number of input genomic regions, and then generate a plot of the DNase-I signal across the input genomic regions. The resultant normalized DNase-I cutting profiles were then smoothed using a LOWESS smoother as implemented in R (package: gplots v3.0.1.1). An offset was applied to the profile by subtracting an average read count (calculated from the first 200 nt positions) from the normalized read count intensity to standardize the baseline of each profile. Profiles for other, related DNase-seq datasets were processed in the same manner and were plotted on the same set of axes to enable direct comparisons across all such plots in Figs. 2, 3 and 5 and Figs. S2, S5 and S6.

Impact of hypophysectomy and GH injection on chromatin opening and closing

MACS2 was used to discover differential DHS peaks in DNase-seq samples prepared from intact male and intact female mouse liver, hypox male and hypox female mouse liver, and from livers of hypox male mice given a single injection of GH and euthanized either 30, 90, or 240 min later (Table S1). The effects of hypophysectomy on chromatin accessibility were determined by comparing hypox liver DNase-seq samples to the corresponding samples prepared from intact liver (control) samples in male liver and separately in female liver. The effects of a single injection of GH were determined by comparing DNase-seq samples from livers of hypox male mice treated with GH to the untreated hypox male liver controls at each time point. For each comparison, a DHS peak union list was generated by merging the MACS2 DHS peak calls from all individual biological replicates for the corresponding treatment group. For example, a single peak union list for the hypox male compared to intact male liver samples was generated by merging all the MACS2 peaks from each of the n=8 individual male mouse liver DNase-seq samples (n=6 intact male control livers and n=2 hypox male livers; Table S1). Genomic regions that showed significantly differential DNase-seq signal between hypox and intact control livers, or between hypox male + GH and hypox male livers were discovered separately for each comparison using diffReps, using the nucleosome option (200 bp window). Significance was based on |fold-change| > 2 and FDR < 0.05 (Benjamini-Hochberg adjusted p-value) for diffReps-normalized signal intensity values. The diffReps-identified differential sites were then filtered by their overlap with the peak union list for all 8 samples to obtain the sets of differential DHS (e.g., 2,142 DHS that open and 4,856 DHS that close in male liver following hypophysectomy; Table S4A). MACS2 peak union DHS that did not overlap a diffReps region were annotated as static DHS peaks (50,055 sites). Concatenation of the diffReps-identified differential DHS with this set of static DHS yielded the full set of 57,053 DHS peak union sites for this dataset. Each of the 70,211 standard reference set liver DHS (see above) was then labeled based on whether it overlapped a diffReps-identified differential DHS that opened following hypophysectomy, a diffReps differential DHS that closed following hypophysectomy, a static DHS, or “none”, for those reference set DHS that did not overlap any of the set of 57,053 hypophysectomy study DHS peak union sites. Corresponding analyses were performed for the intact and hypox female liver samples, and for the hypox male + GH time-course comparisons mentioned above. Table S4 summarizes the numbers of hypophysectomy and hypophysectomy + GH responsive differential DHS for each of these comparisons, and the subsets that overlap the reference set of 70,211 liver DHS.

Sex-biased genes and hypophysectomy response classification

RNA-seq data from intact male and intact female mouse liver was previously used to identify sex-biased genes based on a gene list comprised of 24,197 RefSeq genes [16]. This list was expanded to include both RefSeq and multi-exonic sex-specific lncRNA genes, as follows. Sex-specific RefSeq genes were identified from polyA-selected total liver RNA-seq samples from intact male and intact female liver using the “genebody” method of sequence read counting [16] at a threshold of |fold-change| > 1.5, adjusted p-value < 0.05 and FPKM > 0.25 for the sex with a greater signal intensity, which yielded 387 male-specific and 517 female-specific RefSeq genes. Sex-specific multi-exonic lncRNA genes were identified from polyA-selected liver nuclear RNA-seq samples in male and female mouse liver using the “ExonCollapsed” read counting method [62] at a threshold of |fold-change| > 2, adjusted p-value < 0.05 and FPKM > 0.25 for the sex with a higher signal intensity, which yielded 121 male-specific and 102 female-specific multi-exonic lncRNA genes, for a total of 508 male-specific and 619 female-specific genes (Table S3). A stringent set of sex-independent genes was defined by FPKM > 1 in both male and female liver, |fold-change| < 1.2, and adjusted p-value > 0.1, which yielded a total of 7,253 stringently sex-independent genes, including lncRNA genes (Table S3, column N). Sex-biased genes responsive to hypophysectomy were defined by |fold-change| > 2 and adjusted p-value < 0.05 for hypox mouse liver versus intact mouse liver, determined separately for male liver and for female liver. Using these thresholds, RNA-seq data from intact and hypox male and female liver samples were used to assign sex-biased genes into classes I and II and their corresponding subclasses (IA, IB, IC, IIA, and IIB) [6, 16]. Class I sex-biased genes are those sex-biased genes that are down regulated by hypophysectomy in the sex where they show higher expression in intact mice. Class II sex-biased genes are those that are up regulated by hypophysectomy in the sex where they show lower expression in intact mice (see model in Fig. S4A). Subclasses A, B, and C indicate the response to hypophysectomy in the dominant sex (class II genes) or in the opposite sex (class I genes), as defined in Table 3 of [6]. The number of sex-biased genes in each hypophysectomy response class is shown in Table S3.

STAT5 binding and motif analysis

ChIP-seq analysis of STAT5 binding in male and female mouse liver identified 15,094 merged peaks comprised of male-enriched and female-enriched, and male-female common STAT5 binding sites [15]. A small subset of these STAT5 binding sites, 75 peaks, overlapped ENCODE blacklisted regions and was excluded from downstream analysis. BEDTools overlap analysis of these STAT5 binding sites allowed us to designate each of the 70,211 reference DHS as STAT5-bound if a STAT5 binding site was present, or as ‘not bound’. Table S2A (columns J-N) provides full details, including the STAT5 binding site, the sex-specificity of STAT5 binding (male-enriched, female-enriched, or common to both sexes), the normalized ChIP-seq read counts for STAT5 binding in STAT5-high activity male livers (MH) and in STAT5-high activity female livers (FH), and the average of these two sets of sequence read counts. The STAT5B motif M00459 from the TRANSFAC motif database (Release 2011.1) [63] was used to determine the frequency of STAT5 motif occurrence in each set of DHS sequences. Motifs found in DHS sequences were identified using FIMO (v4.12.0) [64] using the option (--thres 0.0005) to improve detection of short length motifs. The number of STAT5B motif occurrences in each of the 70,211 reference set DHS sequences is shown in Table S2A, column O.

Enrichment analysis

For all enrichment calculations described below, the significance of the enrichment was determined by a Benjamini and Hochberg adjusted Fisher’s Exact test p-value as implemented in R. Enrichments with adjusted p-value < 1E-03 were considered statistically significant. Enrichments of sex-biased DHS for being enhancer, promoter, or insulator regions (e.g., male-biased DHS for being enhancers compared to the background set of sex-independent DHS) were calculated, as follows: Enrichment score = (ratio A) / (ratio B), where: ratio A = the number of sex-biased DHS that are enhancers, divided by the number of sex-biased DHS that are not enhancers; and ratio B = the number of sex-independent DHS that are enhancers, divided by the number of sex-independent DHS that are not enhancers. For example, 2,551 male-biased DHS were classified as enhancers, and 178 other male-biased DHS were not enhancers (2,551/178=14.3), whereas 43,591 sex-independent DHS were classified as enhancers, and 22,525 sex-independent DHS were not enhancers (43,591/22,525=1.93), which gives an enrichment score A/B=14.3/1.93=7.41. The set of 66,116 sex-independent DHS was used as the background for these enrichment calculations.

Enrichments of sex-biased DHS subsets (enhancer, insulator, or promoter) for mapping to sex-biased genes were calculated (e.g. male-biased enhancer DHS mapping to male-specific genes) as follows (Table S2B): Enrichment score = (ratio A) / (ratio B), where: ratio A = the number of male-biased enhancer DHS that map to male-specific genes, divided by the number of male-biased enhancer DHS that map to sex-independent genes; and ratio B = the number sex-independent enhancer DHS that map to male-specific genes, divided by the number of sex-independent enhancer DHS that map to sex-independent genes. For example, among the male-biased DHS classified as enhancers, 404 male-biased enhancer DHS map to male-specific genes, and 525 male-biased enhancer DHS map to sex-independent genes (404/525 = 0.77), whereas 1,495 sex-independent enhancer DHS map to male-specific genes, and 15,457 sex-independent enhancer DHS map to sex-independent genes (1,495/15,457 = 0.097), which gives an enrichment score A/B = 0.77/0.097 = 8.0. The set of 66,116 sex-independent DHS was used as the background for these enrichment calculations.

Enrichments of hypophysectomy-responsive DHS for mapping to class I or class II sex-biased genes (e.g. DHS that close in response to hypophysectomy in male liver and that map to male class I genes) were calculated for male and female liver as follows: Enrichment score = ratio A/ ratio B, where: ratio A = the number of DHS that open (or that close) following hypophysectomy and that map to class I (or to class II) sex-biased genes, divided by the number of DHS that open (or that close) following hypophysectomy and that map to sex-independent genes; and ratio B = the number of hypophysectomy-unresponsive DHS (static DHS) that map to class I or class II sex-biased genes, divided by the number of static DHS that map to sex-independent genes. For example, in male liver, 217 DHS that close following hypophysectomy map to a male class I sex-specific gene, and 1,203 other DHS that close map to sex-independent genes (217/1,203=0.1803), whereas 444 static DHS map to a male class I sex-specific gene, and 14,053 static DHS map to a sex-independent gene (444/14,053 = 0.0316), which gives an enrichment score A/B = 0.1803/0.0316 = 5.7. The static DHS used for these enrichment calculations correspond to the set of 35,562 static DHS in male liver and 30,394 static DHS in female liver.

Enrichments for overlap with genomic regions containing biologically relevant sets of transcription factor binding sites, chromatin marks and combinations of epigenetic features (chromatin states) were calculated for each of the following 4 sets of DHS: (1) dynamic male-biased DHS (834 sites), (2) static male-biased DHS (1,895 sites), (3) dynamic sex-independent DHS (1,532 sites), and (4) static female-biased DHS (1,359 sites) relative to a background set of static sex-independent DHS (64,584 sites) (see Fig. 2F). Briefly, these sets of DHS were identified based on their sex-specificity [19] and were classified as “dynamic” if they overlapped a STAT5-high vs. STAT5-low differential DHS, otherwise they were classified as “static”. Characterization of the 70,211 reference set DHS members as static or dynamic and their sex-specificity designations are listed in Table S2A. Biologically relevant regions annotated include sex-biased transcription factor binding sites for STAT5 and BCL6 [15], CUX2 [30], and FOXA1/FOXA2 [32]. We also examined DHS enrichment for sex-biased chromatin marks and chromatin states defined for male and female liver [18]. FOXA1 and FOXA2 ChIP-seq data [32] was processed and reanalyzed with MAnorm [65] in analyses performed by Gracia Bonilla of this laboratory, which identified sex-biased FOXA1 and sex-biased FOXA2 binding sites in mouse liver. Data from the MAnorm analysis of FOXA1 and FOXA2 ChIP-seq peaks (Table S8) was filtered by log2(fold-change) value (i.e., M-value) to define male-biased and female-biased binding sites. Genomic coordinates for each of the sex-biased transcription factor binding sites, chromatin marks, and chromatin states are provided (Table S5, Table S7).

Enrichments of DHS for overlapping biologically relevant regions were calculated (e.g., dynamic male-biased DHS for containing STAT5-High (Male) binding sites compared to the background set of static sex-independent DHS) as follows: Enrichment score = ratio A/ ratio B, where: ratio A = the percentage of dynamic male-biased DHS with a STAT5 binding site; and ratio B = the percentage of static sex-independent DHS with the same type of binding site. For example, 235 dynamic male-biased DHS contain a STAT5-High (Male) binding site (235/834 = 0.2818), whereas 311 static sex-independent DHS contain a binding site (311/64,584 = 0.0048), which gives an enrichment score A/B = 0.2818/0.0048 = 58.7. The set of 64,584 static sex-independent DHS (Fig. 2F) was used as the background for these enrichment calculations.

Chromatin state map analysis

Chromatin state maps (14 state model), developed for male mouse liver, and separately for female mouse liver, are based on epigenetic data from a panel of six histone marks and DHS data in each sex [18]. These maps were used to identify sex differences in chromatin state and chromatin structure at each DHS region and their relationships to sex-biased gene expression, as follows. BEDTools was used to assign one of the 14 chromatin states to each of the 70,211 reference DHS based on the overlap of the DHS with the male liver chromatin state map, and separately, based on its overlap with the female liver chromatin state map (Table S2A, columns AN-AQ). DHS whose genomic coordinates span two or more different chromatin states were assigned to the state with the largest number of overlapping base pairs, or in case of equal numbers of base pairs, the chromatin state with the smaller genomic coordinates. Enrichments of sets of static and dynamic male-biased DHS for being in one of the 14 chromatin states in male liver (e.g., dynamic male-biased DHS for being in chromatin state E6, which is an enhancer-like state) were calculated as follows: Enrichment score = (ratio A) / (ratio B), where: ratio A = the number of male-biased DHS that are in a particular chromatin state, divided by the number of male-biased DHS not in that chromatin state; and ratio B = the number of static sex-independent DHS in that chromatin state, divided by the number of static sex-independent DHS not in that chromatin state. For example, 604 dynamic male-biased DHS are classified as chromatin state E6, and 230 other dynamic male-biased DHS are not classified as chromatin state E6 (604/230 = 2.626), whereas 24,082 static sex-independent DHS are classified as chromatin state E6, and 40,502 static sex-independent DHS are not classified as chromatin state E6 (24,082/40,502 = 0.5946), which gives an enrichment score A/B = 2.626/0.5946 = 4.4. The set of 64,584 static sex-independent DHS was used as the background set for these enrichment calculations.

DHS peak normalization

DNase-seq data to be visualized in the UCSC Genome browser (https://genome.ucsc.edu/) was normalized using the number of sequence reads in each DHS peak region per million mapped sequence reads (reads-in-peaks-per-million, RiPPM) as a scaling-factor [49]. Normalization was carried out using a comprehensive list of DHS peak regions merged across each dataset (peak union list), obtained by concatenating FASTQ files for biological replicates of each control and treatment group, as described above. DHS peak regions identified in both individual liver samples and by analysis of the combined samples were concatenated into a single list, and then the BEDTools Merge command was used to combine overlapping features to generate a single list of non-overlapping DHS peaks. The fraction of reads in peaks for each sample was then calculated to obtain a scaling factor. Raw sequence read counts were divided by this per-million scaling factor to obtain RiPPM normalized read counts.

Acknowledgements

Supported in part by NIH grant DK121998 (to DJW).

Authors’ contributions

AR and DJW conceived of the study. All animal studies and wet lab analyses including STAT5 EMSA gels, isolation of liver nuclei, DNase-seq analysis and library preparation were performed by JC. Data processing, primary computational analysis, preparation of analysis scripts, and data visualization were carried out by AR. All other data analysis and preparation of datasets and figures for publication were carried out by AR and DJW. AR and DJW jointly prepared a preliminary draft of the manuscript. The final manuscript was written by DJW and was reviewed and approved by all the authors. DJW supervised the overall project and edited the final manuscript for publication.

Supplemental figures

Fig. S1.

Fig. S2.

Fig. S3.

Fig. S4.

Fig. S5.

Fig. S6.

Fig. S7.

Fig. S8.