Functional genomics analysis of developing zebrafish and human endoderm reveals highly conserved cis-regulatory modules acting during vertebrate organogenesis

ATAC-seq reveals CRMs functioning in distinct sox17-expressing lineages during zebrafish embryogenesis

To enrich for endodermal cells during endoderm organogenesis, we exploited sox17:GFP fish (Mizoguchi et al. 2008). Although endogenous sox17 is rapidly downregulated at the end of gastrulation (Alexander and Stainier 1999), GFP protein persists throughout the endoderm for days after the endogenous gene has been silenced (Fig. 1A). However, sox17 expression is also observed in erythroid and endothelial lineages, limiting the utility of sox17:GFP alone for enrichment of endoderm (Chung et al. 2011). We therefore crossed homozygous sox17:GFP fish with fish homozygous for both the gata1a:dsRed and kdrl:mCherry transgenes (Fig. 1B). We then used fluorescence-activated cell sorting (FACS) to separate sox17+ vascular and erythroid mesodermal lineages (termed sox17M) from endoderm-enriched sox17+ lineages (sox17E) and all other sox17- lineages (sox17N) at 28 and 48 h postfertilization (hpf) (Fig. 1C; Supplemental Fig. S1). We chose 28 hpf as the earliest time point to ensure robust detection of all three fluorescent proteins. To verify that our sorting strategy enriched for endodermal cells in the sox17E population, we performed bulk RNA-seq on 28-hpf samples. This revealed strong enrichment for transcripts expressed in specific endoderm cell populations including pancreatic, liver, and intestinal markers (Fig. 1D; Supplemental Fig. S2; Supplemental Table 1). However, we also note enrichment of markers of fin epithelia and notochord, which we attribute to transdifferentiation of the sox17-expressing left–right organizer (KV), which gives rise to posterior cell types once its role is complete (Ikeda et al. 2022; Supplemental Fig. S3). Nevertheless, the sox17E population shows strong enrichment for endoderm markers, thus validating our approach. We therefore proceeded to perform independent duplicate ATAC-seq analysis on the three distinct FACS populations at 28 and 48 hpf (Supplemental Figs. S4–S7).

Figure 1.

Fluorescence-activated cell sorting separates endodermal and mesodermal cell populations arising from sox17+ progenitors. (A) Widefield images of sox17:GFP embryos at 28 and 48 hpf with distinct endodermal structures indicated. (B) Widefield images of coexpression of sox17:GFP with kdrl:mCherry and gata1a:dsRed transgenes in trunk and tail vasculature at 48 hpf. (C) Sorting strategy. Also see Supplemental Figure S1. (D) RNA-seq heatmap of the top 20 significant genes showing strongest fold enrichment in sox17E over sox17N at 28 hpf.

Open in new tabDownload PowerPointLink to figure

ATAC-seq analysis revealed thousands of differentially accessible regions (DARs) distinguishing the three sorted populations both within and between time points (Fig. 2A,B; Supplemental Tables 2, 3). Chromatin accessibility profiles were consistent with the predicted cell identities being enriched in each sorted cell population. For example, the embryonic hemoglobin gene cluster, erythroid marker gata1a, and endothelial marker fli1b show enhanced accessibility in the sox17M population (Fig. 2C; Supplemental Fig. S8), whereas markers of the posterior foregut (PFG; gata6 and foxa3), pancreas and duodenum (pdx1), and liver and intestine (fabp2) show enhanced accessibility in the sox17E population (Fig. 2D; Supplemental Fig. S8). However, we note that many more regions show enhanced accessibility in sox17E over sox17M, compared with sox17E over sox17N (Fig. 2A). This suggests that although there are differences in accessibility profiles that distinguish the sox17E population from sox17M, many regions of accessible chromatin in sox17E are nevertheless shared with cell types in the sox17N population rather than being unique to endoderm.

Figure 2.

Sorted cell populations have distinct chromatin accessibility profiles indicative of constituent cell identities. (A) K-means clustered heatmaps of relative ATAC-seq read densities at differentially accessible regions (DARs; FDR ≤ 0.05) for the comparisons indicated. Clusters are rank-ordered from greatest to least accessibility in the sox17E population. (B) Venn diagrams indicating number of sox17E > sox17N DARs overlapping sox17E > sox17M DARs at 28 and 48 hpf. Venn diagrams indicating overlaps of called peaks used in the analysis are depicted in Supplemental Figure S6 and genomic distribution of DARs in Supplemental Figure S7. (C) Example loci indicating greater accessibility of erythrocyte marker genes (gata1 and hemoglobin gene cluster) and the endothelial marker fli1b in the sox17M population. (D) Example loci indicating greater accessibility of markers of the posterior foregut (gata6 and foxa3), pancreas and duodenum (pdx1), and liver and intestine (fabp2) in the sox17E population compared with sox17M and sox17N. Peak heights in counts per million reads (CPM) are indicated. Significant DAR comparisons in panel D and outlined in red. For zoomed-out tracks rescaled per condition to the strongest local peak, confirming the enrichment shown here, see Supplemental Figure S8.

Open in new tabDownload PowerPointLink to figure

To more globally ascertain whether the ATAC-seq profiles of the sorted populations are consistent with the expected constituent cell populations, we performed anatomical enrichment analysis on promoter regions that showed enhanced accessibility in the sox17E population. We specifically focused on promoter regions to avoid errant annotation of DARs to putative target genes. At the stages studied, the endoderm predominantly consists of epithelial cells of the pharyngeal endoderm and intestinal rod, as well as the developing pancreas and liver primordia. As expected, we found strong enrichment for anatomical terms consistent with epithelial structures, gut, liver, and pancreas, for DARs more accessible in both sox17E > sox17M and sox17E > sox17N, especially at 48 hpf (Fig. 3A; Supplemental Tables 4–7). This is consistent with greater development of endodermal organs by this stage. We also found enrichment for notochord markers in sox17E > sox17M, consistent with transdifferentiated cells from KV incorporating into this structure (Supplemental Figs. S3, S9). Importantly, we find no such enrichment for endodermal terms for sox17N > sox17E or sox17M > sox17E DARs, confirming enrichment of endoderm in the sox17E populations (Supplemental Fig. S10). Conversely, we found strong enrichment for anatomical terms consistent with nervous system and eye development but not endodermal tissues for sox17N > sox17E promoter DARs (Supplemental Table 8). This is highly consistent with these structures being absent from the sox17E population, as expected. However, we also found strong enrichment for markers of neural structures in sox17E > sox17M and the lateral line system in both sox17E > sox17M and sox17E > sox17N (Fig. 3A). This suggests either that cells arising from sox17+ endodermal progenitors and KV in the sox17E population exhibit similar gene accessibility (and consequent potential expression) signatures to the nervous and lateral line systems or that there is hitherto unrecognized sox17:GFP expression in these cell types. A consideration, however, is that promoter accessibility is not necessarily correlated with promoter activity (Starks et al. 2019). For example, a promoter may be accessible but actively silenced by transcriptional repressors. It is therefore possible that accessibility of neural gene promoters is not indicative of neural cells in the sox17E population.

Figure 3.

Promoter accessibility within the sox17E populations is consistent with presence of endodermal cell populations and lateral line neurons. (A) Heatmap of −log₁₀(adjusted P-value) from fishEnrichr Anatomy GeneRIF Predicted Z-score analyses of promoters showing greater accessibility in the comparisons indicated (FDR ≤ 0.05) as annotated by ChIPseeker. The key provides a continuous scale of −log₁₀(adjusted P-value), with minimum and maximum significance indicated, along with the 50th percentile significant value. Selected terms are shown; complete fishEnrichr outputs are in Supplemental Tables 9 and 10. (B) Dorsal widefield images of sox17:GFP embryo at 48 hpf. Nine focal layers have each been artificially colored based on depth from dorsal to ventral as indicated. (C) 3D rendering of confocal images 48 hpf sox17:GFP laterally orientated embryo with head to the left and tail to the right. Scale bar, 100 µm. (OC1) Organ of Corti 1; (plg) posterior lateral line ganglion; (primD) dorsal primordium; (primII) second primordium; (iv) intersegmental vessel; (L) liver; (P) pancreas; (dC) cells within duct of Cuvier; (pe) pharyngeal endoderm.

Open in new tabDownload PowerPointLink to figure

To further explore whether there is GFP expression in the nervous and lateral line systems, we performed detailed imaging of sox17:GFP embryos. This revealed subtle but detectable expression in a small set of neurons exhibiting morphology consistent with being a subset of lateral line neurons within the hindbrain and trunk (Fig. 3B,C). This includes organ of Corti 1 (OC1), posterior lateral line ganglion (plg), dorsal primordium (primD), and second primordium (primII) (Pujol-Martí and López-Schier 2013). We otherwise did not detect GFP expression in developing neural structures. Development of the posterior lateral line ganglion, nerve, and its projections has been shown to be affected by inhibition of retinoic acid (RA) signaling (Begemann et al. 2004). Consistent with expectations, we found that modulation of RA signaling through addition of RA itself, an inhibitor of Aldh1a mediated RA synthesis (diethylaminobenzaldehyde [DEAB]), and an inverse agonist of RA receptor (BMS493) all disrupted the formation of these sox17:GFP+ neurons (Supplemental Fig. S11). We therefore conclude that the sox17:GFP reporter is expressed in a subset of lateral line neurons but not other neural structures. However, to the best of our knowledge, endogenous sox17 expression has not been reported in the zebrafish lateral line, either in literature providing detailed in situ hybridization analysis of sox17 (Chung et al. 2011) or in whole-embryo single-cell RNA-seq time course data sets including Daniocell (Sur et al. 2023). These resources do, however, report endoderm and hematovascular sox17 expression. It is therefore possible that sox17:GFP reporter expression in the lateral line neurons is an artifact caused by the position of transgene integration into the genome rather than reflecting a true sox17 expression domain.

Overall, we conclude that we have profiled chromatin accessibility across two distinct sox17-expressing populations: sox17M containing erythroid progenitors and a subset of endothelial cells, as well as sox17E containing all endoderm, plus KV-derived posterior cell types and lateral line neurons. Given sox17E captures endodermal chromatin accessibility and thus active CRMs (albeit accessible sequences are not always active CRMs), we next wanted to determine which of these CRMs are likely to act during endoderm formation across the 400 million years of evolution separating zebrafish and humans.

Highly conserved CRMs functioning in both human and zebrafish endoderm cluster around genes controlling diverse aspects of endodermal- and vertebrate-specific development

To determine which potential CRMs captured by our ATAC-seq data are recognizably functional in human endoderm, we compared our zebrafish data to existing ChIP-seq data for the marker of active promoters and enhancers, H3K27ac (Creyghton et al. 2010). Specifically, we utilized H3K27ac ChIP-seq data from hESCs that had undergone efficient directed differentiation to either anterior foregut (AFG), PFG, or midgut/hindgut (MHG) (Supplemental Table 11; Loh et al. 2014). AFG principally gives rise to anterior structures including thyroid and lungs in human, PFG to liver and pancreas, and MHG to small and large intestine.

We wanted to avoid eliminating DARs that are potentially functionally important in the endoderm, while enriching for genomic regions that are not likely to be controlling constitutively expressed genes. We therefore primarily focused our attention on DARs showing greater accessibility in sox17E relative to sox17M (sox17E > sox17M DARs). This is because our initial analyses indicate that relatively few DARs distinguish sox17E from the sox17N population but that the majority that do also distinguish sox17E from sox17M (Fig. 2A,B).

To identify zebrafish sox17E DARs corresponding to human H3K27ac ChIP-seq peaks, we compared both data sets to HCNEs exhibiting ≥70% sequence identity across 30 or more alignment columns (Engström et al. 2008). We identified a total of 5000 HCNEs overlapping H3K27ac peaks in at least one of the AFG, PFG, and MHG cell populations (Fig. 4A); 7435 HCNEs overlap sox17E > sox17M DARs, whereas 236 HCNEs overlap sox17E > sox17N DARs, of which 112 are common to both sox17E > sox17M and sox17E > sox17N DARs (Fig. 4B). Of the 5000 HCNEs overlapping human H3K27ac peaks, 1701 were also identified in sox17E > sox17M DARs (Fig. 4C,D; Supplemental Figs. S12 and S13; Supplemental Tables 12, 13). These HCNEs occur proximal to genes with prominent roles in endoderm development such as the regulator of gut development DLL1/dld (Fig. 4E; Pellegrinet et al. 2011; Troll et al. 2018), regulator of pancreas development FOXA2/foxa2 (Fig. 4F; Shin et al. 2008; Lee et al. 2019; Elsayed et al. 2021), and pancreatic and biliary regulator HES1/her6 (Fig. 4G; Spence et al. 2009).

Figure 4.

HCNE CRMs bearing the hallmarks of functionality in human and zebrafish endoderm are proximal to genes with deleterious pancreatic, neural, and craniofacial phenotypes. (A) Number of HCNEs overlapping human H3K27ac peaks in AFG, PFG, and/or MHG patterned endoderm cell populations derived from hESCs. Number of HCNEs not overlapping H3K27ac peaks are at the bottom right. (B) Distribution of HCNEs overlapping sox17E > sox17M DARs versus sox17E > sox17N DARs at 28 and 48 hpf. Number of HCNEs not overlapping DARs are at the bottom right. (C) Distribution of HCNEs overlapping sox17E > sox17M DARs also overlapping human H3K27ac peaks in AFG, PFG, and/or MHG patterned endoderm cell populations derived from hESCs. Number of HCNEs not overlapping H3K27ac peaks or DARs are at the bottom right. (D) Schematic of human HCNE subsets used for functional, anatomical, and phenotypic enrichment analyses using GREAT. (E–G) Example orthologous loci bearing zebrafish HCNE DARs and human H3K27ac HCNE peaks. Normalized peak heights in CPM are indicated. HCNEs overlapping peaks and DARs are displayed below and are color-coded as the tracks. Gray boxes and lines indicate regions with HCNE matches between species. (H) Heatmap of enrichment of single-gene mouse knockout phenotypes for collections of HCNEs as indicated in D. Significance in −log₁₀(HyperFdrQ) is indicated both relative to the whole genome and to all HCNEs to account for nonrandom distribution of HCNEs. Discrete color-coded values in the key indicate minimum and maximum significance, with 25th, 50th, and 75th percentiles between, as indicated on the continuous scale. Terms are color-coded according to anatomical categories as indicated.

Open in new tabDownload PowerPointLink to figure

To more globally analyze the targets of putative CRMs acting in both human and zebrafish endoderm, we used the genomic regions enrichment of annotations tool (GREAT) to perform functional and phenotypic annotation analyses on the common 1701 HCNEs (McLean et al. 2010). Annotation of putative target genes used the “basal plus extension” method in which each gene is assigned a basal regulatory domain of 5 kb upstream of and 1 kb downstream from the transcription start site (regardless of other nearby genes). The gene regulatory domain is extended in both directions to the nearest gene's basal domain but no more than 1000 kb in one direction. To determine whether the HCNEs acting in AFG, PFG, and MHG endoderm show distinct characteristic distribution patterns indicative of alternative target genes, we analyzed subsets unique to each population as well as the complete combined set (Fig. 4D). Furthermore, because HCNEs are known to have nonrandom genomic distribution patterns with specific biases toward trans-dev genes, we tested whether the 1701 HCNEs (hereafter referred to as “endoderm HCNEs”) showed enrichment around specific classes of genes (categorized according to function of gene products and their loss-of-function phenotypes) both overall and compared with all HCNEs (Fig. 4H; Supplemental Figs. S14–S18).

Consistent with expectations, our analyses reveal a highly significant association between endoderm HCNEs and genes encoding DNA-binding and gene regulatory proteins, including HMG domain-containing and chromatin binding factors. Notably, however, the endoderm HCNEs appear to show a greater enrichment of such genes compared with all HCNEs (Supplemental Fig. S14). Furthermore, both biological process Gene Ontology (GO) terms, and mammalian phenotype terms associated with endoderm HCNE subsets from AFG, PFG, and MHG cells are broadly consistent with their anterior–posterior location within the vertebrate body plan. For example, putative target genes of HCNEs shared between AFG H3K27ac peaks and sox17E DARs are significantly associated with GO biological process terms, mouse phenotypes, and human phenotypes focused on anterior structures including the brain, ear, and craniofacial structures (Fig. 4H; Supplemental Figs. S15, S17, S18). Similarly, putative targets of PFG H3K27ac peaks and sox17E DARs are associated with pancreas development and pancreatic abnormalities in knockout mice, as well as cardiac defects in humans (Fig. 4H; Supplemental Fig. S18). Conversely, putative targets of MHG H3K27ac peaks and sox17E DARs appear to be largely associated with formation of tube and ductal structures such as neural tube closure, endolymphatic duct, and kidney formation, as well as vertebral column formation (Fig. 4H; Supplemental Figs. S15, S17, S18). This may be indicative of common gene regulatory programs governing tube formation, including the gut, and CRM accessibility potentially being broadly coordinated regionally along the anterior–posterior axis.

We note that many of the putative targets of HCNE CRMs identified in human and zebrafish endoderm function in nervous system development (Fig. 4H; Supplemental Figs. S15, S17, S18). For example, FOXA2/foxa2 functions both in pancreas development and in the formation of the floor plate of the neural tube (Ang and Rossant 1994; Brand et al. 1996; Shin et al. 2008; Dal-Pra et al. 2011). Similarly, HES1/her6 has been shown to exhibit nervous system defects, including premature neurogenesis, severe neural tube defects, increased numbers of pulmonary neuroendocrine cells, and also pancreatic hypoplasia in mouse knockout models (Ishibashi et al. 1995; Ito et al. 2000; Jensen et al. 2000), and regulates cell proliferation in the hindbrain in zebrafish (Coolen et al. 2012). As we discuss later, it is therefore highly likely that the significant enrichment for terms associated with the nervous system can be attributed to common neuro-pancreatic gene regulatory programs operating across the developing nervous system and pancreas. Overall, we conclude that we have identified a compendium of HCNEs likely to control gene expression during early organogenesis across vertebrate evolution.

HCNEs bearing hallmarks of functionality in human and zebrafish endoderm are enriched for binding sites of endoderm TFs

To determine which TFs are likely to be acting via the endodermal HCNEs, we performed de novo and known TFBS enrichment analysis in the human and zebrafish sequences (Fig. 5A–C). We particularly focused on the PFG population given its role in pancreas development and the strong association between putative HCNE target genes and pancreas development in our previous analyses (Fig. 4H). TFBSs are often degenerate, with multiple TFs potentially able to bind the same site. Indeed, many of the TFBSs we identified by our analysis are highly similar (Supplemental Fig. S19). We therefore also analyzed expression of the candidate TFs corresponding to the enriched TFBSs in RNA-seq data from the human AFG, PFG, and MHG populations (Loh et al. 2014). This revealed expression of some TFs in the same endoderm populations exhibiting enrichment of their corresponding TFBSs in HCNEs (Fig. 5C). For example, TFBSs for SOX21 and MEIS1 show greater enrichment in AFG and PFG HCNEs than in MHG HCNEs, and SOX21 and MEIS1 are also more highly expressed in AFG and PFG. Similarly, HOXA11, HOXD11, and HOXA13 TFBSs show greater enrichment in MHG HCNEs, consistent with greater expression of the TFs in MHG endoderm. However, we also find strong enrichment for TFBSs of multiple TFs known to have key functions in the endoderm despite the TFs showing no/minimal expression in the AFG, PFG, or MHG RNA-seq data sets (Fig. 5C). For example, ASCL1, ASCL2, PTF1A, PDX1, NEUROD1, and NKX6.1 are all known or suggested to have roles in pancreatic endoderm development (Offield et al. 1996; Naya et al. 1997; Krapp et al. 1998; Sander et al. 2000; Yee et al. 2001; Kawaguchi et al. 2002; Dong et al. 2008; Binot et al. 2010; Flasse et al. 2013; Duque et al. 2022; Vanheer et al. 2023; Zhu et al. 2023), and their TFBSs are enriched in PFG HCNEs, although their expression is low or not detected in PFG endoderm. However, we find strong expression of these TF genes in RNA-seq data from hESCs further differentiated beyond the PFG stage to pancreatic progenitors (PPs). This suggests that the HCNE CRMs these TFs act through bear functional marks prior to the expression of these TFs. Notably, multiple SOX family TFs, including SOX2, -3, -4, -6, -9, -15, -17, and -21, show enrichment in endoderm HCNEs and appreciable expression in PFG and PP populations, as do FOXA1 and FOXA2. Given prior evidence that SOX and FOX TFs often have pioneer activity (Kamachi and Kondoh 2013; Iwafuchi-Doi et al. 2016; Julian et al. 2017; Fuglerud et al. 2022), it is tempting to speculate that earlier expression of these TFs may render endodermal HCNEs accessible for subsequent binding by later expressed pancreatic TFs like PTF1A. Indeed, interrogation of published ChIP-seq data sets reveals a subset of endodermal HCNEs show occupancy by FOXA2 at multiple time points prior to PTF1A expression, followed by subsequent PTF1A occupancy. Notably, the summits of ChIP-seq peaks are directly at the HCNEs, strongly suggesting the HCNEs are key to TF recruitment (Supplemental Fig. S20).

Figure 5.

Motif enrichment analysis reveals candidate TFs acting via endoderm HCNE CRMs. (A,B) High-confidence de novo motifs enriched in HCNEs overlapping both human PFG H3K27ac peaks and sox17E > sox17M DARs at 28 and/or 48 hpf for the human (A) and zebrafish (B) HCNE sequences. P-values and closest match vertebrate TFs assigned by HOMER are indicated. (C) Hierarchical clustering of enrichment of known motifs in human and zebrafish endoderm HCNEs overlapping both sox17E > sox17M DARs at 28 and/or 48 hpf, as well as human AFG, PFG, and MHG H3K27ac peaks as indicated. Enrichment in both human and zebrafish HCNE sequences is shown. RNA-seq expression of TF genes corresponding to enriched motifs is indicated as fragments Per Kilobase of transcript per Million mapped reads (FPKM) for both human and zebrafish data sets as indicated. (AFG) Anterior foregut; (PFG) posterior foregut; (MHG) mid/hindgut; (PP) pancreatic progenitors. Zebrafish have multiple orthologs of mammalian TF genes owing to the additional whole genome and other small duplications during teleost evolution. Where present, a, b, and c paralogs are depicted. (*) Although the enriched motif is FOXA2:EBOX FPKM, heatmap indicates FOXA2 expression alone.

Open in new tabDownload PowerPointLink to figure

PFG HCNEs have characteristic patterns of TFBS co-occurrence but show relatively limited evidence for consistent rigid grammatical constraint

Given we identify TFBSs of major regulators of pancreas development in PFG HCNEs identified in both species, we conclude that these highly conserved putative CRMs are likely to be important for pancreas development. That the HCNE sequences are conserved between humans and zebrafish demonstrates that arrangements of putative TFBSs have remained largely consistent across 400 million years of evolution. A constrained TFBS “grammar” (arrangement, spacing, and orientation of TFBS) has historically been suggested to point to the “enhanceosome” model of CRMs, in which a rigid TFBS grammar is required for correct assembly of TF complexes (for reviews, see Long et al. 2016; Jindal and Farley 2021). We wanted to determine whether PFG HCNEs contain consistent sets of grammatically constrained TFBSs, potentially suggestive of consistent TF complexes acting across subsets of HCNEs. We therefore analyzed both the co-occurrence of the identified enriched TFBSs within HCNEs and also whether co-occurring TFBSs showed significantly consistent spacing patterns (Fig. 6; Supplemental Fig. S21).

Figure 6.

Consistent 12 bp spacing of homeodomain–bHLH TFBS pairs is significantly enriched for a small subset of orthologous HCNEs. (A) Matrix of TFBS motifs from Figure 5 and Supplemental Figure S19 showing consistent spacing in SpaMo analysis. Cells in the upper matrices are color-coded according to best orientiation, that is, the strand and position of the secondary motif relative to the primary motif for the most significant spacing. However, significant spacing may be present in other orientations also. Numbers in cells indicate the spacing of motifs in base pairs that is significantly enriched. Numbers in black indicate significant spacing in only one of human and zebrafish; gold numbers indicate motif pairs significant for HCNEs of both species. Cells in the lower matrices are color-coded based on significance of TFBS spacing. Numbers in cells indicate the number of occurrences of that spacing in significant orientations defined by SpaMo. Homeodomain–bHLH pairs are indicated in bold black boxes. (B) Distance and orientation plots of the two TFBS pairs identified as significant in both human and zebrafish HCNEs. Quadrants within individual plots are as labeled on the figure. Red lines within graphs indicate secondary TFBS locations contributing to significant spacing from the primary TFBS (depicted in the center), and gray bars indicate all other instances of the secondary TFBSs. (C) Bar graphs indicating cumulative numbers of HCNEs containing homeodomain–bHLH TFBS pairs exhibiting 12 bp spacing. The TFBSs pairs from A in which the TFBS pairs were identified are annotated on bar segments. The numbers of orthologous HCNEs shared between segments are indicated between the bars. The genomic locations of the six HCNEs indicated are shown in Supplemental Figure S22.

Open in new tabDownload PowerPointLink to figure

There are well-characterized examples of “suboptimization” of TFBSs within developmental CRMs, in which low-affinity TFBSs tune target gene expression preventing deleterious ectopic or overexpression (Farley et al. 2015, 2016). However, position frequency matrices used for most motif analyses are derived from binding experiments (e.g., ChIP-seq) that are most indicative of high-affinity TFBSs. To avoid exclusion of potentially important low-affinity TFBSs in our analyses, we analyzed co-occurrence using both permissive and stringent cutoffs for identification of individual TFBSs. Our results indicate strong pairwise co-occurrence of specific TFs throughout PFG HCNEs, whereas some TFBS combinations are rare (Supplemental Fig. S21). As expected, patterns of TFBS co-occurrence are very similar between human and zebrafish PFG HCNEs. This suggests that the same sets of TFs are likely capable of acting via the same HCNE CRMs in both species. However, this merely indicates whether TFBSs consistently occur in the same HCNEs but not whether they show consistent spacing or orientation. We therefore next tested whether TFBS pairs have consistent arrangement and spacing, potentially indicative of constrained binding of TF complexes. Although we find 16 TF pairs showing consistent arrangement in human PFG HCNEs and 15 in the homologous zebrafish HCNEs (P < 1 × 10⁻⁵), only PTF1A-PBX2 and MEIS1-ASCL2 were identified in both species (Fig. 6A,B). Notably, the significant spacing of PTF1A-PBX2 and MEIS1-ASCL2 TFBSs was the same in both species (12 and 12 bp, respectively). PTF1A and ASCL2 are both basic helix–loop–helix (bHLH) factors with very similar DNA-binding preferences, whereas PBX2 and MEIS1 are both homeodomain TFs that also bind similar sequences to each other (Supplemental Fig. S19). Furthermore, multiple TF pairs showing significantly consistent spacing patterns in only one of the two species are homeodomain–bHLH pairs exhibiting 12 bp spacing (Fig. 6A,B). Identification of the PFG HCNEs for each of the significant 12-bp-spaced homeodomain–bHLH TF pairs reveals substantial overlaps between HCNE sets both within and between species (Fig. 6C; Supplemental Fig. S22). This suggests that subtle differences in highly similar motifs may dictate whether significance is achieved in one or both species, as well as for which of the potential homeodomain–bHLH motif pairs. It is therefore possible that multiple different combinations of homeodomain–bHLH TF pairs may be capable of binding these sites.

To determine whether TFs corresponding to 12-bp-spaced homeodomain–bHLH TFBSs pairs are coexpressed in the developing endoderm, we analyzed published zebrafish single-cell RNA-seq data (Farrell et al. 2018; Sur et al. 2023). This revealed that although the homeodomain TFs are more broadly expressed than bHLH TFs in the developing endoderm, there is nevertheless clear coexpression of multiple homeodomain–bHLH TF pairs in different endoderm subpopulations. This includes meis1a/b with tcf12 in pharyngeal pouch endoderm, ptf1a with meis1a/b in exocrine pancreas progenitors, and a ASCL2 homolog ascl1a with pdx1, pbx2, and meis1a in intestinal secretory cell progenitors (Supplemental Fig. S23). It is therefore possible that these homeodomain–bHLH TF pairs coregulate a subset of target genes in the developing endoderm. However, it is also possible that other TF pairs capable of binding the same sequences act via these sites.

The bHLH TFBSs in question have a core consensus of CAGCTG and are therefore essentially palindromic, whereas the homeodomain TFBSs have more definable directionality (Fig. 6B; Supplemental Fig. S19). Our analyses reveal that the orientation of the homeodomain TFBSs is typically (but not always) consistent relative to the bHLH TFBS (Fig. 6A,B). It is therefore possible that this orderly arrangement of TFBSs is critical to TF complex engagement or competition at these sites. However, to the best of our knowledge, no physical interactions have yet been identified for the homeodomain–bHLH TF pairs showing significantly consistent spacing.

Overall, our analyses reveal consistent patterns of TFBS co-occurrence in PFG HCNEs between humans and zebrafish, potentially indicative of conserved cooperative gene regulatory programs. However, there is only limited evidence for a highly consistent CRM grammar both across sets of HCNEs and between species, with consistent arrangement of TFBSs pairs only apparent for small subsets of HCNEs within each species, and only 12-bp-spaced homeodomain–bHLH TFBS pairs observed between species in six orthologous HCNEs. This potentially suggests either that few of the TFs with individually enriched TFBSs act in complexes with each other or that the complexes are flexible in terms of combined recognition of target TFBSs.

HNF1B HCNE CRMs identified in endoderm drive expression in the developing hindbrain

We next wanted to determine whether the HCNE CRMs we identified in human and zebrafish endoderm can drive expression in the predicted endodermal tissues. Given our observation that many HCNE CRMs identified in zebrafish and human endoderm are proximal to genes that govern both pancreas and nervous system development, we chose to focus on putative HCNE enhancers of such a gene to determine whether the HCNEs are similarly capable of driving endoderm and neural gene expression. Furthermore, increasing evidence suggests that subsets of cases of many disorders including diabetes and pancreatic agenesis are likely to be driven by genetic alteration of CRMs (Weedon et al. 2014; Claringbould and Zaugg 2021; Miguel-Escalada et al. 2022). We therefore also wanted to examine HCNE CRM function at a gene in which coding mutations are already known to cause diabetes. Notably, our analyses indicate HCNE CRMs at monogenic diabetes gene HNF1B and its zebrafish ortholog hnf1ba, both of which are expressed in PFG endoderm (Sun and Hopkins 2001; Loh et al. 2014). HNF1B haploinsufficiency causes maturity onset diabetes of the young (MODY), including pancreas atrophy, whereas homozygous disruption of hnf1ba in zebrafish leads to a failure of pancreas development during embryogenesis (Fajans et al. 2001; Sun and Hopkins 2001; Quilichini et al. 2021). However, although HNF1B haploinsufficiency in humans also causes renal defects, developmental defects in the zebrafish renal system in hnf1ba mutants are largely mitigated owing to coexpression of the hnf1ba paralog hnf1bb in the developing pronephros (Naylor et al. 2013). This is in stark contrast to the developing endoderm, in which hnf1bb is not substantially expressed and does not compensate for loss of hnf1ba (Sun and Hopkins 2001). Our sox17E ATAC-seq data indicate regions of accessibility containing HCNEs in both introns 4 and 5, while the homologous human HCNEs are within PFG H3K27ac peaks in the same introns (Fig. 7A). These HCNEs have not been maintained in hnf1bb after duplication of the ancestral hnf1b gene (Supplemental Fig. S24). Differences in expression at the hnf1ba/b loci may therefore be owing to the function of the HNF1B/hnf1ba HCNEs. Consequently, we prioritized HNF1B/hnf1ba intron 4 and 5 HCNEs for analysis, predicting that they would govern the more prominent expression of HNF1B/hnf1ba in PFG and its derivatives, which hnf1bb lacks. However, hnf1ba also exhibits a much broader hindbrain expression domain than hnf1bb (Choe et al. 2008), potentially suggesting the HCNEs present in hnf1ba but not hnf1bb may also govern hindbrain expression domains. Nevertheless, interrogation of published ChIP-seq data from hESCs differentiated in vitro to PPs indicates the HNF1B intron 5 HCNE is bound by PDX1, a key regulator of pancreatic development, suggesting a likely function in pancreas development (Supplemental Fig. S25). Furthermore, TF genes corresponding to co-occurring TFBSs in the intron 4 HCNEs show marked coexpression with hnf1ba in developing endoderm (Supplemental Fig. S26).

Figure 7.

HNF1B/hnf1ba HCNE CRMs bearing the hallmarks of functionality in zebrafish endoderm and human PFG endoderm are sufficient to drive reporter expression in the developing hindbrain but not endoderm. (A) H3K27ac ChIP-seq and ATAC-seq data at human HNF1B and zebrafish hnf1ba loci, respectively. HCNEs overlapping H3K27ac peaks and DARs are indicated in corresponding colors below the tracks. Normalized peak heights in CPM are indicated. HCNE and putative enhancer sequences cloned for reporter assays are outlined. (B) Lateral view of mCherry reporter expression driven by the i4Enh +6–8 kb enhancer at 48 hpf in sox17:GFP embryos injected with the reporter construct at the one-cell stage. Maximum intensity z-projections of 27 z-slices from confocal images. (C) 3D rendering of confocal images showing 10 z-slices at a time going through the embryo along the z-axis from the left to right side. (D) Dorsal; (R) right; (P) posterior; (V) ventral; (L) left; (A) anterior. Scale bars, 100 µm.

Open in new tabDownload PowerPointLink to figure

In addition to the intron 4 and 5 HCNEs, there is also a sox17E ATAC-seq peak 3 kb upstream of the transcription start site that does not harbor HCNEs. To determine whether these putative enhancer regions and the HCNEs are capable of driving expression in foregut endoderm, we tested their ability to drive mCherry reporter expression in sox17:GFP embryos. We tested the ability of multiple genomic regions to drive expression. For intron 4, this included the entire accessible zebrafish region (i4Enh +6–8 kb), discrete elements on the flanks of the intron 4 accessible region lacking HCNEs (i4Enh +6 kb and i4Enh +8 kb), and just the HCNE cluster (i4zHCNE). We also similarly tested the equivalent human intron 4 HCNE cluster (i4hHCNE). For intron 5, we tested discrete regions containing the zebrafish and human HCNE (i5zHCNE and i5hHCNE, respectively). We also tested the upstream accessible region (Enh –3 kb). Each reporter construct also contained the crystallin, alpha a (cryaa) promoter upstream of gfp. This drives GFP expression in the lens of the eye, thus providing a constitutive marker to control for injection. We imaged the zebrafish regularly up to and including 48 hpf, consistent with the terminal ATAC-seq time point. Except for Enh –3 kb, all constructs drove expression in multiple tissues, with the most consistent mCherry expression seen in hindbrain (Fig. 7B,C; Supplemental Figs. S27, S28; Supplemental Tables 14, 15). None of the reporter constructs yielded mCherry coexpression with sox17:GFP, indicating that they are not individually capable of driving expression in the endoderm (or lateral line neurons). All constructs containing HCNEs drove expression in the developing hindbrain (i4zHCNE, 62.2% of embryos; i4hHCNE, 87.5%; i4Enh +6–8 kb, 28.6%; i5zHCNE, 28.6%; i5hHCNE, 88.9%) (Fig. 7B,C; Supplemental Fig. S28; Supplemental Tables 14, 15). The hindbrain expression domain was typically discrete and did not overlap forebrain marker otx2b:Venus (Supplemental Fig. S29). This hindbrain expression driven by HNF1B/hnf1ba HCNEs is consistent with known expression of the endogenous genes in developing rhombomeres (Pouilhe et al. 2007; Choe et al. 2008). However, other sequences appear to be required for endoderm expression.

We note that hnf1ba intron 5 exhibits chromatin accessibility in the sox17M population, and both i5zHCNE and i5hHCNE consistently drive reporter expression in the zebrafish heart as well as the hindbrain (Fig. 7A; Supplemental Fig. S27; Supplemental Table 15). We previously showed that sox17:GFP is coexpressed with kdrl;mCherry in the developing heart, consistent with known sox17 expression in cardiac precursors in the lateral plate mesoderm (Chung et al. 2011; Johal et al. 2025). Hnf1ba expression has also been reported in the heart (Chi et al. 2008b). The i5HCNE enhancer therefore potentially regulates both cardiac and neural expression domains of hnf1ba but is insufficient for pancreas expression despite bearing the hallmarks of functionality in endoderm.

To explore why the studied reporter constructs were incapable of reproducing the hnf1ba endoderm expression pattern, we further examined both our zebrafish ATAC-seq data and published human endoderm functional genomic data (Lee et al. 2019). This revealed binding of TFs FOXA2, GATA4/6, and HNF1B in non-HCNE genomic regions within HNF1B introns 1, 4, 5, 6, and 8 and also in proximal intergenic regions (Supplemental Fig. S30). We also note FOXA2 and HNF1B binding at a HCNE in intron 1 (Supplemental Fig. S30). However, this HCNE did not correspond to a significant DAR in the zebrafish ATAC-seq data (Fig. 7).

FOXA2 is required for enhancer priming during pancreatic differentiation from hESCs (Lee et al. 2019), whereas haploinsufficiency of GATA6 and GATA4 compromise PP differentiation (Shi et al. 2017), consistent with GATA4 and GATA6 being implicated in neonatal diabetes and pancreatic agenesis (Allen et al. 2012; Shaw-Smith et al. 2014). The presence of these critical TFs at non-HCNE CRMs at HNF1B despite the HCNEs coinciding with H3K27ac, ATAC-seq peaks, and PDX1 binding (Supplemental Fig. S30) suggests that its endodermal expression domain is regulated by additional CRMs that are not highly conserved between zebrafish and humans. Interestingly, zebrafish sox17E ATAC-seq peaks are also observed in a similar location in intron 8 to human FOXA2/GATA4/GATA6 binding, as well as in additional upstream regions (Supplemental Fig. S31). It is possible that the locations of key CRMs have therefore remained roughly conserved, whereas TFBS grammar within these CRMs has not.

Overall, we identify a subset of HCNEs bearing hallmarks of functionality in developing zebrafish and human endoderm during early organogenesis. These HCNEs are strongly enriched for TFBSs of factors critical during endodermal organ development and strongly associated with genes controlling craniofacial and pancreas development. However, the putative HCNE target genes governing pancreas development also have key roles in the developing nervous system, and the HCNE enhancers at the HNF1B/hnf1ba locus are sufficient to drive expression of only the HNF1B/hnf1ba neural and cardiac, but not endodermal, expression domains. The present evidence suggests that the HCNEs bear key hallmarks of functionality during endoderm development, including accessibility and H3K27ac, but must operate within a complex regulatory landscape to influence endodermal expression.

A common neuro-pancreatic gene regulatory and cis-regulatory program?

Many similarities between endocrine cells of the pancreas and neurons have been noted over the years including production of polypeptide hormones, neurotransmitters, and their receptors, as well as common chromatin methylation signatures (Pearse and Polak 1971; Fujita et al. 1980; Gonoi et al. 1994; Maechler and Wollheim 1999; van Arensbergen et al. 2010). Furthermore, many of the TFs required for development of the endocrine pancreas have key roles in the developing nervous system throughout vertebrates, including HNF1B, FOXA2, ISL1, PTF1A, ONECUT1, NEUROD1, HES1, PAX6, SOX4, and MEIS2 (Ang and Rossant 1994; Ishibashi et al. 1995; Pfaff et al. 1996; Miyata et al. 1999; Sun and Hopkins 2001; Haubst et al. 2004; Hutchinson and Eisen 2006; Nolte et al. 2006; Hori et al. 2008; Potzner et al. 2010; Dal-Pra et al. 2011; Coolen et al. 2012; Roy et al. 2012; Machon et al. 2015; Jin and Xiang 2019). This is consistent with HCNE CRMs identified by our analyses of zebrafish and human endoderm being associated with both pancreas and nervous system development. TFBSs significantly enriched within these HCNE CRMs include those of TFs also operating during pancreas and nervous system development including PTF1A, PDX1, FOXA2, ISL1, SOX4, and NEUROD1 (Fig. 5). Given the HNF1B/hnf1ba endoderm HCNE CRMs drive expression in the developing nervous system, it is tempting to speculate that, as well as there being a common set of TFs governing endocrine pancreas and nervous system development (neuro-pancreatic TFs), these TFs also may act together within common CRMs that function in both the pancreas and nervous system. However, a broader analysis of HCNE CRMs performed at scale would be necessary to draw such conclusions. Because any given HCNE CRM is not necessarily individually sufficient to drive neural and pancreatic gene expression, future work will require both reporter assays and genetic deletion studies to test the effects on putative target genes. Combinatorial analysis of CRMs in both reporter assays and deletion studies will also likely be required to fully characterize the roles of HCNE CRMs regulating any given gene.

Common gene regulatory programs spanning pancreas and nervous system development potentially implies co-option of an ancient gene regulatory program from a common ancestral cell type. Indeed, all deuterostomes have a nervous system, but only jawed vertebrates have a true, morphologically distinct pancreas. Thus, logically a common gene regulatory program could have arisen during the evolution of the nervous system and subsequently been co-opted during pancreatic evolution. This is supported by analysis of sea urchin Stongylocentrotus purpuratus larvae, which identified a subset of neurons derived from cells with a putative “prepancreatic” signature consisting of expression of homologs of multiple neuro-pancreatic TFs including SpSoxC (homologous to SOX4) and SpPtf1a (Perillo et al. 2018). However, of the 1701 endodermal HCNE CRMs identified by our analysis of human and zebrafish data, only 24 are also HCNEs between human and lamprey genomes by the same criteria of 70% identity across at least 30 alignment columns. It therefore seems likely that although combined action of neuro-pancreatic TFs in the nervous system predates emergence of the vertebrate lineage, the HCNE CRMs they are predicted to bind are likely to have evolved after the divergence of jawless and jawed vertebrates.

Individual HCNE CRMs appear to have a distinct regulatory logic

Highly conserved CRMs have been suggested to be constrained by the action of TF complexes requiring consistent configuration of TFBSs (Guturu et al. 2013; Long et al. 2016; Jindal and Farley 2021). Our analysis of PFG HCNE CRMs suggests highly conserved patterns of TFBS co-occurrence (Supplemental Fig. S21). However, relatively few TF pairs show significant enrichment for specific TF spacing across different HCNEs within the same species (Fig. 6A). Indeed, although we identify a significant incidence of 12 bp spacing of putative homeodomain–bHLH binding sites, this was only detected in a common set of six HCNEs. This suggests either that the selective pressure maintaining HCNEs is not based on the action of TF complexes requiring rigid TFBS spacing or that TF complexes requiring such rigid grammar do not act widely across HCNE CRMs. It is similarly possible that depending on the configuration of DNA-binding domains within TF complexes, they are robust to variability in TFBS spacing. An alternative model to explain conserved patterns of TFBS co-occurrence without rigid spacing relates to a recently proposed “dependency grammar” (Jindal and Farley 2021). This considers that multiple parameters including TFBS affinity, order, spacing, and orientation are tuned by selective pressures to yield viable patterns and levels of expression from key enhancers. Thus, interplay between the spacing and affinity of sites dictates the functional capabilities of enhancers without either parameter needing to be rigidly set (Farley et al. 2015, 2016). This may be linked to suboptimization of enhancers via hinderance of combined TF binding, preventing overactivation (Farley et al. 2015). This is especially important for trans-dev genes, in which their overactivation or ectopic expression may lead to harmful alterations in cell fate assignment during development. Enhancer configuration with genuine but less rigid grammatical logic may also offer an evolutionary advantage because such enhancers are more likely to be robust to minor sequence changes without evident alterations in function.

Given our results suggest combinatorial rules in terms of TFBS co-occurrence without rigid spacing in PFG HCNE CRMs, this supports the dependency grammar model in which coregulated genes are regulated by similar sets of TFs, but the tuning of grammatical parameters varies to ensure correct functional outputs. This tuning via dependency grammar is likely to be critical, hence the high degree of conservation of the HCNEs.

Overall our analyses reveal a compendium of HCNE CRMs acting across the 400 million years of evolution separating zebrafish and humans. These HCNE CRMs should be prioritized for analysis of commonalities of gene regulatory control programs controlling neural and pancreatic development. Our analyses reveal consistent patterns of TFBS co-occurrence in endoderm HCNE CRMs that should form the basis of future studies of TF combinatorial and competitive binding to reveal conserved mechanisms tuning the correct spatiotemporal expression of trans-dev genes.

Zebrafish strains and transgenics

ABix, Tg(-5.0sox17:GFP)^ha01Tg (Mizoguchi et al. 2008), Tg(kdrl-HsHRAS:mCherry)^s896 (Chi et al. 2008a), and Tg(gata1a:DsRed)^sd2Tg (Traver et al. 2003) fish were reared as previously described (Westerfield 2007). All zebrafish studies fully complied with the United Kingdom Animals (Scientific Procedures) Act of 1986 as implemented by University of Warwick.

Preparation of cells for RNA-seq and ATAC-seq

Embryos were dechorionated using Pronase and dissociated in a collagenase:trypsin blend as described in the Supplemental Materials. Cell sorting was performed using a Becton Dickinson FACSAria fusion cell sorter with a 100 µm nozzle and sheath fluid pressure of 25 pounds per square inch (psi). Debris was eliminated based on side scatter area (SSC-A) and forward scatter area (FSC-A), and singlet cells were selected based on FSC-A and forward scatter width (FSC-W). Next, the single cells were sorted based on fluorescence signal from green fluorescent protein EGFP (B488-530) and red fluorescent proteins mCherry and DsRed (YG561-610). Flow cytometry was operated and analysed using BD FACSDiva software.

RNA-seq

Total RNA was prepared from sorted cells using a Qiagen RNeasy mini kit with on-column DNase I treatment (Qiagen) following the manufacturer's protocol. Ten nanograms of total RNA per sample was used to construct sequencing libraries using the NEBNext Ultra II directional RNA library prep kit for Illumina according to the manufacturer's instructions and was sequenced using a NovaSeq 6000-S4-type flow cell. Reads were mapped onto zebrafish genome, danRer11 using STAR (Dobin et al. 2013). HTSeq was used to count aligned reads per gene (Anders et al. 2015). The raw count matrix was imported into iDEP (Ge et al. 2018), and differentially expressed genes (DEGs) were identified using default parameters. Heatmaps of DEGs were produced using Morpheus (https://software.broadinstitute.org/morpheus).

ATAC-seq

Fifty thousand sorted cells per sample were used for Omni-ATAC-seq using adapted methods outlined by Buenrostro et al. (2015) and Corces et al. (2017). Full details are provided in the Supplemental Materials. Libraries were sequenced by the genomics facility at the University of Warwick with Illumina NextSeq 500 using the high output kit v2.5 (FC-404-2002) and by Novogene using a NovaSeq 6000-S4-type flow cell. Each sample used different barcoded reverse primers, allowing for the samples to be multiplexed for sequencing.

Reads were trimmed from the 3′ end to a uniform 75 bp using seqtk so that mapping of reads was not impacted by sequencing read length. Adapter sequences were trimmed using Trimmomatic version 0.39 (Bolger et al. 2014). Data quality control and confirmation of the adapter removal were performed using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were aligned to the genome reference consortium zebrafish build 11 (GRCz11v97) using Bowtie 2 (Langmead and Salzberg 2012) using “very sensitive” and a 2 kb maximum fragment length for paired-end alignments. SAMtools (Li et al. 2009) was also used to filter out reads that mapped to the mitochondrial genome, duplicated reads, and low-quality reads with a MAPQ score < 22. The numbers of reads per sample used in subsequent analyses are shown in Supplemental Table 16.

Accessible regions (peaks) per sample were called using MACS2 (Zhang et al. 2008) with default parameters except broad peak setting, mappable genome size of 1.4 × 10⁹, false-discovery rate cutoff of 0.05, bandwidth of 300 bp, and a high-confidence fold-enrichment between five and 50. To control for Tn5 sequence-specific signature bias, peaks were called against an ATAC-seq control library produced from purified genomic DNA (Buenrostro et al. 2013).

DARs between the cell populations and time points were identified in R version 4.2.1 (R Core Team 2021) using DiffBind version 3.0.15 (Buenrostro et al. 2013) with the DESeq2 option (Love et al. 2014) on the broad peaks files (.broadPeak) produced by MACS2 without merging of biological replicates. Summits were set to false so that the summit heights (read pileup) and locations were not calculated for each peak, allowing the whole peak size to be considered. DARs identified with FDR ≤ 0.01 DARs were used for subsequent analysis. DARs were assigned to the nearest TSS, and distributions of DARs relative to genes were assessed using ChIPseeker 1.32.1 (Yu et al. 2015) and GRCz11 (danRer11) Ensembl gene annotation version 97.

GO term enrichment for biological processes, molecular functions, and pathways was performed using overrepresentation analysis (ORA) in WebGestalt 2019 (Liao et al. 2019). FishEnrichr (Chen et al. 2013; Kuleshov et al. 2016) was used to identify zebrafish anatomical structures associates with genes lists.

Integrative Genomics Viewer (IGV) version 2.5.0 (Robinson et al. 2011) was used to visualize mapped ATAC-seq and ChIP-seq reads, DARs, and HCNEs. ATAC-seq peaks from biological replicates were merged using SAMtools (Li et al. 2009) to show representative peaks. Tracks were normalized using the CPM function in deepTools bamCoverage (Ramírez et al. 2016).

To produce heatmaps of ATAC-seq data, BAM files were downsampled to equalize read numbers between conditions using Picard. Heatmaps were produced using seqMINER v1.3.4 using the KMeans enrichment linear clustering normalization method (Ye et al. 2011).

ChIP-seq data analysis

Previously published H3K27ac, FOXA2, GATA4, GATA6, HNF1B, PDX1, and PTF1A ChIP-seq data were downloaded from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE52658, GSE183672, and GSE114102 (Loh et al. 2014; Lee et al. 2019; Miguel-Escalada et al. 2022). Reads were aligned to hg38 using Bowtie 2 (Langmead and Salzberg 2012) with default parameters apart from -N 1. ChIP-seq peaks were called using MACS2 (Zhang et al. 2008) using default settings against cognate controls. A more stringent Q-value of <10⁻⁸, instead of the default <0.05, was qualitatively determined as an appropriate cutoff for the called peaks by looking at the peaks around key genes known to be expressed in endodermal cell populations, as well as genes expressed specifically in nonendodermal cell populations.

Analysis of HCNEs overlapping ATAC-seq and ChIP-seq data sets

HCNEs comparing human (hg38) and zebrafish (danRer10) genome builds using a window size of 30 bp and 70%–100% identity threshold were downloaded from ANCORA (Engström et al. 2008; http://ancora.genereg.net). danRer10/hg38 HCNEs were converted to danRer11/hg38 using liftOver (Hinrichs et al. 2006). BEDTools intersect (Quinlan and Hall 2010) was used to determine if danRer11/hg38 HCNEs overlap with zebrafish danRer11 DARs from ATAC-seq and human hg38 ChIP-seq peaks.

Functional and anatomical enrichment analyses were performed for hg38 HCNE sets using GREAT version 4.0.4 with the default settings (McLean et al. 2010; Tanigawa et al. 2022). Dual analyses were performed with background regions set to the whole genome or a list of all HCNEs. Heatmaps of the top 20 enriched terms per set of HCNEs were selected based on hypergeometric false discovery rate corrected q-values (HyperFdrQ).

Motif enrichment analysis

De novo and known motif enrichment analysis was performed using HOMER v4.11.1 (Heinz et al. 2010). For known motif enrichment analysis, the top 30 most significant motifs in each HCNE set (P < 1 × 10⁻⁸) were considered, leading to a combined set of 60 motifs. A heatmap indicating enrichment of these motifs in each HCNE set was produced and clustered using Morpheus (https://software.broadinstitute.org/morpheus). To annotate expression of TFs corresponding to these motif's RPKM values from relevant human endoderm, RNA-seq data sets were downloaded from GEO under accession numbers GSE52658 (AFG, PFG, and MHG endoderm) (Loh et al. 2014) and GSE216266 (passage zero PPs) (Jarc et al. 2024). Co-occurrence of selected motifs within HCNEs were analyzed using the paired motif enrichment tool (PMET; implemented at https://pmet.online and https://github.com/duocang/PMET-Shiny-App) (Rich-Griffin et al. 2020) on genomic intervals using default parameters. PMET output was used to produce a matrix using the acast function in reshape2 in R prior to producing a heatmap using conditional formatting in Microsoft Excel. Pairwise analysis of spacing of motifs from HOMER was performed using SpaMo in the MEME suite (Whitington et al. 2011; Bailey et al. 2015). All position frequency matrices from significant de novo and known motif enrichment analysis using HOMER were provided as primary and secondary motifs for pairwise comparison of spacing.

Cloning of reporter constructs for reporter assays

Putative enhancer and HCNE elements of interest were PCR-amplified from zebrafish or human genomic DNA using custom primers containing attB4 and attB1 sequences and Q5 high-fidelity DNA polymerase (NEB) and cloned by Gateway recombination into pDONR-P4-P1R (Thermo Fisher Scientific). All primers used are listed in Supplemental Table 17. To allow visualization of reporter activity in Tg(sox17:GFP) zebrafish, constructs were generated to express mCherry downstream from the putative enhancer. To do this, published plasmids pENTRbasEGFP (Addgene 22453), pENTREGFP2 (Addgene 22450), and pDESTtol2pACrymCherry (Addgene 64023) (Villefranc et al. 2007; Berger and Currie 2013) were modified using NEBuilder HiFi DNA assembly master mix (NEB) to switch EGFP and mCherry between the plasmids. To use the E1b promoter as the basal promoter in reporter assays, it was amplified from gata2a-i4-E1b-GFP-Tol2 (kindly gifted by Dr Rui Monteiro) and inserted into pENTR mCherry by HiFi cloning. Following sequence verification, hnf1ba putative enhancer/HCNE, pENTR E1bP:mCherry, p3E-mcs1 (Addgene 49004) (Moore et al. 2013), and pDESTtol2pACryEGFP constructs were recombined by Gateway recombination (Thermo Fisher Scientific). This generated plasmids with putative hnf1ba enhancer elements upstream of the E1b promoter and mCherry, as well as the cryaa promoter upstream of EGFP.

Embryo compound treatments, reporter assays, and imaging

Embryos for imaging were incubated from 22 hpf with 0.003% 1-phenyl 2-thiourea (PTU) to block pigmentation and improve optical transparency. Embryos were immobilized for imaging using 30 µg/mL tricaine (ethyl 3-aminobenzoate methanesulfonate E10521-10G Merck) and mounted and orientated in 1% agarose molds created using a 3D-printed stamp (Kleinhans and Lecaudey 2019). Embryos were imaged and images processed as described in the Supplemental Materials.

To analyze the effects of RA signaling, embryos were treated from 5 hpf with either RA (Sigma-Aldrich R2625), DEAB (Sigma-Aldrich D86256), BMS493 (APExBIO B7415), or vehicle (DMSO 0.1%) as a control. Following treatment, embryos were incubated in the dark to prevent degradation of the light-sensitive chemicals.

For reporter assays, Tg(sox17:GFP) embryos were microinjected at the one cell stage with 25 pg reporter constructs combined with 25 pg Tol2 transposase-capped mRNA produced from pCS2-Tol2. Embryos were screened using widefield fluorescent microscopes for expression of cry:EGFP in the lens at 48 hpf for confirmation of reporter construct integration into the genome. The location of mCherry expression in these embryos was observed using widefield fluorescent microscopes and scored. Embryos of interest were imaged further using confocal microscopes.