HNF-1β primarily binds gene promoters and intragenic regions in the developing kidney
To identify HNF-1β binding sites in the developing kidney, we performed chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq). Kidneys (metanephroi) were dissected from wild-type embryos at E14.5 and pooled. Kidneys from two independent pooled litters were analyzed separately (two biological replicates). Chromatin was extracted from the kidneys, cross-linked, and immunoprecipitated with an anti-HNF-1β antibody. Cross-linking was reversed, and DNA libraries were generated and sequenced. Read depth and ChIP-seq quality measures can be found in Supplementary Table 1. Mapping to the mouse genome (mm10) identified 8,284 HNF-1β binding sites that were common in the two biological replicates. Peaks were found at known HNF-1β target genes, including Cdh16 and Pcbd1 (Fig. 1a). The locations of the peaks matched the locations of previously identified HNF-1β binding sites [23, 24], which validated the ChIP-seq methodology in embryonic kidneys. De novo motif analysis was performed on the top 500 HNF-1β peaks after filtering out low-complexity regions and simple repeats (Fig. 1b). The HNF-1β consensus sequence was highly represented in these peaks (p = 2.9x10− 13).
Of the HNF-1β binding sites, 26% were located in gene promoters, which were defined as 3 kb upstream from the transcription start site (TSS); 37% were intragenic, defined as the 5′-UTR, 3′-UTR, introns, exons, and 300 bp downstream from the transcription end site (TES); and 37% were located in distal intergenic regions, defined as > 300 bp downstream of the TES or > 3 kb upstream from the TSS (Fig. 1c). HNF-1β binding was visualized across the gene body by aligning the TSS and TES of each protein-coding gene and creating a heatmap of HNF-1β ChIP-seq reads. HNF-1β showed the highest binding in promoters immediately upstream from the TSS, and lower levels of binding were observed within the body of genes. Another small peak of enrichment was found immediately downstream from the TES (Fig. 1d).
Of the total HNF-1β binding sites, 5,068 sites were located within or in closest proximity to protein-coding genes. Most of these sites were located in gene promoters (30%) or were intragenic (51%). Only 19% were distal intergenic. A total of 982 sites mapped to long non-coding RNAs, of which 17% were in the promoter region, 27% were intragenic, and 56% were distal intergenic. A total of 134 sites mapped to microRNAs, of which 29% were in the promoter region and 71% were distal intergenic. The remaining gene classifications of HNF-1β peaks can be found in Supplementary Table 2. Taken together, the HNF-1β consensus sequence highly represented in the E14.5 ChIP-seq, as well as known HNF-1β binding sites. The majority of HNF-1β peaks in protein-coding genes were within the promoter and gene body, whereas most peaks in non-coding genes were located in distal intergenic regions.
HNF-1β binding sites in the developing kidney colocalize with open chromatin and histone modifications
To determine if the HNF-1β binding sites were functional, we examined the relationship between HNF-1β binding and chromatin accessibility. Transcriptional activation is associated with open chromatin, which can be measured using ATAC-seq (Assay for Transposase Accessible Chromatin) [25]. Moreover, HNF-1β itself functions as a pioneer transcription factor that remains bound to DNA through mitosis and can activate transcription by promoting histone modifications that result in increased chromatin accessibility [26]. To determine whether HNF-1β binding in embryonic mouse kidneys is associated with increased chromatin accessibility, we compared the locations of HNF-1β ChIP-seq peaks and ATAC-seq peaks. Here, we took advantage of an existing ATAC-seq dataset from the ENCODE database [27, 28] obtained from wild-type E14.5 mouse kidney [29], the same tissue source used for HNF-1β ChIP-seq. We identified 2,263 HNF-1β ChIP-seq peaks that colocalized with ATAC-seq peaks, which represented 27% of the total HNF-1β peaks. To validate this colocalization with HNF-1β activity, binding sites at known HNF-1β targets were examined. The HNF-1β binding sites at the promoters of Cdh16 and Pcbd1 overlapped with open chromatin (Fig. S1). The association between HNF-1β binding and ATAC-seq peaks was statistically significant (p < 2x10− 16, odds ratio = 16.5, hypergeometric test).
Next, we examined the colocalization of ATAC-seq and HNF-1β peaks in different genomic regions. Chromatin accessibility was greater at HNF-1β binding sites in promoter regions compared to intragenic and distal intergenic regions (Fig. 2a). The HNF-1β enrichment was equal in peaks located in the promoter, intragenic, and distal intergenic regions, so the difference in ATAC-seq signal was not simply due to differences in HNF-1β binding. Instead, these findings may suggest that binding of HNF-1β preferentially increases chromatin accessibility in promoters. One possible mechanism involves histone acetylation. Acetylation of Lys27 on histone H3 (H3K27ac) reduces the electrostatic interaction with DNA and thereby increases chromatin accessibility. Moreover, the C-terminal domain of HNF-1β, which is required for transcriptional activation, has previously been shown to interact with CBP and P/CAF, coactivators that have intrinsic histone acetylase activity [7, 8]. Therefore, we tested the relationship between HNF-1β binding and histone acetylation in chromatin from embryonic kidneys. Data from a ChIP-seq analysis of H3K27 acetylation in wild-type E14.5 mouse kidney were downloaded from the ENCODE database [29]. Peaks of H3K27ac enrichment were identified, and the locations were compared with HNF-1β binding peaks. We identified enrichment of H3K27ac at 1,970 HNF-1β binding peaks, representing 24% of all HNF-1β binding peaks. HNF-1β binding sites at the promoters of Cdh16 and Pcbd1 overlapped with H3K27ac sites (Fig. S1). The association between HNF-1β binding and H3K27 acetylation was statistically significant (p < 2x10− 16, odds ratio = 19.6, hypergeometric test). Similar to the ATAC-seq results, H3K27ac enrichment was higher in promoters compared to intragenic and distal intergenic regions (Fig. 2a).
To further confirm the functional significance of HNF-1β binding, we examined the relationship between HNF-1β binding and histone H3 Lys4 trimethylation (H3K4me3), an epigenetic mark associated with gene activation [30]. Using a similar approach as above, we found enrichment of H3K4me3 at 1,002 HNF-1β ChIP-seq peaks in E14.5 mouse kidney. HNF-1β binding sites at the promoters of Cdh16 and Pcbd1 overlapped with H3K4me3 sites (Fig. S1). The association between HNF-1β binding and H3K4me3 was statistically significant (p < 2x10− 16, odds ratio = 9.2, hypergeometric test). Enrichment of H3K4me3 was highest at HNF-1β binding peaks in promoter regions and was not detected at HNF-1β binding sites in intragenic and distal intergenic regions (Fig. 2a). Next, we examined the relationship between HNF-1β binding and histone H3 Lys4 monomethylation (H3K4me1), a histone modification associated with primed enhancers [30, 31]. We found 1205 H3K4me1 ChIP-seq peaks that overlapped with HNF-1β ChIP-seq peaks. The association between HNF-1β binding and H3K4me1 was statistically significant (p < 2x10− 16, odds ratio = 7.6, hypergeometric test). In contrast to H3K4me3, H3K4me1 colocalized with HNF-1β binding peaks in intragenic and distal intergenic regions in addition to promoter regions (Fig. 2a). Intragenic HNF-1β binding sites at Gfra1 and Wnt9b overlapped with open chromatin, H3K27ac, and H3K4me1, suggesting these are enhancer sites (Fig. S1). The distribution of H3K4me3, H3K27ac, and H3K4me1 at HNF-1β binding sites was bimodal, consistent with the sliding of nucleosomes due to transcription factor binding [32]. The overlap between HNF-1β binding, ATAC-seq, and each histone modification is summarized in Supplementary Fig. 2.
Taken together, HNF-1β binding sites in E14.5 kidneys are associated with regions of open chromatin and activating histone modifications. The colocalization of HNF-1β binding and H3K27ac correlated with ATAC-seq and showed highest enrichment in gene promoters and lower signal in intragenic and distal intergenic regions. Consistent with this finding, colocalization with H3K4me3 was also highest in gene promoters, indicating that HNF-1β binding is associated with promoter activation. In contrast, colocalization with H3K4me3 was absent in intragenic and distal intergenic regions. Instead, HNF-1 peaks in these regions colocalized with H3K4me1, indicating that they corresponded to active enhancers.
HNF-1β regulates gene expression in purified UB cells
To investigate the relationship between HNF-1β binding and HNF-1β-dependent gene regulation, we performed RNA-seq analysis on UB cells purified from E14.5 mouse kidneys. To determine the dependence on HNF-1β, studies were performed in HNF-1β mutant mice and wild-type mice. Using appropriate genetic crosses, we generated Hoxb7/Cre;Hnf1bfl/fl;tdT (HNF-1β mutant) mice in which Cre recombinase was expressed under the control of the Hoxb7 promoter resulting in inactivation of HNF-1β and activation of a tdTomato reporter gene exclusively in the ureteric bud. Hoxb7/Cre;Hnf1bfl/+;tdT mice were used as controls. Consistent with previous studies [11], HNF-1β mutant kidneys were hypoplastic and had reduced UB branching compared to control kidneys (Fig. S3-S4). TdTomato was expressed exclusively in the ureteric bud in both mutant and control kidneys, which enabled purification of the UB cells by flow cytometry. Kidneys were dissected from five HNF-1β mutant kidneys and five control kidneys at E14.5, and single cell suspensions were prepared by treatment with collagenase. TdTomato-positive UB cells were purified by FACS. Representative analyses are shown in Supplementary Figs. 3 and 4. In total, UB cells accounted for 10.1 ± 1.5% of live cells in wild-type kidneys and 4.6 ± 1.2% of live cells in mutant kidneys. RNA was extracted from the purified UB cells, and cDNA libraries were prepared and sequenced. Each library contained more than 16 million reads (Table S3). All differentially expressed genes had at least 25 mean read counts in one or both genotypes.
UB cells isolated from control mice expressed known markers of the ureteric bud, including Hoxb7, Pax2, Pax8, Ret, Gfra1, and Sox9, which validated the preparation of the UB cells (Table S4). UB cells isolated from HNF-1β mutant mice showed a 94% reduction in HNF-1β mRNA levels (Table S4). Comparison of the transcriptomic profiles from HNF-1β mutants and controls identified 1,632 genes that were more highly expressed in control UB cells compared to HNF-1β mutant UB cells (log2fc≤-0.5, FDR < 0.05, negative binomial test with edgeR). We defined these genes as being activated by HNF-1β. Of these genes, 485 (29.7%) contained nearby HNF-1β binding sites, indicating that they were likely directly activated by HNF-1β. As a positive control, we identified genes that had previously been shown to be activated by HNF-1β, including Wnt9b, Cys1, Pkhd1, Pcsk9, Kif12, and Cdh16. Conversely, 2,223 genes had significantly lower expression in control UB cells compared to HNF-1β mutant UB cells (log2fc ≥ 0.5, FDR < 0.05, negative binomial test with edgeR). We defined these genes as being repressed by HNF-1β. Of these genes, 526 (24.1%) contained nearby HNF-1β binding sites, indicating that they were likely directly repressed by HNF-1β (Fig. 2b). Included in this list were Socs3 and Ccdc80, which we previously showed were repressed by HNF-1β. Fold changes in gene expression and read counts for each sample can be found in Supplementary Table 4. Taken together, these results identified a subset of genes that were directly regulated by HNF-1β in the UB in vivo.
Histone modifications and increased chromatin accessibility are associated with HNF-1β-dependent gene regulation
Next, we examined whether genes that are activated by HNF-1β have differences in chromatin accessibility or histone modifications compared with repressed genes. Comparison of the ATAC-seq signal at HNF-1β peaks showed that both activated and repressed genes associated with open chromatin. However, the enrichment of open chromatin was 30% higher at intragenic HNF-1β binding sites and 21% higher at distal intergenic sites in activated genes compared with repressed genes (Fig. 2c, Table S5). These differences were statistically significant (p = 2.3x10− 09 and 1.5x10− 03, respectively, Wilcoxon test). Similarly, the enrichment of H3K27ac was 37% higher at intragenic HNF-1β binding sites and 33% higher at distal intergenic sites in activated genes compared with repressed genes (p = 2.6x10− 19 and 3.5x10− 04, Wilcoxon test). The enrichment of H3K4me1 was 25% higher at intragenic sites and 15% higher at distal intergenic sites in activated genes compared with repressed genes (p = 3.9x10− 15 and 2.3x10− 02, Wilcoxon test). There were no significant differences between activated and repressed genes in ATAC-seq, H3K27ac, and H3K4me1 enrichment at HNF-1β binding sites in promotor regions (Fig. S4, Table S5).
The relationship between two repressive histone marks, H3K27me3 and H3K9me3, and HNF-1β binding was also examined. Only 136 H3K27me3 peaks overlapped with HNF-1β peaks, which was statistically significant (p < 2x10− 16, hypergeometric test) but had a lower odds ratio (1.97) and much less overlap than ATAC-seq and activating histone modifications. Twenty H3K9me3 peaks overlapped with HNF-1β peaks, which was not statistically significant (p = 1, odds ratio = 0.54, hypergeometric test). There were no differences in H3K27me3 and H3K9me3 enrichment at HNF-1β binding sites in activated and repressed genes compared to genes with no change in expression (Fig. S5). Taken together, HNF-1β intragenic and distal intergenic peaks showed increased chromatin accessibility and activating histone marks in activated genes compared to repressed genes.
HNF-1β peaks are associated with greater chromatin accessibility in UB-derived cell types in the developing kidney
Next, we determined whether the relationship between HNF-1β binding and chromatin accessibility seen in the embryonic kidney was observed in the UB. We took advantage of a single nuclear (sn)ATAC-seq analysis performed on postnatal day (P)0 mouse kidney [33], an age at which kidney development is still ongoing. The enrichment of snATAC-seq signal from the UB-derived principal cells (PC) and intercalated cells (IC) was compared with immune cells and nephron progenitors (NP), which do not endogenously express HNF-1β. The enrichment of open chromatin at HNF-1β binding sites in PC was 47% higher than immune cells (p = 1.3x10− 54, Wilcoxon test) and 41% higher than NP (p = 3.5x10− 13, Wilcoxon test) (Fig. 3a). The enrichment of open chromatin at HNF-1β binding sites was 47% higher in IC than immune cells (p = 1.2x10− 29, Wilcoxon test) and 41% higher than NP (p = 0.0023, Wilcoxon test) (Fig. 3b, Table S6).
We then compared the snATAC-seq signal at the HNF-1β binding sites in activated and repressed genes. The enrichment of open chromatin at HNF-1β binding sites in PC was 48% higher in activated genes than repressed genes (p = 8.5x10− 26, Wilcoxon test). The enrichment of open chromatin at HNF-1β binding sites in IC was 43% higher in activated genes than repressed genes (p = 3.6x10− 24, Wilcoxon test). Taken together, the association of ATAC-seq signal with HNF-1β binding sites seen in embryonic kidney was also observed in cells derived from the UB. Moreover, there was greater enrichment of open chromatin in genes that were activated by HNF-1β compared to repressed genes.
HNF-1β regulates axon guidance genes in the developing kidney
To identify novel pathways that are regulated by HNF-1β in the developing mouse kidney, we performed KEGG analysis on the genes that were in closest proximity to the HNF-1β binding sites in E14.5 kidneys (Fig. 4a). Axon guidance was identified as a significant canonical pathway (adjusted p-value = 7.13x10− 09). The related Rap1 signaling pathway, which is involved in intracellular signaling for plexin axon guidance receptors [34], was also identified. Axon guidance is a repulsive and attractant chemotaxis pathway that was first described in axon growth cones. More recently, axon guidance has been implicated in branching morphogenesis in epithelial organs, including the kidney [35]. We found that HNF-1β binds near 68 axon guidance genes in the developing kidney. Of these genes, 34 encode ligands or receptors in the Semaphorin/Neuropilin/Plexin, Ephrin/Ephrin receptor, Slit/Robo, and Netrin/Unc families. Pathway analysis of RNA-seq data showed that axon guidance was also highly represented among the genes that were directly activated or repressed by HNF-1β (Fig. S6). Gene set enrichment analysis (GSEA) showed that the GO Axon Guidance pathway was enriched in the RNA-seq data from E14.5 UB cells (Enrichment score=-0.45, FDR q-value = 0.0) (Fig. 4b). RNA-seq counts were ranked from the genes with the highest to lowest expression in HNF-1β mutant UB cells compared to wild-type UB cells. The enrichment score of -0.45 indicated that the genes contributing most to the enrichment of the axon guidance pathway were genes that had lower expression in HNF-1β deficient cells, i.e. activated by HNF-1β. To confirm these findings, pathway analysis was performed on genes that were differentially expressed in HNF-1β-deficient mIMCD3 cells compared to wild-type mIMCD3 cells [36]. Mutant mIMCD3 cells contained the same deletion of Hnf1b exon 1 as HNF-1β mutant mice. Pathway analysis showed that axon guidance was the second most highly enriched pathway in HNF-1β mutant cells compared to wild-type cells (Table S7). These results indicate that HNF-1β regulates the transcription of axon guidance genes both in vitro in mIMCD3 cells and in vivo in developing kidney.
The changes in expression of HNF-1β gene targets in the ephrin, semaphorin, and Slit/Robo signaling pathways were investigated in greater detail. Nrp1, Sema3c, Sema3d, Sema6a and Slit2 were significantly activated by HNF-1β, whereas Efna1, Epha3, Epha4, Epha7, Ntn4, Plxna2, Sema3a, Sema4b, Slit3, Srgap1, Unc5c and Unc5d were significantly repressed by HNF-1β (Fig. 4c). To validate RNA-seq results and investigate spatial expression of axon guidance genes bound by HNF-1β, fluorescent RNAscope in situ hybridization was performed on sagittal sections E14.5 wild-type and HNF-1β mutant kidneys. Kidney sections were co-stained with antibodies against RFP (TdTomato) to visualize the UB. Efna1 was expressed in the wild-type UB and surrounding MM, and its expression increased in the mutant UB and MM (Fig. 5a). Nrp1 was expressed in the wild-type UB and MM. In the HNF-1β mutant kidney, expression of Nrp1 decreased in the UB but remained unchanged in the MM, reflecting the UB-specific knockout of HNF-1β (Fig. 5b). Sema3c, Sema3d, and Slit2 were highly expressed in the wild-type UB and weakly expressed in the surrounding MM. In the HNF-1β mutant kidney, expression was decreased in the UB but remained unchanged in the MM (Figs. 5c, 5d, 5f). Sema6a was expressed in wild-type UB and MM; expression decreased in the mutant UB and increased in the mutant MM (Fig. 5e).
Next, we examined the relationship between HNF-1β binding, ATAC-seq, and histone modifications at Efna1, Nrp1, Sema3c, Sema3d, Sema6a, and Slit2. Efna1, which was repressed by HNF-1β, contained an intragenic HNF-1β peak that did not overlap with any of the histone modifications associated with activation (Fig. S7). Nrp1, which was activated by HNF-1β, contained an intragenic HNF-1β peak in intron 13 colocalized with H3K27ac, H3K4me1, and ATAC-seq enrichment (Fig. S7). Sema6a contained an intragenic HNF-1β binding site in its first intron that overlapped with H3K27ac and H3K4me1 peaks (Fig. S7). Sema3c did not contain an HNF-1β peak in the promoter, but there was an HNF-1β motif near the promoter and ATAC-seq, H3K4me3, and H3K27ac peaks were present near the promoter (Fig. S7). Sema3d did not contain an HNF-1β peak in the promoter, but ATAC-seq, H3K4me3, and H3K27ac peaks were present near the promoter (Fig. S7). Slit2 did not have an HNF-1β binding site in the promoter, but ATAC-seq and H3K4me3 peaks were present near the promoter (Fig. S7). Taken together, HNF-1β regulates the expression of axon guidance genes in the UB and can either activate or repress gene transcription. Promoters of activated axon guidance genes also overlapped with ATAC-seq, H3K4me3 and H3K27ac peaks. Some intragenic HNF-1β peaks overlapped with both H3K27ac and H3K4me1 peaks.