Construction of an atlas of transcription factor binding during mouse development identifies popular regulatory regions

Transcriptional regulation during embryonic development relies on the tissue-specific deployment of combinatorial systems of DNA-binding transcription factors. Transcription factors physically associate with regulatory regions along the genome, recruiting non-DNA-binding co-factors and chromatin remodelers. Their combined action generates chromatin domains that are permissive for the activity of RNA Polymerase II, leading to transcription (Klemm et al., 2019). Characterizing the genome-wide binding profile of all the different classes of proteins that are involved in this process, hereafter referred to as gene regulators (GRs), is therefore a necessary step for understanding how the genome is regulated in each cell during development (Merika and Thanos, 2001). As technologies that map the genome-wide physical behavior of these proteins are not trivial to employ, classic experimental approaches typically aim at targeting one or a limited number of such players. Here, we propose that, to truly understand how the genome is differentially regulated across cell types during development, one should measure the activity of all gene regulators at different developmental time points. Vertebrate genomes, however, encode for ∼1600 transcription factors and a comparable number of other regulators (Ng et al., 2021). The number and combinatorial activity of GRs render their comprehensive study difficult experimentally. State-of-the-art technologies now allow us to approach this endeavor.

Among the technologies that have been developed to profile the DNA-binding behavior of GRs, cleavage under targets and release using nuclease (CUT&RUN; hereafter C&R) (Skene and Henikoff, 2017; Skene et al., 2018) and CUT&Tagmentation (CUT&Tag; hereafter C&T) are emerging as methods of choice. C&R, in contrast to chromatin immunoprecipitation assays with sequencing (ChIP-seq), is conducted under native conditions without cross-linking agents, which drastically reduces the generation of experimental artifacts (Landt et al., 2012). Moreover, C&R and C&T permit scaling down in cell numbers, with some applications even yielding results from single cells (Bartosovic et al., 2021). Although both C&R and C&T reliably work for the detection of abundant and stable post-translational modifications of histones, they sometimes ‘suffer’ when targeting less stable chromatin regulators such as lowly expressed transcription factors and non-DNA-binding GRs (Meers et al., 2019). This drawback led us to develop C&R-low volume-urea (LoV-U), which robustly detects all types of GRs, and facilitates intermediate-scale parallel profiling thanks to low sample volumes (Zambanini et al., 2022).

Here, by using C&R-LoV-U (generally referred to as C&R), we initiated the systematic mapping of GRs during mouse organogenesis. We produced a dataset of >100 individual C&R tracks of GRs acting downstream of developmental signaling pathways in four different mouse tissues (forelimbs, hindlimbs, branchial arches and liver), at two stages of organogenesis (10.5 and 11.5 days post-coitum). While we realize that this is only the beginning of a larger enterprise, this core dataset allowed us to observe key regulatory mechanisms of development, and gives us the impetus and basis from which to scale-up. Our findings include measuring the extent to which each GR has shared or tissue-specific activity, how GR deployment changes over time, and the convergence – or crosstalk – of developmental signaling pathways on the genome.

An overview is needed to understand how GRs coordinate the complexity of gene expression during development. To initiate the simultaneous profiling of the genome-wide binding behavior of GRs in developing mouse tissues, we decided to perform C&R on tissues dissected from embryos at embryonic day (E) 11.5, when organogenesis is established and ongoing (Cao et al., 2019). We selected four tissues for analysis: forelimb buds (FL), hindlimb buds (HL), liver (L), a prominent hematopoietic site at this stage (Jagannathan-Bogdan and Zon, 2013), and branchial arches (BA), through which neural crest cells migrate to contribute to facial and trunk processes (Hari et al., 2012) (Fig. 1A).

Crucial for achieving experimental success in C&R is the functionality of the antibody employed. We tested a collection of antibodies (often multiple per target) against >50 known developmental GRs. On average, one in four antibodies yielded reliable genome-wide signal, underlying the intrinsic difficulty of detecting lowly abundant nuclear proteins in their native conditions by C&R using commercially available antibodies. Here, we present a curated set of genome-wide binding profiles (made up of a merge of two biological replicate tracks) of GRs that surpassed our experimental and computational quality controls across the four selected tissues (Fig. 1B; Figs S1, S2). The profiled targets include transcription factors and co-factors chosen to monitor the activity of key developmental pathways [ETV4 tracking Fgf signaling (Mao et al., 2009); STAT5A tracking Jak/Stat signaling (Smith et al., 2023); DLX5 for Bmp signaling (Holleville et al., 2003); SMAD4 for Tgfβ signaling (Guglielmi et al., 2021); GLI3 for Shh signaling (Büscher et al., 1997); TCF7 and LEF1 for Wnt signaling (Valenta et al., 2012)]; chromatin regulators such as the histone deacetylase enzyme HDAC1 (Brunmeir et al., 2009) and the SWI/SNF component SMARCA4 (also known as BRG1) (Hota et al., 2022); and post-translational modifications on histones associated with active transcription such as H3K4me3, which signals active promoters (Schneider et al., 2004), H3K4me1 for active and poised enhancers (Bae and Lesch, 2020) and H3K27ac, which marks promoters, enhancers and super-enhancers (Calo and Wysocka, 2013).

Producing these numerous binding profiles allowed us to assess (1) the tissue-specific genomic activity of developmental signaling pathways along with the resident chromatin dynamics, (2) time specificity of gene regulator binding and (3) the existence of convergent regulation or crosstalk of pathways (Fig. 1A). The analyses to measure (1), (2) and (3) represent the bulk of this study.

Spearman's correlation analyses followed by hierarchical clustering across all the datasets immediately revealed interesting patterns (Fig. 1B; binding profiles are obtained by merged replicates; single replicate tracks are shown in Fig. S1A). For example, we found groups of binding profiles that clustered based on one individual target, indicating instances in which that target possessed a relatively similar binding activity across tissues, as for H3K27ac, H3K4me1 and – perhaps less expected – SMAD4 and DLX5 (Fig. 1B, black boxes). In other cases, clustering was observed within tissues, where several factors appeared to co-regulate the same groups of genes in the liver (Fig. 1B, red boxes). Notably, these groups were also marked by H3K27ac, likely reflecting the permissive chromatin state of the underlying loci. IgG negative controls clustered aside from all biological targets (Fig. 1B, gray box), conferring overall credibility to the identified peaks. Consistent with this, principal component analysis (PCA) separated the samples mostly by target along PC1, which explained ∼58.7% of the variance with histone marks farthest to the left and transcription factors towards the right (Fig. 1B, right; ellipses shaded by factor color code; see Fig. S1B for individual replicate track PCA). Both patterns – overlap of GR activity within each tissue or similarity of GR binding profile across tissues – were visible upon close inspection of our tracks across exemplifying genomic loci: while several loci presented a relatively homogeneous signal across tissues (Fig. 1C, black boxes on the left), many others displayed a high degree of tissue specificity. These included binding sites within the locus of the neural crest regulator Pitx2 mostly and expectedly occurring in BA, the tissue containing cells from the neural crest (Kioussi et al., 2002) (Fig. 1C, center left), the HoxA gene cluster being prominently regulated in limbs (Nemec et al., 2017) (Fig. 1C, center right), and the hematopoietic master gene Scl (Tal1) in L, the embryonic hematopoietic niche (Real et al., 2012) (Fig. 1C, far right). Individual replicate tracks at additional loci are shown in Fig. S1C. The expected tissue-specific regulation of representative loci assured us of the reliability of our experimental setup and consequent analytical output.

Most of our binding profiles represent novel datasets: aside from GLI3 in the FL, our profiles are the first C&R of these particular targets in vivo during development, and with exception of the histone modifications and a few targets in the FL, we could also not find any existing ChIP-seq data. Regardless, we compared our datasets to the closest matching published data we could identify for histone marks (Fig. S3) and other GRs (Fig. S4) (Arase et al., 2017; Attanasio et al., 2014; Gamart et al., 2021; Gorkin et al., 2020; Guo et al., 2021; Hu et al., 2023 preprint; Lee et al., 2022; Lex et al., 2020). Consistent with previous literature (Hu et al., 2023 preprint), C&R generally identified fewer peaks than the corresponding ChIP-seq data, but with a high degree of overlap, depending on the degree to which we were able to match the tissue and time point. Different technologies and peak-calling strategies might explain these differences, and typically ChIP-seq leads to more artifacts (Boyd et al., 2021; Hu et al., 2023 preprint). However, it is worthwhile considering that, thanks to using C&R, we started each experiment from the equivalent material extracted from three limbs per sample, which is roughly 1/100 of the material required in a typical ChIP-seq (Zambanini et al., 2022; Zimmerli et al., 2020), granting our approach increased scalability, decreased material costs, and enormously improved ethical considerations.

GRs have pleiotropic activity in different tissues and their compound interactions are central to development (Andrey et al., 2017). We therefore decided to measure the differences and overlaps of the genome-wide binding across tissues for each targeted factor, a previously under-investigated phenomenon during development, now made possible with our datasets. We calculated how many peaks were identified in only one, two or three tissues or across all four tissues (overlapping by at least 1 bp), and plotted the resulting ‘decay’ curves (Fig. 2A). Some targets, such as H3K4me3 and H3K4me1, presented a relatively low slope with a great fraction of peaks still shared among two, three and four tissues (Fig. 2A, top left). This indicates a large overlap of active regulatory regions in the different developmental contexts. Most other targets had steeper slopes – the number of peaks was highest in a single tissue – indicating higher tissue-specific actions (Fig. 2A). For downstream analyses, we focused on one example of high overlap/low decay (H3K4me3), one with an intermediate slope (HDAC1) and one presenting a steep slope (GLI3) (Fig. 2B). All displayed peaks present across all the tissues, and varying fractions of tissue-specific peaks. For H3K4me3, these peaks included known organ-specific genes (Fig. 2B).

The number of shared peaks in the HDAC1 dataset was proportionally smaller than that observed for H3K4me3 and H3K4me1, and the examples of tissue-specific regulation were more numerous (Fig. 2C, left; signal intensity heatmaps of shared and tissue-specific peaks, with HL showing the most tissue-specific peaks). As expected, sequences underlying unique HDAC1 peaks were enriched in signatures of tissue-specific factors, likely explaining the tissue specificity, including, for example, GATA and SCL in the L (Fig. 2C, right). Interestingly, FL, HL and BA all presented motifs for PITX1 binding in non-overlapping peaks (Fig. 2C, right). PITX1 has been shown to be a HL-specific factor, acting on a developmental program shared between FL and HL (Nemec et al., 2017). Our data suggest that PITX1 is not limited to the core regulation of limb formation, but is active in other tissues, such as the BA, where it binds to different regulatory regions (again, the peaks presenting PITX1 motif are non-shared across tissues). GLI3 activity seemed to be predominantly tissue-specific with only few instances of peaks shared by at least two tissues (Fig. 2D, left; signal intensity heatmaps of shared and specific peaks). As GLI3 is a DNA-binding protein, we could confirm the presence of its cognate motif in the DNA sequence underlying the identified peaks (Fig. 2D, right). GLI3 motifs were highly enriched in BA and FL, emphasizing its function in these tissues (Akiyama et al., 2015; Elliott et al., 2020). Interestingly, while GLI3 drives digit number specification both in FL and HL (Lopez-Rios et al., 2012), our data suggest that it might do so by different mechanisms. The absence of GLI3 motif enrichment in HL and L points to its ability to cooperate with other DNA-binding transcription factors as a co-factor (Fig. 2D, right), a bimodal modality of action similar to that which we previously observed for TBX3 (Amaia et al., 2023 preprint). These differential modes of action of GLI3 could provide a new rationale for a mechanistic understanding of the multitude of human phenotypes caused by GLI3 variants (Démurger et al., 2015).

Gene regulatory networks are often dynamic and change along developmental trajectories (Tsankov et al., 2015; Zuniga and Zeller, 2020). We investigated the dynamics of GR binding by comparing two developmental time points in the liver – the site of embryonic hematopoiesis. We mapped the genome-wide binding profile of all the targets as above at the additional earlier time point of E10.5, then compared it with the dataset obtained at E11.5 (Fig. 3A). Comparison of E10.5 and E11.5 identified that transcription factors and co-factors tend to cluster within time points (Fig. 3B, left, red boxes), whereas chromatin marks are grouped depending on the type of modification, regardless of time (Fig. 3B, left, black boxes). We considered this interesting: the common localization of GRs downstream of diverse signaling pathways at a specific developmental stage might reflect their convergence onto broad genomic regulatory regions such as, as previously reported, super-enhancers (Hnisz et al., 2015). In the PCA plot, this pattern was less clear, although PC2 (representing 10.9% of the variance) might reflect the time points as all targets moved down on this axis from E10.5 to E11.5 (Fig. 3B, right). Our data indicate that these regulatory regions co-regulated by GRs are different depending on the developmental time. Post-translational modifications of histones, by contrast, seemed more stable over time, possibly indicating the consolidation of genomic identity that underlies the acquisition of a defined cell fate. Signal intensity profiles across all peaks supported this interpretation by showing how well intensity is conserved between E10.5 and E11.5 in the case of H3K4me3 and H3K4me1, whereas it varies for other factors, in particular GLI3 and STAT5 (Fig. 3C). Fluctuating signal over genomic regions, we reasoned, might reflect how dynamically these factors ‘jump’ to new loci for the regulation of ontogenetically subsequent gene expression programs.

The E11.5 datasets for LEF1 and SMAD4 were not included in the previous analysis as they failed our quality control metrics, and we never obtained successful DLX5 binding profiles in L at either time point. Multiple experimental attempts to obtain binding patterns of LEF1 and SMAD4 at E11.5 failed, despite the fact that they were easily obtained at E10.5, and other targets worked well at E11.5 in parallel. These targets, in published RNA-seq (Cardoso-Moreira et al., 2019) exhibit an increase in expression at E10.5 followed by a drop in expression at E12.5 (Fig. 3D). It is therefore plausible that these GRs are detectable only at an earlier time point. Examples of more detailed analyses on the time-specific dynamics of our targets to exemplify the utility of our datasets are included in Fig. S4.

Measuring how cell communication signals are integrated on the genome is fundamental to understanding cell differentiation during organ formation (Reményi et al., 2004). Our setup, which captures multiple transcriptional effectors of extracellular signals, allowed us to determine the interactions of these signaling pathways at the chromatin level. Here, we focused on the FL, as a great wealth of data concerning the underlying gene regulatory networks exists for this tissue (Andrey et al., 2017; Zuniga and Zeller, 2020). Spearman's correlation of all mapped factors at E11.5 in FL unearthed several interesting observations. For example, we considered notable the proximity, within the clustering, between the chromatin remodelers SMARCA4, HDAC1 and H3K27ac-marked regions (Fig. 4A). SMARCA4 is part of the SWI/SNF complex, which hydrolyzes ATP to remodel chromatin through nucleosome sliding and eviction, and it is typically associated with rendering the chromatin permissive for transcription (Centore et al., 2020). By contrast, HDAC1, by removing acetyl groups deposited by CBP/p300 on H3K27 residues, represses transcription and prevents premature expression of developmental genes (Dsilva et al., 2023; Wang et al., 2022). Hence, the colocalization of HDAC1 and SMARCA4 on H3K27ac-positive regions was surprising. However, others have found that HDAC1 can be recruited within super-enhancers, which display H3K27ac, to activate transcription (Lai et al., 2023), or that the balance of p300 and HDAC1 activities controls nucleosome eviction by the SMARCA4-containing SWI/SNF complex in macrophages (Pietrzak et al., 2019). Moreover, our recent work in the context of Wnt signaling implied that HDAC inhibition was responsible for the inability to open chromatin for transcription (Pagella et al., 2023), which could be interpreted as a function of HDAC enzymes in gene activation rather than repression. The tight association between these factors in FL cells might indicate that, in this context, SMARCA4 and HDAC1 interact to keep the chromatin open. Accordingly, most of the overlapping peaks between SMARCA4 and HDAC1 are either on H3K27ac- or H3K4me1/me3-positive regions, which are both associated with active and transcription-permissive chromatin (Fig. 4B). Concerning HDAC1, although its role in limb organogenesis has not been addressed, our data indicate a requirement of this factor for correct FL morphogenesis on a molecular level, complementing other phenotypic reports (Paradis and Hales, 2015). Our dataset also identifies idiosyncratic binding behaviors whereby some genomic regions are mostly occupied by individual factors, such as SMARCA4 and STAT5 (Fig. S4).

One limitation of our bulk approach is that it does not allow instances of locus co-regulation by GRs in the same cell to be distinguished from different factors being active on the same regulatory regions but in different cells. To disentangle this, we used existing single-cell RNA-sequencing (scRNAseq) data (He et al., 2020) and plotted the cell type-specific expression of each of our C&R targets (Fig. 4D). While some genes, such as Hdac1 and Smarca4, had a broad expression pattern encompassing many cell types, certain targets showed more of a tissue-specific distribution. For example, Tcf7 and Lef1 expression was highest in the different mesenchyme and fibroblast clusters, respectively, indicating that our TCF7/LEF1-specific C&R reads should be traced to these cell types (Fig. 4D). Stat5a was lowly expressed overall but had detectable expression in macrophages, consistent with its known role in the immune system, as well as in chondrocytes and specific types of muscle cells (Fig. 4D). Therefore, we conclude that STAT5A regulates the identified targets in this cell type. This simple data integration shows that scRNAseq datasets, and by extension also spatial transcriptomics results, can be proficiently used to extrapolate from which cell population the C&R signal derives.

Our FL tracks combined identified 36,717 peaks called in at least one dataset – FL-specific regulatory regions – many of which (∼18,000) are bound by only one of the factors that we tested (Fig. 5A). We generated lists of peaks that were bound by more than one GR and compared them with those bound by at least two, three, four – and so on – factors and noticed that the number of regulatory regions contacted by an increasing number of GRs drops rapidly (Fig. 5A). This possibly underscores that combinations of GRs and DNA regulatory elements are the most common means to drive the expression of differential gene expression programs (as opposed to, for example, having different GRs binding on the same regulatory element or the same GR binding on different regulatory elements). We were surprised, however, to note that the groups of peaks presenting a higher number of bound GRs, which we colloquially refer to as ‘popular’ regulatory regions, were characterized by increased evolutionary conservation as assessed with phastCons, which measures the degree of negative selection of a region on a scale of 0 (no conservation) to 1 (full conservation) (Fig. 5B). The trend became statistically significant when we divided all regulatory regions into four groups (1-2, called in one or two datasets; 3-5, called in three to five datasets; 6-8, called in six to eight datasets; 9-12, called in 9 to 12 datasets), indicating that the regions called as peaks in a greater number of datasets are evolutionarily conserved to a higher degree and thus likely essential for gene regulation (Fig. 5C). We focused on the group with the highest statistical score, which included 99 peaks presenting the concurrent presence of at least 75% (9/12) of our targets (H3K4me1, H3K4me3, H3K27ac, HDAC1, SMARCA4, LEF1, TCF7, GLI3, SMAD4, DLX5, STAT5, ETV4). All these 99 genomic regions presented a high signal-to-noise ratio for all factors, while signal was absent in the IgG control samples (Fig. S5A). The 99 regions were located overwhelmingly in close proximity to promoters, with over 80% of region–gene associations being found within 5 kb of the transcription start site (Fig. 5D). According to ChromHMM data in the E11.5 limb (Gorkin et al., 2020), popular regions overlap with 147 annotated regions, mostly in promoters or proximal regions (53% promoters, 14% promoter-proximal enhancers, and the rest no significant signals or multiple annotations). This group of genomic elements, we hypothesized, might represent an ensemble of important hubs where several pathways converge to regulate the expression of essential genes for limb patterning. Consistent with this, motif analysis revealed enrichment of known determinants of limb morphogenesis (Fig. 5E), such as the LIM homeobox factor LHX9, which has been shown to integrate signaling pathways to control the three-dimensional patterning of growing vertebrate limbs (Tzchori et al., 2009), and engrailed 1 (EN1), one of the initial inducers of the apical ectodermal ridge (Loomis et al., 1996).

Surprisingly, we did not detect enrichment of motifs relating to our target factors within these popular regions. Since CUT&RUN is capable of detecting both the bound regions and often regions in proximity of the 3D space, we hypothesized that the different factors might indeed be binding to different enhancers in different cell types, which all might converge onto promoter regions to regulate these genes. To test this possibility, we integrated HiChIP data targeting H3K4me3 in mouse limbs at E12.5 (Yu et al., 2024) and found that 66 of the 99 popular regions were associated to DNA loops. Since the time point is imperfectly matched, it is possible that some of the remaining 33 popular regions would display looping if examined at E11.5. These 66 popular regions were connected to a total of 371 loops. In comparison, 99 randomly selected regions from H3K4me3 peaks, shuffled and randomly selected five times, were connected to an average of only 285 loops, which is ∼25% less. This suggests the popular regions may be more highly interconnected with other regulatory regions, and multiple enhancers might partially explain the apparent high occupancy of these regions.

The relevance of these regulatory regions seemed to go beyond limb patterning: STRING showed a high interconnection between the 99 hub regions and the targets themselves (Fig. 5F), and GREAT annotated these 99 regulatory regions to genes that, when mutated, cause early embryonic lethality (Fig. 5F, gray bubbles; obs/exp 3.24, FDR 0.028). Among these are the ubiquitously functional DNA repair protein BRCA1 (Hakem et al., 1998) and the histone methyltransferase PRMT1, absence of which leads to stalled development in vertebrate model systems (Shibata et al., 2020) (Fig. 5G, example tracks). We examined how highly bound these regions were across tissues, and found that also in the HL, BA and L, these regions showed high signal in many targets (Fig. S5B). Interestingly, when we plotted Spearman's correlation of all the datasets within the hub regions, we found that target-specific idiosyncrasies allowed for clustering that was based mainly on the target protein (Fig. S5C).

To explore the possible connotations and functions of these regions within the human genome, we converted the coordinates from mm10 to hg38 (minimum reMap ratio 0.5). Of the 99 regions, 88 successfully mapped onto hg38; in comparison, of 99 regions randomly selected from all of the FL peaks, only 66 successfully converted. In humans, these regions generally belong to CpG islands and coincide with ENCODE candidate cis-regulatory elements, both marked as enhancers and promoters, although they show little overlap with conversion of known VISTA enhancers (Fig. S5D). A high level of ChIP-seq reMap signal in most regions indicated that, like in our mouse datasets, these are ‘popular’ regions bound by a multitude of GRs (Fig. S5D). JARVIS scores of negative selection in humans were also high for these regions, suggesting that they retain their functional relevance for human development (Fig. S5D).

The central question of developmental biology is how the cells of an organism, despite having inherited the same genome, can build an embryo by activating the right genes at the right time to ensure the desired differentiation programs. Decades of research in molecular biology have clarified that this problem could be solved by determining the activity of transcription factors, co-factors, and other GRs, which bind the DNA at gene regulatory regions and contribute to gene activation or repression (Dickel et al., 2013; Dunham et al., 2012). Several technologies, including ChIP-seq and now the more proficient C&R, allow determination of where a given GR binds across the genome. However, ChIP-seq, C&R, and all the related technologies are limited to detecting one target at the time, as they rely on the use of protein-specific antibodies. This limitation has caused researchers to focus on one, or on a limited numbers of, favorite GR candidates. However, GRs do not act alone. Rather, they are needed in a combinatorial fashion to functionally engage regulatory regions and modulate gene transcription (Reményi et al., 2004). Hence, the study of gene expression will greatly benefit from the parallel measurement of all GRs acting in the context of interest.

Here, we initiated the assessment of the genome-wide binding profiles of multiple GRs in selected developmental contexts as an attempt to show the feasibility of a larger-scale venture. Certainly, the sheer number of GRs encoded in the vertebrate genome – ∼1600 transcription factors and likely a comparable number of co-factors and chromatin remodelers – might be discouraging and render this endeavor unattainable. Big consortia, such as ENCODE, successfully approached this enormous task by mapping the genome-wide binding site of many GRs using ChIP-seq (Landt et al., 2012; Luo et al., 2020). However, most of this work is conducted in human cell lines (Partridge et al., 2020), and it is not always clear to what extent they are representative of the tissue or disease of origin (Ben-David et al., 2018; Lopes-Ramos et al., 2017). This is particularly relevant for studies of physiological mechanisms occurring during normal embryonic development. With the current study, we started and support the next step: the monitoring of all GRs in relevant in vivo models of developmental processes. The resulting information concerning the combinatorial action of GRs could be used to reveal how the genome is deployed in different organs during development to generate cell type-specific gene expression programs.

We advocate for the gradual construction of an atlas describing the action of virtually every GR acting during different moments of organogenesis. Our C&R-LoV-U protocol now allows for parallel processing of up to ∼60 samples per experimental run, maintaining reasonable time frames and sample integrity with only one scientist at the bench for less than 2 days of work. We plan to use cultured cells that are not limiting in number, such as MEF cells, to test antibodies against GRs, and then once suitable ones are identified, turn back to in vivo experiments. This study used approximately 100,000 cells per sample, meaning that in the future we could profile, for example, 15 targets in four tissues from a single litter of E11.5 embryos. This will allow us, and others, to scale-up this initiative while maintaining ethical practices with regard to animal usage and without skyrocketing experimental costs. This initiative, which we refer to as transcription factor binding-omics (TFBomics), will constitute a colossal repository whereby researchers interested in gene regulation during development, can access, download, or simply consult the data instead of conducting sophisticated experiments such as ChIP-seq and C&R. Our hope is that other interested researchers in the field of genome regulation during development will assist us in gathering further resources and experimental power to profile and characterize all GRs.

Despite its limiting size, the current study is already enough to allow analytical possibilities that permit new discoveries, among the many that this atlas could offer. Here, we will briefly mention three findings that we consider the most relevant. Our first discovery derives from the exploration of how the activity of select GRs varies across tissues at the same developmental time point (Fig. 2). Several GRs possess ample tissue-specific activity – that is, peaks identified in one tissue and not in the others – but also common actions, whereby they appear to regulate the same DNA elements in different tissues. This provides a clear-cut molecular rationale for gene pleiotropy. It is interesting to note that, when the binding signal of a given GR is compared between tissues, we observed instances of variability in its intensity (see, for example, HDAC1 in Fig. 2C) indicating that GRs are present but bind with lower affinity or in varying fractions of the cell population depending on the tissue, as well as more clear-cut cases (see, for example, the HL-only peaks of GLI3 in Fig. 2D) in which sets of tissue-specific targets seem to be exclusive of a particular developmental trajectory. Our second discovery is the investigation of how GRs act over time. Domcke and Shendure recently advocated for the construction of cell trees rather than atlases, as these will better explain functional gene expression during ontogeny (Domcke and Shendure, 2023). We profoundly agree, and have the ambition to make the TFBomics dataset grow and include as many organ-specific differentiation time points as possible. Two stages of liver development, for example, were enough to discover that the two time points (E10.5 and E11.5), despite being relatively close, present a conspicuous rewiring of the chromatin landscape (see Fig. 3C). This constitutes a rationale for raising the awareness that generalizing the binding pattern of a GRs to other models, tissues, and even time points, when only one context has been assessed, might be incorrect, and provides the impetus to perform large-scale profiling of GRs while considering time-course trajectories. The third novel observation, and perhaps the most thought-provoking one, concerns the identification of ‘popular’ regulatory regions, frequented by a conspicuous number of GRs. These are lower in number than other regulatory regions bound by one or only a few GRs (Fig. 5A). However, they present a statistically higher evolutionary conservation. Interestingly, the popular regulatory regions we identified were not necessarily super-enhancers, which were previously suggested to be regulatory hubs where signaling cascade converge (Hnisz et al., 2015), and are instead likely promoters according to their proximity to transcription start sites and chromatin-state annotations (Fig. 5D). We propose that different signaling cascades evolved to regulate distinct groups of genes depending on positional information whilst all cooperating to maintain the expression of cell-essential genes by acting on conserved genomic elements that can simultaneously accommodate a multitude of factors, possibly via 3D connections with cell type-specific enhancers. One should keep in mind that a deep form of activity conservation that goes beyond sequence conservation, but concerns a regulatory syntax whereby enhancers can be interpreted by the available GRs, has been shown to exist (Wong et al., 2020). Finally, whether these regions are crowded – that is, GRs bind or occupy the space simultaneously – or bound in different cell types or alternatively along shorter timescales, remains to be determined.

One main limitation currently exists in the construction of the TFBomics. C&R-LoV-U is a bulk-cell population approach. Although single-cell C&T exists (Bartosovic et al., 2021) and few reports of success are available, from extensive contacts with the community and experience in our laboratory we understand that most laboratories still rely on ChIP-seq for mapping dynamic transcription factors or non-DNA-binding co-factors, for which both C&R and C&T are often unsuccessful. While ChIP-seq might still present advantages and result in a higher number of binding events given by the large number of cells used and the application of chemical cross-linking, we favor C&R and C&T as they allow the reduction of material needed of one or two orders of magnitude, and do not require cross-linking, thereby eliminating the artifacts typically introduced by this procedure. Moreover, our recent C&R-LoV-U protocol has been successful with almost all GRs that we have tested so far, even when targeting dynamic, non-DNA-binding proteins (Zambanini et al., 2022). Using C&R-LoV-U is, in our opinion, the current best choice for the large-scale parallel profiling of GRs. scRNAseq and spatial transcriptomics technologies can be synergistically combined with C&R to determine the precise expression of each GRs, thereby increasing the spatial resolution of the C&R findings.

One limitation of our approach and of the current study is that not all datasets, even if they passed quality control, display the same signal-to-noise ratio. This variation is inherent to the technique employed and the quality and specificity of the reagents (e.g. the antibodies) may influence conclusions and data interpretation.

Human embryonic kidney 293T cells (HEK293T) were cultured at 37°C in a humidified incubator with 5% CO₂ in high-glucose Dulbecco's Modified Eagle Medium (41965039, Gibco) supplemented with 10% bovine calf serum (1233C, Sigma-Aldrich) and 1× penicillin-streptomycin (15140148, Gibco).

Animal housing and experimentations were performed under the ethical animal work license obtained by C.C. at Jordbruksverket (Dnr 2456-2019 and 03741-2024), and all procedures and housing were according to the Swedish laws and guidelines. JAX Swiss Outbred mice (The Jackson Laboratory, strain 034608) were used. Housing was in Allentown NexGen individually ventilated cages, floor area 500 cm² with a maximum of four mice per cage. Cages had aspen wood shavings and two types of shredded paper for nesting, and paper tubes were also provided. Temperature was 21±2°C, humidity 45-65% and light cycle 12 h/12 h (07.00 h/19.00 h). Mice had unrestricted access to sterilized drinking water, and ad libitum access to a pelleted and extruded mouse diet. Mice were housed in a barrier-protected specific pathogen-free unit, and this status was monitored and confirmed according to FELASA guidelines by a sentinel program. The mice were free of all viral, bacterial and parasitic pathogens listed in FELASA recommendations (Mähler et al., 2014). Embryo age was determined based on timed mating and vaginal plug observation (E0.5) and confirmed by morphological criteria. Experiments were performed in the afternoons (females euthanized at approximately 13.00 h). Pregnant females were euthanized by cervical dislocation and E10.5 or E11.5 embryos were surgically removed. Tissues were dissected under a dissection stereomicroscope (SZ61, Olympus) and pooled from each litter.

Tissue pools were dissociated to single cells by incubation in TrypLE Express Enzyme (12-604-013, Thermo Fisher Scientific) for 15 min at 37°C on a shaker. For liver tissues, dissociation was achieved with manual pipetting. The cell suspension was resuspended in ice-cold PBS, filtered through a 40 µm cell strainer (KKE3.1, Carl Roth) and further processed for C&R. Three limbs were used per sample, and the equivalent of one liver or one set of branchial arches were used per sample. For HEK293T cells, 500,000 cells were used per sample.

C&R LoV-U was performed as described by Zambanini et al. (2022). Filtered cells were washed three times in nuclear extraction (NE) buffer [20 mM HEPES-KOH (pH 8.2), 10 mM KCl, 0.5 mM spermidine, 0.05% IGEPAL, 20% glycerol, Roche Complete Protease Inhibitor EDTA-Free], then resuspended in 40 µl NE buffer per sample and bound to 20 µl magnetic ConA agarose beads (ABIN6952467, antibodies-online) previously equilibrated in 44 µl binding buffer per sample. After incubation, nuclei and beads were resuspended for 5 min in wash buffer [20 mM HEPES (pH 7.5), 150 mM NaCl, 0.5 mM spermidine in Roche Complete Protease Inhibitor EDTA-Free (COEDTAFRO)] with 2 mM EDTA and split into PCR tubes. Antibody incubation was performed in 200 µl wash buffer with antibody diluted at 1:100 (for antibody and batch information, see Table S1) overnight at 4°C on a rotator. After overnight incubation, samples were washed five times in wash buffer and resuspended in 200 µl of pAG-MN buffer (wash buffer with pAG-MN 120 ng/sample) and incubated for 30 min at 4°C on a rotator. After five washes, digestion was performed for 30 min on ice: started by resuspending the samples in 50 µl of wash buffer with 2 mM CaCl₂. After 30 min, the digestion buffer was removed and moved to tubes containing 3 µl 250 mM EDTA and 250 mM EGTA in H₂O. The digestion reaction on the beads was stopped with 50 µl of 1× Urea STOP buffer (100 mM NaCl, 2 mM EDTA, 2 mM EGTA, 0.5% IGEPAL, 8.5 M urea) and the samples were incubated 1 h at 4°C. Beads were collected on the magnet and liquid transferred to a PCR tube containing digestion buffer. DNA was purified by two subsequent rounds of bead-based cleanup using Mag-Bind TotalPure NGS beads (M1327, Omega Bio-Tek) at 2×, and then resuspended in 20 µl Tris-HCl (pH 7.5).

Library preparation was performed with KAPA Hyper Prep Kit for Illumina platforms (KK8504, KAPA Biosystems) according to the manufacturer's guidelines with modifications: End repair and A-tailing was performed in 0.4× volume reactions with 20 µl of purified DNA. The thermocycler conditions were set to 12°C for 15 min, 37°C for 15 min and 58°C for 25 min to prevent thermal degradation of the shortest fragments. Adapter ligation was performed in 0.4× volume reactions. KAPA Dual Indexed adapters were used at 0.15 µM. A post-ligation clean-up was performed with Mag-Bind TotalPure NGS beads at 1.2× the sample volume. Resuspension was carried out in 10 mM Tris-HCl (pH 8.0). Library amplification was performed in 0.5× volume reactions. The cycling conditions were set as follows: 45 s initial denaturation at 98°C, 15 s denaturation at 98°C, 10 s annealing/elongation at 60°C, 1 min final extension at 72°C, hold at 4°C, with 13 cycles. After amplification, a post-amplification clean-up was performed using NGS beads at 1.2× sample volume. Libraries were then run on an E-Gel EX 2% agarose gel (G402022, Invitrogen) for 10 min using the E-Gel Power Snap Electrophoresis System (Invitrogen). Bands of interest between 150 and 500 bp were excised and purified using the GeneJET Gel Extraction Kit (K0691, Thermo Fisher Scientific) according to the manufacturer's instructions. Libraries were quantified with a Qubit fluorometer (Thermo Fisher Scientific) using a high-sensitivity DNA kit (Q32854, Thermo Fisher Scientific), pooled and sequenced 36 bp pair-end on a NextSeq 550 (Illumina) using the Illumina NextSeq 500/550 High Output Kit v.2.5 (75 cycles) (20024906, Illumina) to an approximate depth of 5-10 million reads per sample.

Quality of reads was assessed using fastqc (version 0.11.9; https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Raw reads were trimmed with bbmap bbduk (version 38.18; Bushnell et al., 2017) to remove adapters, artifacts, [AT]18, [TA]18, and poly G or C repeats. Reads were aligned to the mm10 or hg38 genome with bowtie2 (version 2.4.5; Langmead and Salzberg, 2012) with the options –local –very-sensitive-local –no-unal –no-mixed -no-discordant –phred33 –dovetail -I 0 -X 400. SAMtools (version 1.11; Li et al., 2009) was used to create bam files, fix mate pairing, and for deduplication. Bam files were then filtered to remove reads falling within C&R suspect list regions (Nordin et al., 2023).

Individual replicates were assessed based on duplication rate after trimming (<10%), sequencing depth of mapped, deduplicated, and suspect list filtered bam files (typically 4-10 million reads), fragment profiles (transcription factors show mostly sub- and mono-nucleosome fragments; histone marks typically show mono- and di-nucleosomes), signal at positive control loci (known binding sites for the target), and concordance between biological replicates with genome-wide binning. For example, most replicates cluster together on correlation plots or PCA (Fig. S1A,B). We did not have a hard cutoff regarding replicate similarity before merging, and we did not, for example, call peaks and require a certain amount of overlap, due to our previous investigation into peak overlap making us cautious to use this as a metric of success (Nordin et al., 2024), especially considering dynamic targets such as transcription factors during development.

For peak calling, visualization and signal graphs, replicate bam files were merged with SAMtools into a single file. Bedgraphs were created with BEDtools (version 2.23.0; Quinlan and Hall, 2010) genomecov on pair-end mode and default settings. Peaks were called using MACS2 (version 2.2.6; Zhang et al., 2008) with the options -f BAMPE and -p 1e-3 for narrowPeaks. deepTools (version 3.5.1-0; Ramírez et al., 2014) was used to convert bam files to normalized BigWig files (bamCoverage using -RPGC option and -e to extend reads to fragment length). Peaks were overlapped using BEDtools, with a minimum overlap of 1 bp. Popular regions were defined as regions which overlapped simultaneously 9/12 datasets. Bigwigs were visualized in IGV; tracks shown in the figures are a merge of two biological replicates, normalized to reads per genome coverage, and scaled by factor.

deepTools was used to make signal intensity heatmaps, correlation plots and PCAs. Peak overlaps were performed using BEDtools. Venn diagrams and Upset plots were created using Intervene (version 0.6.4; Khan and Mathelier, 2017). Peak-gene annotation was performed using GREAT (version 4.0.4; McLean et al., 2010) on default settings. Motif analysis was done with HOMER (version 4.11; Heinz et al., 2010) with the option -size given. Gene ontology analysis was performed using Enrichr (Kuleshov et al., 2016). For motifs and gene ontology visualization, color represents significance (P-value) and bubble size fraction of peaks (motifs) or odds ratio (gene ontology). Top five results for each peak set are shown regardless of significance: no q-value cutoff was applied. For conservation analysis, the phastCons 60-way bigwig track for mm10 was downloaded from UCSC genome browser (http://genome.ucsc.edu; Nassar et al., 2023). Peaks were stratified based on how many datasets they were present in. Peaks in each group were treated as bins and the average conservation score within each peak-bin was calculated using deepTools computeMatrix, and then the output matrix was used for plotting and the values for statistical analysis, which was performed using a one-way ANOVA followed by Tukey's post-hoc adjustment. UCSC genome browser was also used for the visualizations shown in Fig. S5D.

Peak regions and, if applicable, bigwig files for ChIP-seq and CUT&RUN datasets were downloaded from the references given. If necessary, peak lists were converted from mm9 to mm10 using the UCSC LiftOver (Hinrichs et al., 2006) with default settings. RNA-seq bulk RPKM data was downloaded from the reference, and these values were plotted for each target. The preprocessed scRNAseq dataset from the mouse embryonic limb atlas (10x_limbdataset) (He et al., 2020), was obtained from the UCSC single-cell browser (https://mouse-limb.cells.ucsc.edu/) as a scanpy-compatible 10x.h5ad file. The Anndata objects were analyzed using Python packages Scanpy (v.1.9.5) (Wolf et al., 2018) and Anndata (v.0.9.2) (https://github.com/scverse/anndata), inspired by Seurat (Butler et al., 2018). The dataset was filtered for stage E11, containing 17,574 cells, and the gene expression of interest was visualized using uniform manifold approximation and projection coordinates matching the original paper. Cell-type cluster assignments were based on cell labels and marker genes from the original study.

The computations and data handling were enabled by resources provided by the National Supercomputer Centre (NSC), funded by Linköping University. Peter Münger at the National Supercomputer Centre is acknowledged for assistance concerning technical and implementational aspects in making the codes run on the Sigma resource. Sequencing was performed by scientists of the Cantù lab operating the Cell Biology Core Facility Unit at the Linköping University.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.204311.reviewer-comments.pdf