Bicoid-nucleosome competition sets a concentration threshold for transcription constrained by genome replication

Abstract

Transcription factors (TFs) regulate gene expression despite constraints from chromatin structure and the cell cycle. Here we examine the concentration-dependent regulation of hunchback by the Bicoid morphogen through a combination of quantitative imaging, mathematical modeling and epigenomics in Drosophila embryos. By live imaging of MS2 reporters, we find that, following mitosis, the timing of transcriptional activation driven by the hunchback P2 (hb P2) enhancer directly reflects Bicoid concentration. We build a stochastic model that can explain in vivo onset time distributions by accounting for both the competition between Bicoid and nucleosomes at hb P2 and a negative influence of DNA replication on transcriptional elongation. Experimental modulation of nucleosome stability alters onset time distributions and the posterior boundary of hunchback expression. We conclude that TF-nucleosome competition is the molecular mechanism whereby the Bicoid morphogen gradient specifies the posterior boundary of hunchback expression.

Figure 1:. Nucleosome positioning at the hb P2 locus is sensitive to Bicoid concentration and cell cycle timing.

A) Bcd binds in a concentration-dependent manner to hb regulatory elements. Shown is a 15 kb region flanking the hb locus including the P2, shadow, and stripe regulatory elements (indicated at top), and counts-per-million (CPM) normalized ChIP-seq for Bcd protein. Mean binding measurements for wild type embryos (top row) is shown in comparison to three mutant lines that express Bcd uniformly across the AP axis at low, medium, or high concentrations. Data from three independent biological replicates were averaged for each genotype. Plotted ChIP-seq coverage ranges from 0-120 CPM for all conditions shown. The genomic locations of the minimal P2 enhancer and extended P2+promoter regions are shown in green (bottom). Scale bar = 1 kb. B) Bcd is required to establish accessible chromatin at the P2 enhancer element. Shown is a 0.6 kb detail of the region upstream of the transcription start site of the hb P2 isoform including the entire minimal P2 enhancer element (green, bottom). ATAC-seq data corresponding to accessible regions and nucleosome dyad positions predicted by NucleoATAC (nuc probability) are shown for five genotypes: wild type, bcd^E1, and low, medium, and high uniform Bcd. Predicted dyad positions referred to in the text are highlighted with arrowheads. Data from three independent biological replicates were averaged for each genotype. Plotted open ATAC coverage ranges from 0-20 CPM for all conditions shown. Scale bar = 50 bp. C) ChIP-Nexus data supports the presence of at least nine Bcd binding sites in the minimal P2 enhancer. Shown is the same genomic region as in (B). Positive and negative contributions to Bcd binding are plotted (blue, Brennan, et al. 2023) and the positions are shown of six binding sites identified through DNase footprinting (A1-3, X1-3), in addition to three extra motifs with PWM scores greater than 80% (E1-3). Two footprints are observed overlapping motifs for Zld (“z”). It is unclear if this represents Bcd binding at these sites, or footprinting of co-immunoprecipitated Zld proteins. The lower portion of the panel shows the positions of two nucleosomes predicted by NuPoP, and a schematic of the P2 enhancer with nine Bcd binding sites. D) Schematic of a minimal hbP2-MS2 reporter: the minimal hb P2 enhancer (green) drives the expression of a chimeric MS2(24)- yellow intron -Gal4 reporter gene through the heterologous Drosophila Synthetic Core Promoter (DSCP, grey). Sequences encoding 24 MS2(v5) hairpins (blue) were introduced to the 5’ end of the yellow intron (yellow) upstream of the gal4 coding sequence (orange). E) Expression of hbP2-MS2 initiates after nuclear import of Bcd and expected replication initiation. Shown is the mean reporter activity of hbP2-MS2 (blue) and mean nuclear intensity of EGFP-Bcd (green) measured over Nuclear Cycle 13 (NC13). All three measurements were performed in separate embryos. The timeline (top) highlights the cell cycle landmarks relative to the period of NC13, timed from the beginning of anaphase 12. F) Nucleosome occupancy over the P2 element fluctuates over the cell cycle. ATAC data over the same region as in (B) is plotted, but for wild type embryos collected at mitotic metaphase (top, NC13+0’) and late interphase (bottom, NC13 + 15’, Blythe and Wieschaus, 2016). Plotted ATAC coverage ranges from 0-20 CPM for all conditions shown. Scale bar = 50 bp.

Figure 2:. Transcriptional onset time of the hbP2-MS2 reporter varies as a function of AP position

A) A heatmap of hbP2-MS2::MCP-GFP fluorescence intensity profiles for ~1400 individual nuclei (n = 19 independent movies) over NC13 is shown. The x-axis corresponds to minutes following anaphase 12, normalized to 20 minutes. The y-axis represents individual nuclei measured and scored as positive for focal expression of MCP-GFP, and sorted by scored onset time. The fluorescence intensity is depicted according to the colorbar on the right. B) Shown are the MS2::MCP-GFP measurements for two representative nuclei, one from the anterior (dark blue) and another from the posterior (light blue) of the hb P2 expression domain. Indicated in red are the portions of the data that yield the onset time and loading rate measurements. While onset times show high variance across the hb P2 expression domain, loading rates and maximal intensities show less variation. C) The fraction of nuclei actively transcribing the reporter does not vary over the region of greatest change in Bcd concentration. The fraction of active nuclei per 2.5% AP bin was calculated and the mean value between movies was plotted (blue ± standard deviation). The gray box in this plot and the following from 60-80% AP position indicates that no movies were taken within this region. D) Reporter onset times show an inverse correlation with Bcd concentration. The onset time of hbP2-MS2 expression was calculated and plotted as a function of AP position (gray points). A rolling average of onset times (n = 100) is shown (blue ± standard deviation) across the 95% C.I. of active nuclei AP positions. E) The average Pol II loading rate of active nuclei does not vary with Bcd concentration. The loading rate of hbP2-MS2 was calculated and plotted as a function of AP position (gray points). A rolling average of loading rates (n = 100) is shown (blue ± standard deviation) across the 95% C.I. of active nuclei AP positions.

Figure 3:. Transcriptional onset time of the hbP2-MS2 reporter is a direct reflection of Bicoid concentration

A) High levels of uniformly expressed Bcd produce uniform Bcd concentrations typically observed at ~38% egg length. Wild type, graded EGFP-Bcd (2x EGFP-Bcd, n = 6 embryos, gray) or uniform EGFP-Bcd (1x tub>uEGFP-Bcd, n = 8 embryos, green) was live-imaged and nuclear fluorescence intensity during NC13 was quantified and plotted as a function of AP position. A rolling average of fluorescence within a 10% AP window is shown ± standard error of the mean. The dotted line indicates the approximate AP position where uniform Bcd expression coincides with graded Bcd (38%). B) Uniform Bicoid drives uniform expression of hbP2-MS2 across the entire AP axis. At left is shown EGFP-Bcd fluorescence (gray) in representative live embryos expressing either 2x EGFP-Bcd (top) or 1x tub>uEGFP-Bcd (bottom). EGFP emissions from a midsagittal confocal section over the last three minutes of NC13 interphase were summed to produce the images shown. Shown at right are hybridization chain reaction images for hbP2-MS2-Gal4 (cyan) in either 2xEGFP-Bcd (top) or 1x tub>uEGFP-Bcd (bottom) fixed embryos. Nuclear DNA (DAPI) is shown in gray. The hbP2-MS2-Gal4 signal extends across the entire AP axis in embryos expressing uniform Bicoid. C) Flattening the Bcd gradient flattens the onset time distributions of hbP2-MS2 to reflect the mean onset time observed at 37% EL. Shown are onset time distributions for hbP2-MS2 measured in uniform Bcd (blue) and graded Bcd (red) as described in panel A. The vertical dashed line indicates the computed AP axis position where the uniform matches the distribution of graded Bcd hbP2-MS2 onset times. Plotted gray points indicate individual measurements from the uniform Bcd experiment only, lines represent rolling averages of onset times (n = 50 ± standard deviation) over the 95% C.I. of AP positions with active nuclei. D) Kr does not regulate the posterior boundary of hbP2-MS2. Gene expression in embryos from a cross between hbP2-MS2-gal4 Kr¹/SM6a individuals was measured by hybridization chain reaction. Left panels show representative hbP2-MS2 reporter activity (cyan) in a wild type (Kr/+) embryo (top) compared with reporter activity in a homozygous Kr mutant embryo (bottom), both staged at early NC14. Right panels show Kr transcript expression (red) in the same embryos pictured at left. The Kr¹ allele is RNA null. DAPI staining is shown (gray) in both panels to facilitate visualization of the specimens.

Figure 4:. A nucleosome-aware stochastic model for onset time prediction produces a posterior boundary for simulated hbP2-MS2 activity.

A) Cartoon schematic of the initial simulation strategy demonstrated for two example nuclei at NC13. Two nuclei (A & B, lower left) in the embryo encounter both spatial and temporal differences in Bcd concentration throughout the cell cycle (top left), which can be expressed as Bcd(x, t). The transition to an “on” state is defined by k_on and in this simulation depends solely on whether Bcd binds successfully to the target. Beginning at anaphase of NC12, the value of Bcd(x, t) is evaluated and k_on is posited for each nucleus at timepoint t_i until a nucleus is scored as “on”. If, at timepoint t, k_on yields a wait-time for binding less than the time-step (Δt), the nucleus is scored as ‘on’ at time t (t_onset). If the calculated wait-time is greater than Δt, the timepoint is incremented to t+Δt and the process repeats. This process yields a t_onset for each nucleus in the field x. B) Modeling Bcd binding with an open chromatin model reflecting Michaelis-Menten kinetics fails to produce spatiotemporal onset time patterns seen in vivo. All plotted data are modeled onset times. Gray points represent individual modeled onset times. The red trace is the rolling mean ± standard deviation of modeled onset times over the 95% C.I. of AP positions with modeled active nuclei. C) Over its physiological range of expression, Bcd is predicted to compete effectively with nucleosomes for occupancy of hb P2 DNA. Shown is the calculated fractional Bcd occupancy predicted by a TF-nucleosome competition model over the region covered by the proximal (brown) and distal (yellow) nucleosomes. The TF-nucleosome competition model was applied using KO_O = 9 nM, K_N = 2610 nM, and L = 700 for each nucleosome, and using n = 4 and n = 5 for the proximal and distal nucleosomes respectively. D) Invoking a nucleosome-aware model for Bcd binding to hb P2 predicts several features of the in vivo onset time distributions. This model captures the distinct posterior boundary to hb P2 expression. All plotted data are modeled onset times. Gray points represent individual modeled onset times. The red trace is the rolling mean ± standard deviation of modeled onset times over the 95% C.I. of AP positions with modeled active nuclei. While the mean and variance of simulated onset times trend towards in vivo measurements, there are notable deviations. This simulation was performed using K_O = 9 nM, K_N = 2610 nM, L = 700, and V_max = 1.

Figure 5:. Zelda-dependent changes in reporter activity are approximated by altering modeled nucleosome stability.

A) Loss of zelda has a moderate impact on chromatin organization at hb P2. Shown are ATAC measurements of accessibility and modeled nucleosome dyad positions for wild type (top) or zelda mutant embryos (bottom). Wild type data is re-plotted from Figure 1B for comparison. Compared with wild type, zelda mutants show a modest reduction in chromatin accessibility over hb P2 and a slight increase in the occupancy of the proximal nucleosome (orange arrowhead). B) Schematic of hb P2 minimal elements highlighting the position of the single endogenous Zelda binding site (orange) within the domain of the proximal nucleosome (below), and relative to Bcd binding sites (blue). To experimentally increase the influence of Zelda on this locus, we designed a mutant reporter that includes one additional Zelda site at the indicated position ([hbP2 + 1x Zld] MS2). C) Alteration of Zelda activity at the hb P2 reporter shifts the posterior boundary of hbP2-MS2 expression. Shown is the fraction of active nuclei in 2.5% egg length bins for hbP2-MS2 in wild type embryos (light blue, n = 19) and zelda-RNAi embryos (dark blue, n = 15), as well as for [hbP2 + 1x Zld] MS2 in wild type embryos (green, n = 17). Dots mark the EC50 of each of the fraction of active nuclei profiles determined by Hill equation fits (Supplemental Information). Decreasing Zelda activity with RNAi shifts the posterior boundary of the reporter anteriorly, while increasing Zelda activity at the reporter by adding a Zelda binding site to it shifts the posterior boundary posteriorly. D) The effect of experimental modulation of Zelda activity on the measured posterior boundaries can be modeled by changing the nucleosome stability parameter L. Shown are the posterior boundary positions (x axis) that are simulated by the model with specific L values (y axis), as determined by the EC50’s of Hill fits to the simulated fraction of active nuclei (Supplemental Information). The solid black dot indicates the maximal AP position achievable by the model when L is set to the lower bound estimate of 50. Colored dots mark the measured posterior boundaries of of hbP2-MS2 in wild type (light blue ± standard deviation), hbP2-MS2 in zelda-RNAi (dark blue ± standard deviation) and [hbP2 + 1x Zld] MS2 in wild type (green ± standard deviation). L = 640, L = 900, and L = 280 allow for prediction of these measurements respectively.

Figure 6:. Accounting for the influence of DNA replication on transcription improves the nucleosome-aware model’s predictions.

A) Cartoon schematic of the revised simulation strategy demonstrated on two example nuclei. As described in the text, the cartoon schematic shows, at left, a single genomic locus in two different nuclei for time points spanning the period of DNA replication, initiated locally at two randomly positioned origins (ori₁ and ori₂, after Blumenthal et al 1974). A consequence of random origin spacing is that a fixed genomic locus completes replication at different times between nuclei. Due to the dependency of transcriptional elongation on completion of DNA replication, random origin spacing is expected to impact observed onset times. At right is a revised model schema, which now requires the enhancer both to bind Bcd and to complete replication prior to switching to the ON state. B) A measurement of mean origin spacing allows for modeling inter-origin distances. As described in the text, defining the parameters of a Gamma distribution using a 9.7 kb average inter-origin spacing and origin frequency of 2 origins per 9.7 kb allows for estimation of measured inter-origin distances, plotted as a histogram. Measured inter-origin distances were manually measured from Blumenthal et al. 1974 Figure 4 and reproduced here. C) Estimation of inter-origin distances allows for prediction of delays in transcriptional onset. The modeled inter-origin distances were converted to replication delays by dividing the distances by 5.3 kb/minute and adding approximately 4 minutes to account for the estimated delay between anaphase and replication onset, plotted as a histogram. D) Incorporating replication delays into the nucleosome-aware model alters the onset time distribution predictions. Shown is one simulation of onset times, plotted as a function of AP position. A rolling average of simulated onset times (n = 50) is shown (blue ± standard deviation) across the 95% C.I. of the AP positions of simulated active nuclei. The simulation was performed with K_O = 9 nM, K_N = 2810 nM, and L = 600. E) A nucleosome- and replication-aware stochastic model accurately predicts the onset time distribution and posterior border of hbP2-MS2 activity. Shown is a set of sampled hbP2-MS2 onset times as a function of AP position, with the rolling average of observed onset times (blue ± standard deviation) and the average of 10,000 onset time simulations (red, dotted) using the Nucleosome + Replication model. For comparison, 10,000 simulations with the Nucleosome-only model is shown (black, dotted). Observed and simulated onset time averages are plotted across the 95% C.I. of the AP positions of active nuclei (observed and simulated, respectively). Simulations were performed with K_O = 9 nM, K_N = 2810 nM, and L = 600 for both the replication-dependent and independent models. F) The nucleosome- and replication-aware stochastic model better accounts for the observed variance in onset time distributions. Plotted are the standard deviations of onset times as a function of AP position for in vivo measurements (blue), and the average standard deviations of 10,000 simulations of the Nucleosome + Replication (red) or Nucleosome-only (black) models. G) Onset times of hbP2-MS2 expression are limited either by replication timing or Bcd concentration. Shown is one instance of the simulation (the same as in D), with points indicating the simulated onset times as a function of AP position. The points have been color coded to indicate whether the onset time was limited by either replication time (blue) or Bcd concentration (red).

Figure 7:. Initiated RNA Pol II waits for the completion of replication and Bcd binding to enter into productive elongation of hb.

A) Inhibition of DNA replication does not impact Bcd binding at the hb locus, but reduces Pol II occupancy over the gene body. Shown is ChIP-seq data over the hb locus for RNA Pol II (dark red), Bcd (blue), and control (gray) comparing occupancy at NC13 in control (untreated) and HU-treated embryos. Data is the average of two independent biological replicates, normalized to read depth (counts per million reads, CPM). The y-axis range for each plot is shown at left (brackets). B) Inhibition of DNA replication does not significantly affect genome-wide Bcd binding. Bcd ChIP-seq data for control (dark blue) and HU-treated (light blue) embryos was standardized and averaged over 1,026 peak regions. The plot shows the average standardized Bcd binding over these peaks. C) Inhibition of DNA replication in NC12-13 primarily affects the entry of RNA Pol II into productive elongation. Upper panels show heatmap representation of standardized RNA Pol II ChIP-seq signal over a set of genes transcribed before large-scale ZGA (pre-MBT genes, Chen et al, 2013). Compared with the control, HU treatment results in loss of RNA Pol II within gene bodies, but retention of a peak near the TSS. The lower panel (C’) is the average of the heatmap representation, with HU-treatment (red) demonstrating reduced Pol II in the gene body but retention of signal over the TSS compared with control (dark blue). D) RNA Pol II at hb is maintained in a paused state in the absence of Bcd. ChIP-seq for initiated RNA Pol II (pSer5) or control was performed on wild type or bcd osk tsl blastoderm stage embryos. In the absence of Bcd, RNA Pol II forms a distinct peak at the hb P2 TSS. Data is an average of three independent biological replicates, normalized to read depth (counts per million reads, CPM). The y-axis range for each plot is shown at right (brackets)..