Bulk-level maps of pioneer factor binding dynamics during the Drosophila maternal-to-zygotic transition Open Access

In Drosophila, zygotic genome activation (ZGA) begins with the transcription of a handful of genes during its minor wave, followed by a major wave when thousands of genes are transcribed (Harrison et al., 2023; Pritchard and Schubiger, 1996; Schulz and Harrison, 2019; Vastenhouw et al., 2019). The transcription factor (TF) Zelda (Zld) has the ability to bind nucleosomal DNA (Garcia et al., 2019; McDaniel et al., 2019) and subsequently to facilitate the binding of other TFs (Satija and Bradley, 2012; Li et al., 2014; Foo et al., 2014; Sun et al., 2015; Schulz et al., 2015; Yáñez-Cuna et al., 2012; Xu et al., 2014), which are the two defining features of a special class of TFs known as pioneer factors (Iwafuchi-Doi and Zaret, 2014; Zaret and Carroll, 2011; Zaret and Mango, 2016). The maternally encoded TF GAGA factor (GAF) also possesses pioneer-like properties (Blythe and Wieschaus, 2015, 2016; Gaskill et al., 2021; Judd et al., 2021; Schulz et al., 2015; Sun et al., 2015).

The relationship between the concentration of TFs such as these pioneer factors and gene expression remains an open question (Ay and Arnosti, 2011; Kim and O'Shea, 2008). Fluorescent imaging of live and fixed tissues can give relative concentrations of TFs, which are often used to predict TF activity by assuming some nonlinear relationship between concentration and activity (Jaeger et al., 2004; Kanodia et al., 2012; Manu et al., 2009; O'Connell and Reeves, 2015; Papatsenko and Levine, 2005, 2011; Zinzen et al., 2006). However, validation of these relationships using dynamic quantitative data in live cells is currently lacking. Furthermore, the binding of the TFs to DNA – which is required for transcriptional regulation – is dependent not only on their nuclear concentration, but also on factors such as chromatin accessibility and saturation kinetics. On the other hand, recent studies have identified hubs of TF binding at DNA sites, which may provide more direct measures of TF activity (Dufourt et al., 2018; Mir et al., 2017, 2018; Yamada et al., 2019). Thus, knowing only total concentration is inadequate to predict TF binding, and direct measurements of TF hubs lack generalizability; an input/output map between the two, based on quantitative measurements of concentration and binding, is required.

To bridge this gap, we used raster image correlation spectroscopy (RICS) to quantify the concentration and binding of the pioneer-like factors Zld and GAF in live embryos over the course of nuclear cycles (ncs) 10-14 (Al Asafen et al., 2024; Brown et al., 2008; Digman and Gratton, 2009; Digman et al., 2005a,b; Schloop et al., 2025). These measurements allowed us to construct dose/response relationships suggesting that the binding of Zld to transcriptionally active sites is dependent on factors other than its concentrations alone, whereas the GAF concentration is the primary driver of its binding. Furthermore, we found that GAF must bind and saturate its sites in the inactive regions of the DNA before it can bind to active DNA sites, resulting in a delay in its role in ZGA. Similar approaches can be used to obtain a comprehensive quantitative picture of dynamics of other TFs for gene regulation studies.

Raster-scanned confocal images have a fast-scanning direction (pixel-to-pixel) and a slow-scanning direction (line retracing; Fig. 1A). RICS analysis uses intensity fluctuations of GFP-tagged molecules, correlated in time and space, to build an autocorrelation function (ACF) in two dimensions, corresponding to the fast (Δx) and slow (Δy) scanning directions, respectively (Fig. 1B; see Materials and Methods). The amplitude of the ACF, A, is inversely proportional to the GFP concentration (Fig. 1C), while the shape of the ACF, especially in the slow direction, is determined by the fraction of GFP that is freely diffusible versus that which is immobile (Fig. 1D). Furthermore, cross-correlations in the fluctuations between GFP and RFP (tagged to His2Av) allow us to calculate the fraction of GFP-tagged molecules that are bound to the same structure as the His2Av-RFP (likely to be DNA; Bacia and Schwille, 2007; Weidemann et al., 2002; Fig. 1E). Therefore, RICS allows us to quantify not only the absolute concentration of GFP-tagged molecules, but also the fraction that is freely diffusible, the immobile fraction, and the fraction correlated to His2Av (Fig. 1F).

We performed RICS analysis on the nuclei of blastoderm-stage embryos expressing sfGFP-Zld (Hamm et al., 2017) (Fig. 2A and Movie 1) to measure the dynamics of the nuclear concentration of Zld over time. We found the total nuclear concentration of Zld increased from one nuclear cycle to the next during nc 10 to 14 (Fig. 2B). Increase in Zld levels has been reported previously using immunoblotting (Harrison et al., 2010; Nien et al., 2011). Furthermore, the nuclear concentrations of Zld vary significantly within each nc, as more sfGFP-Zld enters nuclei after mitosis and fills it during the interphase. The longer duration of nc 14 allows Zld levels to reach a steady state, unlike the other, shorter, ncs. In contrast, we observed that the nuclear Zld concentration drastically decreases from the end of one interphase to the beginning of the next (during mitosis), in agreement with previous observations that, unlike most pioneer factors, Zld is not mitotically retained on the chromosomes (Dufourt et al., 2018).

RICS analysis also allowed us to determine φ, the fraction of Zld that is immobile (or nearly so) due to either binding to immobile structures (such as DNA) or to forming large aggregates with very low diffusivity (Eqn 4). The immobile fraction of Zld remains nearly the same from nc 10 to nc 13, after which it reaches a steady state at a lower value of ∼0.5 at nc 14 (Fig. 2C). Fifty percent of the Zld population was found to be immobile in single-molecule imaging (Mir et al., 2018).

In each nc, at the beginning of interphase, the immobile concentration of Zld, calculated as the product of total nuclear concentration and the immobile fraction, increases, perhaps because the chromatin decondenses and Zld binds the replicated DNA (Antonin and Neumann, 2016; Blythe and Wieschaus, 2016; McDaniel et al., 2019; Yuan et al., 2014), and it starts decreasing after reaching a maximum as chromatin starts condensing before mitosis (Fig. S1A). This is consistent with the rapid formation of dynamic hubs of Zld after mitosis (Mir et al., 2018), the hubs likely being a population included in our immobile concentration measurements. It should be noted that the presence of hubs results in Zld particles with higher brightness than the Zld monomer. As such, the presence of hubs, in principle, could alter our measurements of total Zld nuclear concentration and of the immobile fraction. We have quantified the average particle brightness and found it to be roughly constant over time, suggesting that hubs of Zld have little effect on our measurements (see Materials and Methods and Fig. S2).

We also used cross-correlation RICS (ccRICS) to determine ψ, the fraction of the pioneer factors that correlate with His2Av-RFP (Eqn 9 in the Materials and Methods). In Drosophila, the enrichment of the histone variant His2Av has been found to correlate with transcriptional potential (Leach et al., 2000; Mavrich et al., 2008), and therefore, the pool correlated with His2Av-RFP (hereafter referred to as the active pool) is likely responsible for pioneering gene activation. As the immobile fraction always measured as larger than the active fraction (Fig. 2C,D), we inferred that the immobile fraction must be composed of two pools: one that correlates with His2Av-RFP and one that does not. Hereafter, the pool that does not correlate with His2Av but is immobile, either due to binding to inactive chromatin regions or due to low diffusivity from the formation of large clusters, is referred to as the inactive pool. We found the active fraction decreases from nuclear cycle to nuclear cycle (Fig. 2D). Consequently, the active concentration of Zld, which is the product of the total nuclear concentration and the active fraction of Zld, also decreases, albeit weakly, due to the increasing total nuclear concentration (Fig. S1B). Our results agree with the slight reduction observed in the number of Zld peaks from nc 13 to nc 14 in ChIP-Seq (Harrison et al., 2011). The Zld ChIP-Seq peaks likely represent the binding in the active regions as tightly packed inactive regions are less likely to be identified (Nakato and Shirahige, 2016). Roughly 5 min into nc 14, the concentration of active Zld reaches a steady value around 20-25 nM (Fig. S1B).

Our measurements lend themselves to a general analysis in which the relationship between Zld total nuclear concentration and its binding could be constructed. According to standard thermodynamic equilibrium binding models (see Materials and Methods), an increase in the total Zld concentration should lead to an increase in the bound Zld concentration. However, our data suggest that the active concentration of Zld is a weakly decreasing function of total Zld (Fig. 2E) in which nc 10-12 has a high active concentration of Zld despite the low total Zld concentration, while in nc 14, there is a high total Zld concentration and a slightly lower active concentration. As such, our results suggest that Zld binding is not influenced by Zld concentration alone; additional factors, such as chromatin structure, must be taken into account. Using our data and ChIP-seq results reported by Harrison et al. (2011), the dissociation constant for the active population (K_D,a) for Zld can be estimated to be ∼17 nM (see Materials and Methods). Fitting a Hill function to our data suggests that the Zld sites accessible for binding reduce over time (Fig. 2F), which might result from the chromatin becoming more densely packed during the later cycles (Harrison et al., 2011; Li et al., 2014; Lowenhaupt et al., 1983). On the other hand, the inactive concentration of Zld increases with time (Fig. S1C), suggesting the inactive fraction may result from binding that is saturating in nature. Fitting a Hill function with a Hill coefficient of 2, due to the sigmoidal nature of the data points, yielded the dissociation constant for the inactive population (K_D,i) of ∼24 nM and a max concentration of 34 nM (Fig. 2E). These values of K_D suggest that the binding of Zld to DNA or to immobile structures including aggregates is moderately strong. Overall, the dose/response relationship between total and bound Zld concentration is consistent with the previously observed hubs of Zld bound to DNA. Zld molecules within these hubs have short residence times on DNA (τ∼5 s; Mir et al., 2018). This is lower than our estimation (see Materials and Methods), which might result from the high local concentration of Zld within the hub (Mir et al., 2018) or from the local Zld diffusivity (Mirny et al., 2009).

Our data suggest that Zld may not retain its ability to pioneer chromatin accessibility throughout ZGA. Consistent with this, it was recently shown that, in larval type II neuroblasts, Zld binding is influenced by chromatin accessibility (Larson et al., 2021). Thus, our observation may imply that Zld primarily regulates transcription during the minor wave. Indeed, the removal of Zld activity affects early gene expression patterns, resulting in a delay in, but not complete loss of, the transcription of patterning genes (Liang et al., 2008; Nien et al., 2011). Therefore, our results suggest that, during the minor wave of ZGA (nc 10-13), a high concentration of bound Zld might be needed to facilitate the binding of patterning factors present at low levels. On the other hand, during the major wave at nc 14, either the high level of patterning factors, or the presence of other pioneer-like factors, may be sufficient to drive the patterning factor binding and target gene expression when the bound Zld levels have decreased. This raises the question of whether other factors are present during the major wave that may continue to facilitate the binding of other developmental TFs.

GAF, which is encoded by the Trithorax-like (Trl) gene, plays an essential role in ZGA, along with Zld (Farkas et al., 1994; Gaskill et al., 2021). GAF motifs along with Zld motifs are enriched in the highly occupied target (HOT) regions characterized by open chromatin and bound by many TFs (The modENCODE Consortium et al., 2011; Kvon et al., 2012; Slattery et al., 2014). GAF possesses the properties of a pioneer factor, such as the ability to bind nucleosomal DNA and to create regions of chromatin accessibility by functioning with chromatin remodelers to facilitate the binding of other TFs (Fuda et al., 2015; Gaskill et al., 2021; Judd et al., 2021; Moshe and Kaplan, 2017; Tang et al., 2022; Tsukiyama et al., 1994).

To measure the dynamics of the nuclear concentration of GAF over time, we performed RICS analysis on the nuclei of blastoderm-stage embryos expressing GAF-sfGFP (Gaskill et al., 2021) (Fig. 3A and Movie 2). The results suggest that the total nuclear concentration of GAF remained nearly constant and very low (∼5-10 nM) from nc 10-13, then showed a sudden increase during nc 14 to ∼30 nM (Fig. 3B). We saw that the nuclear concentration at the beginning of an interphase was similar to that at the end of the previous interphase (Fig. 3B). This observation is consistent with previous work showing GAF is mitotically retained on the chromosomes (Gaskill et al., 2021), and is in contrast to Zld (Fig. 2B), which is not mitotically retained (Dufourt et al., 2018).

The immobile fraction of GAF remained nearly the same during nc 10-13, then decreased slightly during nc 14 (Fig. 3C). The sudden increase in GAF total nuclear concentration, together with only a slight decrease in its immobile fraction, resulted in an increase in the immobile concentration of GAF during nc 14 (Fig. S3A). Although the active fraction of GAF, the pool expected to be bound near the actively transcribed genes, remained approximately the same, on average, from nc 10-14, it varied significantly within each nc (Fig. 3D).

As with Zld, we sought to use our data to map the relationship between GAF nuclear concentration and binding. We noted that, unlike Zld, both active and inactive concentrations of GAF increase with an increase in the total concentration of GAF (Fig. 3E). Our results agree with the increase in GAF ChIP-seq peaks from nc 9 to 14 observed previously (Gaskill et al., 2021). Both pools appeared to be saturating functions of total GAF concentration, with the inactive pool saturating quickly at a low overall concentration (Fig. 3E and Fig. S3B,C). In contrast, the active pool required higher levels of total GAF to saturate and had a higher capacity. Fitting Hill functions to the data bore out these observations: the inactive pool had a K_D,i of 0.4 nM and a max concentration of 9 nM, while the active pool had a K_D,a of 5 nM and a max concentration of 35 nM. The K_D values represent a strong affinity of GAF for the binding sites, roughly 4-fold and 60-fold greater than for Zld, consistent with previous reports of stable GAF/DNA binding (Bellec et al., 2022; Espinás et al., 1999; Tang et al., 2022). It appears that, as GAF concentration slowly increases, the majority is apportioned to the inactive pool due to its high affinity. However, because of the low capacity of the inactive pool, upon entering nc 14, the active pool is suddenly able to increase. If the active pool corresponds to transcriptionally active regions of the DNA (i.e. euchromatin), the relative affinities and capacities ensure that GAF acts as a pioneer factor solely during the later, major, wave of ZGA. Furthermore, it has been shown that GAF associates with the GA/CT-rich repeats in the heterochromatin regions throughout cell cycles (Gaskill et al., 2023; Raff et al., 1994), potentially driving transcriptional silencing and the euchromatin association activation during ZGA (Gaskill et al., 2023). As such, the inactive pool may represent GAF bound to heterochromatin. As with Zld, GAF hubs with higher brightness could affect our measurements. Quantification of the brightness suggested small corrections to our data, which did not affect our overall conclusions (see Materials and Methods and Fig. S4).

Our measurements of the biophysical parameters of the pioneer factors and their dynamics allowed the construction of input/output maps between total concentration of the factors and their DNA binding. While a subset of the pioneer factor bound regions has been found to stay inaccessible (Freund et al., 2024; Gibson et al., 2024), binding of pioneer factors to DNA is an indispensable step for establishing and maintaining chromatin accessibility. Following this, binding of a correct combination of TFs can lead to the activation of the accessible regions (Brennan et al., 2023). Using the input/output maps, we can infer that Zld acts as a global activator of early genes enabled by the high active concentration of Zld during the minor wave of ZGA; GAF comes into play during the major wave, as indicated by the high active concentration of GAF during this time (Fig. 4). This is supported by the fact that the regions that gain accessibility early during ZGA are enriched for Zld binding, whereas those that gain accessibility late are enriched for GAF binding (Blythe and Wieschaus, 2016). A similar quantitative approach, in which accurate measurements of biophysical parameters of TFs and their dynamics enable the construction of input/output map between TF concentration and transcriptional activity, is expected to aid other gene expression studies.

The fly strain used for the Zld imaging is sfGFP-zld; His2Av-RFP (II). The sfGFP-zld mutant allele was generated using Cas9-mediated genome engineering by Hamm et al. (2017). The fly strain used for the GAF imaging is His2Av-RFP (II); GAF-sfGFP(C) (III). The GAF-sfGFP(C) (III) mutant allele was generated using Cas9-mediated genome engineering by Gaskill et al. (2021).

Flies were raised on standard cornmeal-molasses-yeast medium at 25°C. Fly cages were prepared with desired fly strains and kept at room temperature for 2 days before imaging. Grape juice agar plates streaked with yeast paste were placed onto the bottoms of the cages for oviposition. For imaging, the flies were allowed to lay eggs for 1 h after which the plates were removed for embryo collection. The embryos were transferred from the plates to mesh baskets, dechorionated using bleach for 30 s and washed with deionized water to remove residual bleach (Carrell and Reeves, 2015).

After dechorionation, the embryos were mounted in 1% low melting point agarose (IBI Scientific, IB70051) in deionized water on a glass bottom Petri dish (MatTek, P35G-1.5-20-C) and deionized water was poured into the Petri dish over the solidified low melting point agarose to cover the samples. The low melting point agarose rendered mechanical stability to the embryos while remaining transparent submerged in water (Huisken and Stainier, 2009; Kaufmann et al., 2012). All images were collected on a Zeiss LSM 900 confocal laser scanning microscope. For the image acquisitions, C-Apochromat 40×/1.2 water immersion Korr objective, 488 nm laser for sfGFP, 561 nm laser for RFP and GaAsP-PMT detector were used. The detector was operated at 650 V with 1× gain and 0% offset. Emission was detected in the range of 410-546 nm for sfGFP and 595-700 nm for RFP. All the images had 1024×1024 pixel resolution. The images were collected at 5× zoom, resulting in pixel size of 31.95 nm and a pixel dwell time of 2.06 µs. This corresponds to a frame time of 5.06 s and a line time of ∼5 ms (ratio of frame time and number of rows in the image). Image acquisition was started when the embryos were at nc 10 and continued until gastrulation, as indicated by the nuclear morphology.

Raster image correlation spectroscopy (RICS) analysis, a derivative of fluorescence correlation spectroscopy, entails constructing autocorrelation functions (ACFs) from imaging data and fitting ACF models to these data-derived ACFs. We performed these analyses according to previous protocols (Al Asafen et al., 2024; Schloop et al., 2025). In brief, live imaging time courses were background subtracted and divided into groups of 7-12 frames. The frames within a group were averaged together to create an averaged frame for that group, which was used for two purposes. First, the nuclei in the averaged frame were segmented according to a watershed algorithm. This segmented nuclear mask was used for each frame in the group. Second, the immobile variation within each frame was then removed by subtracting by the averaged frame on a pixel-by-pixel basis, then the scalar average intensity value of the averaged frame was added back (Digman et al., 2005a,b). Two dimensional (2D) RICS autocorrelation functions of the nuclear fraction of each frame were built using a fast Fourier transform protocol in Matlab, and these ACFs were averaged together for all frames in a given group. The result was a time series of 2D ACFs, each of which corresponded to a given grouping of 7-12 frames. Background subtraction was performed on the fly by examining the histogram of intensities and fitting a Gaussian to the lowest intensity pixels in the image.

We mounted diffraction-limited TetraSpeck beads (0.1 µm, T7279, Invitrogen) in 1% low melting point agarose on a glass bottom Petri dish following the same protocol as the embryos. Similarly, the agarose was immersed in deionized water upon solidification. The microscope parameters, such as the pixel size and frame time, were identical as the RICS acquisitions. We acquired z-stacks with 0.05 µm distance between the slices. We acquired six z-stacks: three different locations in the agarose and two different batches of beads mounted in the agarose.

The beads were then detected in the average image of the z-stack. A 3D Gaussian was fitted to the intensity of the beads using fmincon. The centroid coordinates of each detected bead and the highest intensity z-plane for the corresponding bead were used as the initial guesses for the centers in x, y and z directions. This was carried out separately for the two channels (Fig. S5A,B). The difference between the centers (x, y and z) for the two channels, as determined from the fit, were calculated for each bead (Fig. S5C-D). The full width at half maximum (FWHM) in the z-direction was calculated separately for the two channels (Fig. S5E).

The equation was solved and adjustable parameters (dissociation constant and binding site concentrations) were determined using the least-squares fitting algorithm lsqcurvefit and the global optimization solver MultiStart in Matlab.

Pooled data from nc 8, 13 and 14 identified 12,135 Zld ChIP-Seq peaks (Harrison et al., 2011). We assumed that, during the early ncs, the chromatin is fairly open and all the ChIP-Seq identified sites are accessible for binding, such that the number of Zld sites accessible for binding during nc 11 is also assumed to be 12,135. Given that the nuclear volume during nc 11 is roughly 435 μm³ (Foe and Alberts, 1985), the concentration of Zld sites accessible for binding, averaged over the whole nucleus, is roughly, ⁠. However, given that Drosophila is diploid, we set L=2×46.34=92.68 nM. Using these estimates, Eqn 10 and our data for nc 11, we arrive at K_D,a=17 nM. Using this K_D and our data in Eqn 10, the change in Zld sites accessible for binding over time was estimated.

Dissociation constants are typically defined as ⁠; where k_on and k_off are the on and off rates for the TF-immobile structure interaction. For protein-ligand binding events in which the on-rate is limited by three dimensional diffusion, 10⁻⁴ nM⁻¹ s⁻¹<k_on<10⁻³ nM⁻¹ s⁻¹ (Wittrup et al., 2019). Using this rule of thumb for the k_on value and K_D∼20 nM (valid for both the active and inactive pools of Zld), the residence time, ⁠, can be estimated to be τ∼8 min−50 s. However, if the diffusivity is higher, as would be the case of one-dimensional sliding along the DNA (Mirny et al., 2009), k_on is closer to 10⁻¹ nM⁻¹ s⁻¹, resulting in τ∼0.5 s. Owing to this variability using the rule of thumb, we used the Smoluchowski equation (Mirny et al., 2009) for the Zld-DNA interaction. The diffusivity of Zld obtained from our analysis (D) ∼2 µm²/s. According to the Smoluchowski equation (Mirny et al., 2009): k_on=4πDba, where cross-section of the binding reaction (b)=0.34 nm and the fraction of the molecular surface of the protein that contains the reactive binding interface (a) ∼0.2-0.5. Therefore, k_on∼0.001-0.003 nM⁻¹ s⁻¹, resulting in τ∼19−49 s, which agrees with the lower limit of the value calculated using the rule of thumb.

In practice, the shot noise is equal to the RICS ACF amplitude times the intensity squared, AI², subtracted from the variance described above, ⁠. The shot noise and intensity at each time point in the time course was used to compute the best-fit line (see Fig. S2B for Zld and S4B for GAF). The equation used was ⁠, where b is a free parameter added for robustness of fit and was found generally to be close to zero. For Zld, S=363; for GAF, S=413.

From these, corrected values of all concentrations, the free fraction and the uncorrelated fraction were calculated. These corrections were performed at two levels. In the first level, we assumed all variations within nc 14 were within the range of being roughly constant. In other words, only variations outside of this range were corrected for, meaning that Q_min was set as the maximum value observed during nc 14 (see Fig. S4C). This resulted in corrections made for the increase in brightness in the final point of nc 12 and the final four points of nc 13. While these corrections had a quantitative effect on our results, the corrected relationships between bound and free GAF are maintained (Fig. S4D,E).

We thank Melissa M. Harrison for kindly providing fly stocks and for helpful discussion of the manuscript.

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