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 (hbP2) 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 hbP2 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.

Keywords: CP: Developmental biology; CP: Molecular biology; DNA replication; Drosophila; Embryonic patterning; MS2 reporter; chromatin; computational modeling; morphogen; nucleosome; transcription factor.

Figure 1.. Nucleosome positioning at the hbP2 locus is sensitive to Bcd concentration and cell cycle timing

(A) 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 the Bcd protein. Mean binding measurements for wild-type embryos (top row) is shown in comparison to three lines that express Bcd uniformly across the AP axis at low, medium, or high concentrations ($n=3$ averaged per genotype). The locations of the minimal P2 enhancer and extended P2+promoter regions are shown in green (bottom). Bracketed numbers indicate displayed y-axis range in CPM. (B) A 0.6-kb detail of the region upstream of the transcription start site of the hbP2 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 ($n=3$ averaged per genotype). Predicted dyad positions referred to in the text are highlighted with arrowheads. (C) Bcd ChIP-nexus data shown for the same genomic region as in (B). Positive and negative contributions to Bcd binding are plotted (blue), 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 hbP2 enhancer (green) drives the expression of a chimeric MS2(24)-yellow intron-Gal4 reporter gene through the heterologous Drosophila synthetic core promoter (DSCP, gray). 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) Mean reporter activity of hbP2-MS2 (blue) and mean nuclear intensity of EGFP-Bcd (green) measured over 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) ATAC data over the same region as in (B) but for wild-type embryos collected at mitotic metaphase (top, NC13 + 0′) and late interphase (bottom, NC13 + 15′). See also Table S1.

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

(A) A heatmap of hbP2-MS2::MCP-GFP fluorescence intensity profiles for ~1400 individual nuclei scored as positive for hbP2-MS2 ($n=19$ independent movies) over NC13. Nuclei are sorted by onset time. (B) The MS2::MCP-GFP measurements for two representative nuclei, one from the anterior (dark blue) and another from the posterior (light blue) of the hbP2 expression domain. The portions of the data that yield onset time and loading rate measurements are indicated in red. (C) Mean fraction of hbP2-MS2 active nuclei per 2.5% AP bin (blue ± SD). The gray box in this plot and in Figure 2 D–E indicates that no movies were taken within this region. (D) The onset time of hbP2-MS2 as a function of AP position (gray points). A rolling average of onset times ($n=100$) is shown (blue ± SD) across the 95% CI of active nuclei AP positions. (E) The loading rate of hbP2-MS2 as a function of AP position (gray points). A rolling average of the measurements is shown as in Figure 2D. See also Figures S1 and S2.

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

(A) NC13 nuclear fluorescence of graded EGFP-Bcd (2× EGFP-Bcd, $n=6$ embryos, gray) and uniform EGFP-Bcd (1× tub>uEGFP-Bcd, $n=8$ embryos, green). A rolling average of fluorescence within a 10% AP window is shown ± SEM. The dotted line indicates the approximate AP position where uniform Bcd expression coincides with graded Bcd (38%). (B) Left: midsagittal NC13 EGFP-Bcd fluorescence (gray) in representative live embryos expressing either 2× EGFP-Bcd (top) or 1× tub>uEGFP-Bcd (bottom). Right: hybridization chain reaction images for hbP2-MS2-Gal4 (cyan) in either 2× EGFP-Bcd (top) or 1× tub>uEGFP-Bcd (bottom) fixed embryos. Nuclear DNA (DAPI) is shown in gray. Embryo lengths are approximately 500 μm. (C) Onset time distributions for hbP2-MS2 measured in uniform Bcd (blue, $n=949$ nuclei) and graded Bcd (red, $n=413$ nuclei) as described in (A). The vertical dashed line indicates the computed AP axis position where the uniform matches the distribution of graded Bcd hbP2-MS2 onset times (37% EL). Plotted gray points indicate individual measurements from the uniform Bcd experiment only; lines represent rolling averages of onset times ($n=50\pm \mathrm{S}\mathrm{D}$) over the 95% CI of AP positions with active nuclei. (D) Left: representative hbP2-MS2 reporter activity measured with hybridization chain reaction (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: Kr transcript expression (red) in the same embryos pictured at left. DAPI staining is shown (gray). See also Figure S3. Embryo lengths are approximately 500 μm.

Figure 4.. An allosteric 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 and B, lower left) encounter spatiotemporal differences in Bcd concentration throughout the cell cycle (top left, $\mathrm{B}\mathrm{c}\mathrm{d}(x,t)$). The transition to an “on” state is defined by $k_{on}$ and depends solely on whether Bcd binds to the target. Beginning at anaphase 12, $\mathrm{B}\mathrm{c}\mathrm{d}(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 ($\Delta t$), the nucleus is scored as “on” at time $t(t_{onset}$). If the calculated wait-time is greater than $\Delta t$, the timepoint is incremented to $t+\Delta t$, and the process repeats. This process yields a $t_{onset}$ for each nucleus in the field $x$. (B) Onset times modeled with an open chromatin mechanism reflecting Michaelis-Menten kinetics. Gray points represent individual modeled onset times. The red trace is the rolling mean ± SD of modeled onset times over the 95% CI of AP positions with modeled active nuclei. (C) 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 $K_O=9\mathrm{n}\mathrm{M},K_N=2610\mathrm{n}\mathrm{M}$, and $L=700$ for each nucleosome and using $n=4$ and $n=5$ for the proximal and distal nucleosomes respectively. (D) Onset times modeled with an allosteric mechanism for Bcd binding to hbP2. Gray points represent individual modeled onset times. The red trace is the rolling mean ± SD of modeled onset times over the 95% CI of AP positions with modeled active nuclei. This simulation was performed using $K_O=9\mathrm{n}\mathrm{M},K_N=2610\mathrm{n}\mathrm{M},L=700$, and $V_{max}=1$. See also Figure S4 and Data S1.

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

(A) ATAC measurements of accessibility and modeled nucleosome dyad positions for wild-type (top) or zelda mutant embryos (bottom). Wild-type data are replotted from Figure 1B for comparison. The orange arrowhead marks a slight increase in the occupancy of the proximal nucleosome. (B) Schematic of hbP2 minimal elements highlighting the positions of Zelda binding sites (orange) within the domain of the proximal nucleosome (below) and relative to Bcd binding sites (blue). Top: enhancer of the hbP2-MS2 reporter. Bottom: enhancer of the [hbP2 + 1x Zld] MS2 reporter. (C) The fraction of active nuclei in 2.5% egg length bins for hbP2-MS2 in wild-type embryos (light blue ± SD, $n=19$) and zelda-RNAi embryos (dark blue ± SD, $n=15$), as well as for [hbP2 + 1x Zld] MS2 in wild-type embryos (green ± SD, $n=17$). Dots mark the EC50 of each of the fraction of active nuclei profiles determined by Hill equation fits (Data S1). (D) The posterior boundary positions (x axis) that are simulated by the model with specific $L$ values (y axis), as determined by the EC50s of Hill fits to the simulated fraction of active nuclei (Data S1). 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 hbP2-MS2 in wild-type (light blue ± SD), hbP2-MS2 in zelda-RNAi (dark blue ± SD), and [hbP2 + 1x Zld] MS2 in wild-type (green ± SD). $L=640,L=900$, and $L=280$ allow for prediction of these measurements, respectively. See also Figure S5 and supplemental information: Data S1.

Figure 6.. Accounting for the influence of DNA replication on transcription improves the allosteric model’s predictions

(A) Cartoon schematic of the revised simulation strategy demonstrated on two example nuclei. The cartoon on the left shows a single genomic locus in two different nuclei for time points spanning DNA replication, initiated locally at two randomly positioned origins (ori₁ and ori₂, after Blumenthal et al.). Due to the dependency of transcriptional elongation on completion of DNA replication, random origin spacing is expected to impact observed onset times. On the 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 Gamma distribution with a 9.7-kb average inter-origin spacing and origin frequency of 2 origins per 9.7 kb, plotted as a histogram (gray). Measured inter-origin distances were manually measured from Figure 4 of Blumenthal et al. and reproduced here (purple). (C) Modeled replication delays calculated from the inter-origin distance model. (D) Onset times simulated by the allosteric + replication model. A rolling average of simulated onset times ($n=50$) is shown (blue ± SD) across the 95% CI of the AP positions of simulated active nuclei. The simulation was performed with $K_O=9\mathrm{n}\mathrm{M},K_N=2810\mathrm{n}\mathrm{M}$, and $L=600$. (E) A set of sampled hbP2-MS2 onset times as a function of AP position, with the rolling average of observed onset times (blue ± SD), the average of 10000 onset time simulations (red, dotted) using the allosteric + replication model, and the average of 10000 onset time simulations with the allosteric-only model (black, dotted). Observed and simulated onset time averages are plotted across the 95% CI of the AP positions of active nuclei (observed and simulated, respectively). Simulations were performed with $K_O=9\mathrm{n}\mathrm{M},K_N=2810\mathrm{n}\mathrm{M}$, and $L=600$ for both the replication-dependent and independent models. (F) The SD of onset times as a function of AP position for in vivo measurements (blue), and the average SD of 10000 simulations of the allosteric + replication (red) or allosteric-only (black) models. (G) One instance of the simulation (the same as in (D)), with points 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) 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 are the average of two independent biological replicates, normalized to read depth (CPM). (B) The average standardized Bcd ChIP-seq data for control (dark blue) and HU-treated (light blue) embryos over 1026 peak regions. (C) Upper panels: heatmap representation of standardized RNA Pol II ChIP-seq signal over a set of genes transcribed before large-scale ZGA (pre-MBT genes). Lower panel (C′): the HU-treatment (red) and control (dark blue) averages of the heatmap representation. (D) ChIP-seq data (average of $n=3$ biological replicates, normalized to CPM) for initiated RNA Pol II (pSer5) and control in wild-type and bcd osk tsl blastoderm-stage embryos. See also Figures S6 and S7.