Comprehensive Single Molecule View of Transcriptional Dynamics in Development

Comprehensive Single Molecule View of Transcriptional Dynamics in Development | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Comprehensive Single Molecule View of Transcriptional Dynamics in Development View ORCID Profile Thomas W Tullius , Tohn Borjigin , Michael Levine doi: https://doi.org/10.1101/2025.10.28.685222 Thomas W Tullius 1 Lewis-Sigler Institute Department of Molecular Biology Princeton University Princeton , NJ 08544 Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Thomas W Tullius Tohn Borjigin 1 Lewis-Sigler Institute Department of Molecular Biology Princeton University Princeton , NJ 08544 Find this author on Google Scholar Find this author on PubMed Search for this author on this site Michael Levine 1 Lewis-Sigler Institute Department of Molecular Biology Princeton University Princeton , NJ 08544 Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: msl2{at}princeton.edu Abstract Full Text Info/History Metrics Preview PDF Abstract Dynamic patterns of gene expression drive animal development. These transcriptional dynamics depend on the interplay of diverse developmental enhancers within complex regulatory landscapes. Here, we employ a single-molecule genomic method, Fiber-Seq, to capture the simultaneous transcriptional and regulatory states of millions of 20-30 kb chromatin fibers in the 2-4hr Drosophila embryo. These snapshots of genes and their surrounding regulatory context reveal that the basic unit of transcription at highly transcribed genes is a convoy of 2-14 tightly spaced polymerases (Pol II). Further, a subset of templates exhibits hyperbursting, whereby most or all nucleosomes are evicted from the gene body facilitating maximum rates of transcription. We demonstrate that hyperbursting is achieved through cooperation of numerous cis-regulatory elements, and that the resulting nucleosome eviction appears to trigger a novel silencer in the ftz 3’ UTR. We anticipate that similar mechanisms are used by vertebrate processes requiring intense transcription like somitogenesis and erythropoiesis. Introduction Most previous studies of gene regulation in development focus on interactions of individual enhancers with their target promoters. There are numerous examples, including long-range regulation of the Sonic Hedgehog ( Shh ) gene by the distal ZRS enhancer in developing limb buds of mammalian embryos ( 1 ). However, developmental control genes such as Shh contain complex regulatory landscapes containing numerous enhancers ( 2 ). There are several documented examples of cis-regulatory interactions influencing the timing of gene expression, such as limb patterning by the Hoxd locus in mice ( 3 ). There are also examples in the Drosophila embryo, including enhancer switching in the Brinker locus ( 4 ) and the timing of segmentation stripes of even-skipped expression ( 5 ). These examples of co-dependencies of cis-regulatory elements in development are anecdotal as they depend on painstaking functional analyses in transgenic embryos. Here, we employ a high-throughput single-molecule method, Fiber-Seq ( 6 , 7 ), to examine the relationship of enhancer activities and transcriptional states for hundreds of developmental control genes on millions of individual chromatin fibers in the early Drosophila embryo. Computational analyses of these comprehensive datasets identified “hyperbursts”, an extreme state of transcriptional activation whereby convoys of 2-14 tightly spaced Pol II complexes evict most or all nucleosomes across gene bodies. Hyperbursts depend on the cooperation of multiple cis-regulatory elements, including those dedicated to transcriptional bursting. A newly identified silencer element in the ftz 3’ UTR appears to be triggered by hyperbursts and nucleosome depletion. This element is the source of an anti-sense RNA that is expressed in regions of the embryo where ftz is known to be repressed ( 5 ). We discuss the implications of these findings with respect to mammalian developmental and disease processes requiring intense periods of transcription such as somitogenesis ( 8 ) and erythropoiesis ( 9 ). Results We employed a high throughput single molecule method, Fiber-Seq, to explore the extreme transcriptional dynamics seen in the early Drosophila embryo ( Fig. 1A ) . Our eforts centered on 2-4 hr embryos (nuclear cycle, NC, 13, 14, and early gastrulation) since this is the time when the fate map of the adult fly is established through the regulation of over 100 patterning genes by several hundred developmental enhancers ( 10 , 11 , 12 ). We obtained ∼4000-fold coverage of the entire genome with single molecule read lengths averaging ∼23 kb. This coverage is among the highest of single molecule genomic datasets ever produced. Download figure Open in new tab Figure 1: Transcription at highly expressed genes is driven by Pol II convoys A. Schematic of Fiber-seq and FiberHMM. B. Example Fiber-seq reads at the gt locus. Reads are split between reads with (top) single and (bottom) convoys of Pol II footprints. C. Seven Fiber-seq reads at the prd locus showing Pol II convoys. D. Six Fiber-seq reads at the tll locus showing Pol II convoys. E. Histogram of the basepair distance between convoy Pol II footprints. A pair of secondary x-axes above estimate time from elongation rates. F. Pol II footprints within diferent size convoys, aligned to the 5’ end of the first Pol II footprint. Secondary x-axes above estimate time from elongation rates. Plots show 500 sampled fibers. G. Histogram of percent reads with a Pol II convoy across all expressed genes (n=4329). FiberHMM was used to visualize the chromatin-bound proteome of every individual fiber ( 7 ). We incorporated previously published PRO-Seq assays ( 13 , 14 ) to facilitate the identification of different transcriptional states, including Pol II pausing and PIC occupancy ( 7 , 15 , 16 ) (Fig. S1A,B) , as well as Pol II footprints that extend beyond the 3’ UTR (Fig. S1C) . Numerous Pol II footprints were also identified within gene bodies correlating with PRO-seq elongation signals ( Fig. 1B ). There is a strong correspondence between the fraction of fibers showing Pol II occupancy at a given gene and the overall fraction of cells expressing that gene in scRNA-seq from NC14 embryos (Figure S1D-H) ( 17 ). Pol II convoys at highly transcribed genes The patterns of elongating Pol II footprints revealed an unexpected mode of transcription enriched at highly transcribed developmental patterning genes. Chromatin fibers exhibiting low or moderate levels of transcription possess individual elongating Pol II footprints sporadically distributed across the gene body. By contrast, strongly transcribed fibers contain “convoys” of tightly spaced Pol II footprints ( Fig. 1B-D , Fig. S1I-K). Within these convoys, polymerases were spaced an average of 60 bp (center-to-center), leaving only 5-10 bp between footprints, far tighter than typical internucleosomal spacing ( Fig. 1E , Fig. S1J). This minimal spacing corresponds to the release of Pol II every 2-3 seconds ( 18 , 19 ), approaching the theoretical limit imposed by polymerase sterics ( Fig. 1E ). These release rates are consistent with the large Pol I convoys seen at rRNA genes (Fig. S1M ; 20) and previous studies of Pol II release ( 21 , 19 ). Convoys most frequently consist of 2-4 polymerases but rarely extend up to 14, representing continuous Pol II release events for upwards of 30 seconds ( Fig. 1F ). Strikingly, while these convoys are prevalent at highly transcribed developmental patterning genes, they are also observed for an average of ∼7% of all chromatin fibers at “typical” genes, suggesting a general mechanism for the release of small Pol II clusters during intense periods of transcription ( Fig. 1G , Fig. S1N) . Hyperbursts evict nucleosomes Individual Pol II footprints on fibers with low to moderate transcription exhibit a regular spacing periodicity of ∼147 bp, being found within internucleosomal linker regions (Fig. S2A). By contrast, Pol II convoys found at highly transcribed genes appear to transiently displace nucleosomes in their path ( Fig. 1B , S2A) , with a single nucleosome evicted for every ∼2.5 Pol II footprints within a convoy (Fig. S2B). Most dramatically, 0.5% to 7% of the fibers comprising 285 highly active genes display “hyperbursting”, a state of intense transcription whereby chromatin fibers lack at least 50% of the nucleosomes normally distributed across the gene body ( Fig. 2A-C ). Compared to fibers without nucleosome eviction, hyperburst fibers contain a 3-fold average increase in the number of Pol II convoys consisting of 3 or more polymerases and a 5-fold average increase in the number of fibers containing at least two Pol II convoys within the gene body ( Fig. 2D-F , Fig. S2C-E). Developmental patterning genes are strongly enriched among the genes with the highest level of hyperbursting, including all 5 genes with the highest fraction of hyperburst reads: ftz, sna, gt, tll, and eve . Consistent with a transcription-driven mechanism for nucleosome eviction we note that nucleosome densities are lower upstream of Pol II convoys than downstream ( Fig. 2G ). Download figure Open in new tab Figure 2: Hyperbursts evict nucleosomes from gene bodies A. Fiber-seq reads at zen with >50% nucleosome depletion split by active (Top) and refractory (Bottom) fibers. B. Fiber-seq reads at eve with >50% nucleosome depletion split by active (Top) and refractory (Bottom) fibers. C. Fiber-seq reads at slp1 with >50% nucleosome depletion split by active (Top) and refractory (Bottom) fibers. D. Percentage of Pol II footprints in diferent convoy sizes in non- and hyperburst fibers. E. Percentage of larger than 3 Pol II convoys in non- and hyperburst fibers. F. Percentage of fibers with more than 1 Pol II convoy in non- and hyperburst fibers. G. Nucleosome occupancy 300 bp upstream (excluding promoter) and downstream of convoys for the 100 most hyperburst-enriched genes. Error bars, 95% CI from 1000× bootstrap; ***p < 10⁻⁴ (rank-sum) H. Percent accessibility of the 10kb upstream of TSSs on non- and hyperburst fibers. Errorbars indicate 95% confidence interval from 1000x bootstrapped resamplings. Significance is indicated between bars; ***: p < 10 -4 , * p .05 (rank-sum). I. Percent PIC occupancy of non- and hyperburst fibers. Errorbars indicate 95% confidence interval from 1000x bootstrapped resamplings. Significance is indicated between bars; ***: p < 10 -4 (rank-sum). J. Scatterplot of on-rate (% reads with ≥1 convoy) of fibers binned by centiles of percent nucleosome eviction at slp1 . Grey points are individual active or inactive fibers. K. Percentage of hyperburst reads for ten patterning genes with the highest levels of hyperbursting. L. Mean Pol II footprints per read for 100 most hyperburst-enriched genes, binned by nucleosome eviction; ***p < 10⁻⁴, **p 0.05 (Mann-Whitney U). M. On-rate for 100 most hyperburst-enriched genes, binned by nucleosome eviction; ***p < 10⁻⁴, **p 0.05 (Mann-Whitney U). A subset of fibers exhibiting nucleosome depletion lack Pol II footprints ( Fig. 2A-C ) . We believe that these “refractory fibers” represent off periods between transcriptional bursts. Several arguments support this view. In addition to exhibiting extensive nucleosome depletion, putative refractory fibers possess upstream cis-regulatory element (CRE) accessibility patterns that are far more similar to those seen for hyperburst fibers than weakly transcribed fibers ( Fig. 2H , Fig. S2F). Further, these refractory fibers show a notable depletion of PIC footprints, which is consistent with a transcriptional “off” state ( Fig. 2I ) . The ratio of active to refractory hyperburst fibers provided an opportunity to capture the transcriptional kinetics of intensely active developmental genes during periods of maximum transcription ( Fig. 2J , S2G). Across all fibers, Fiber-seq derived estimates of on-rates are similar to those seen across all nuclei for live-imaging of ftz transcription ( 5 ), with estimates of 14.4% and 16%, respectively. While the fraction of fibers in the hyperburst state varied widely from gene to gene ( Fig. 2K ) , the on-rates and amplitude of transcription plateaued at high levels of nucleosome eviction ( Fig. 2L,M ) . At >80% eviction hyperbursting genes reached an average on-rates of ∼68%, whereas weak or moderately active fibers exhibiting less than 40% eviction possess on-rates of just ∼5% ( Fig. 2L,M ) . This constancy of hyperactive on-rates suggest that developmental patterning genes reach comparable states of maximal transcription and is consistent with the two-state model for transcriptional bursting based on MS2/PP7 live-imaging methods ( 22 ). Dynamic CREs for different states of transcription Hyperburst fibers exhibit enhanced accessibility of upstream regulatory regions as compared with those fibers displaying moderate or low levels of transcription, suggesting coordinated remodeling of regulatory landscapes to achieve maximum transcription ( Fig. 2G ) . To test if hyperbursting harnesses unique regulatory mechanisms, we employed Fiber-Seq Inference of Regulatory Elements (FIRE) to identify putative cis regulatory elements (CREs) within the first 10 kb upstream of all TSSs ( 23 ). Known and putative CREs exhibit highly variable accessibility from fiber to fiber ( Fig. 3A ) , illustrating a dynamic regulatory landscape driven by extensive chromatin remodeling. Download figure Open in new tab Figure 3: Hyperburst-associated regulatory elements function cooperatively at sna A. Percent accessibility for FIRE-identified CREs <9kb upstream of an expressed gene. B. Maximum ΔAUC for an individual CRE upstream of each expressed gene. C. Simplex plot showing regression coeficients predicting transcriptional states (vertices) from CRE accessibility (points); colors indicate strongest coeficient. D. Fiber-seq reads at sna colored by local accessibility density. CREs schematized above; red boxes, active hyperburst-associated CREs. Right, same reads with Pol II footprints in white. E. Fiber-seq at sna grouped by number of accessible hyperburst-associated CREs (count indicated on left); 50 sampled reads shown for 0–1 active elements. F. Count of elongating Pol II footprints per read for combinations of hyperburst-associated CREs, sorted by mean count. Combinations of elements are schematized below, with black boxes showing accessibility. Bars are colored by accessible element counts. Errorbars represent the 95% confidence interval from 1000x bootstrapped resamplings of the fibers. G. Synergy scores from logistic-dropout models showing cooperative efects of sna hyperburst-associated CRE combinations (diference between sum of single-dropouts and combination dropout on model AUC). H. Synergy scores from logistic-dropout models showing cooperative efects of CRE combinations (diference between sum of single-dropouts and combination dropout on model AUC). We then tested the relationship between CRE accessibility and transcription by directly associating dynamic accessibility patterns with different transcriptional states on the same individual chromatin fibers, training logistic regression models with activities of CREs as predictors and transcriptional state assignments as responses. Transcriptional states were chosen to capture distinct stages of gene activation: poised (accessible promoter and paused Pol II), active transcription (non-hyperburst elongation and hyperburst), and attenuated (Pol II footprints only within the 3’ termination region). Overall, 73% of genes had at least one CRE that was predictive of a transcriptional state, and for 68% of genes, the overall pattern of CRE activity was predictive of at least two different transcriptional states ( Fig. 3B , Fig. S3A) . These observations support specialized roles for individual CREs in different phases of the transcription cycle ( Fig. 3C , Fig. S3B,C). Other than poised promoters, combined CRE activities were particularly predictive of hyperbursts (Fig. S3D-I) . Hyperburst regulatory elements at the snail locus To explore the regulatory basis for hyperbursts we turned to snail ( sna ), which encodes a key determinant of epithelial-mesenchyme transitions (EMT) including mesoderm invagination at the onset of gastrulation ( 24 , 25 ). Classical studies identified two developmental enhancers, a distal, “shadow enhancer” (herein D) which is the primary activator, and a proximal enhancer (herein P) that has a weak activating effect and a potential role in attenuation ( 26 , 27 , 28 ). FIRE identified seven distinct modules in the sna 5’ regulatory region, including three modules in the D enhancer and another two in P. A third, bipartite region was also identified in central regions (C1 and C2) that was not previously implicated in sna regulation. All seven modules appear to possess specialized functions with respect to states of transcription: hyperburst-associated (D1, D2, C2, P2), pause-associated (D3, C1), and attenuation-associated (P1) elements ( Fig. 3C,D ). This combinatorial logic suggests inter-connected regulation that transcends independent action of separate enhancer modules. These putative redundancies could be explained either by overlapping activities or distinct contributions of different modules. Consistent with the former model, the burst-associated modules (D1, D2, C2, P2) frequently showed co-accessibility: 30% of fibers had one module active while 17% had two or more ( Fig. 3E ) . Pol II footprint counts and gene-body accessibility scaled with the number of active modules ( Fig. 3F , Fig. S3J) . This scaling was more than additive, with the overall contributions of multiple elements (in particular, D1, D2, and C2) being higher than expected when dropped out of regression training, supporting a model wherein CREs cooperatively tune transcriptional output, rather than functioning redundantly and separately ( Fig. 3G ) . This cooperativity is broadly detected, with hyperburst-associated elements generally showing greater synergy as compared with CREs associated with pausing, elongation, or termination ( Fig. 3H , S3K) . Developmental trajectories at the ftz locus We next sought to determine the role of hyperbursting within the developmental trajectory of a key pair-rule gene, ftz , which exhibits remarkable transcriptional dynamics during NC14 ( 29 ). Live-imaging assays reveal activation, bursting, attenuation, and silencing within a span of just one hour ( 5 ). Classical gene fusion studies identified two developmental enhancers for ftz regulation, a proximal Zebra element (Z) located adjacent to the promoter and a distal autoregulatory element (UPS) located ∼5 kb upstream ( 30 , 31 , 32 , 33 ). Z initiates ftz expression in early embryos, while the UPS maintains and refines these stripes in older embryos ( 29 ). As seen for sna , ftz hyperburst fibers contain maximally accessible CRE modules, in this case the Z enhancer and UPS ( Fig. 4A ). Download figure Open in new tab Figure 4: Hyperburst-associated regulatory elements attenuate ftz expression A. (Top) Hyperburst fiber-seq reads at the ftz locus. (Bottom) The same Fiber-seq reads, ordered identically but (Left) at the UPS and (Right) at the Zebra element. B. Density of nuclei expressing ftz sense intronic RNA and ftz R asRNA across the anterior-posterior axis measured via smFISH (R-element: n = 356, intron: n = 663). C. Fiber-seq reads overlapping the ftz locus and regulatory region, divided into steps based on the top-scoring trajectory of transcriptional activation and attenuation. Reads with no accessible elements are sampled to 150 total. D. Enrichment of activator and repressor motifs at the UPS, with each CRE module indicated via colored blocks. Motifs are identified with a permissive log-odds threshold (≥ 2.0) and filtered by overlapping subnucleosomal Fiber-seq footprints in at least 5% of total fibers. E. Boxplot showing the ratio of activator and repressor motif enrichment for non- and hyperburst-associated CREs. Motifs are identified with a permissive log-odds threshold (≥ 2.0) and filtered by overlapping subnucleosomal Fiber-seq footprints in at least 5% of total fibers. FIRE revealed unexpected complexity within the UPS, having identified three distinct modules (U1, U2, U3), each associated with distinct transcriptional states (Fig. S4A) . U1 is associated with poised, elongating, and hyperbursting states, whereas U2 is an extension of U1 and is associated primarily with hyperbursting. U3 is associated with both hyperbursting and attenuation (Fig. S4A). This analysis also identified a novel CRE, “R”, that is associated both with hyperbursting and attenuation ( Fig. 4A , S4A,B) . In support of a role as an attenuator, the R element produces an antisense RNA that is expressed in interstripe regions where ftz is repressed during stripe refinement ( Fig. 4B , S4C,D) . R is also enriched for autoregulatory and repressive binding motifs (Fig. S4E). Only 17 of the 32 possible accessibility combinations of these elements were observed in more than 1% of fibers, indicating a constrained regulatory logic. The arrangement of ftz fibers that we present is based on the highest scoring route from a graph-based trajectory analysis, where nodes represent the observed combinations of regulatory states and edges represent single-element transitions prioritizing those supported by the maximum number of fibers ( Fig. 4C , Fig. S4F,G). The primary activation trajectory proceeded through six ordered stages: (i) Promoter opening (102 fibers), the promoter becomes accessible while enhancers remain closed; (ii) Early activation (110 fibers), Zebra or U1 become accessible and initiate low levels of transcription; (iii) Burst amplification (62 fibers), U2 and U3 modules open and transcription intensifies; (iv) Hyperbursting (46 fibers), R element opens during maximal transcription; (v) Activation shutdown (67 fibers), Zebra/U1/U2 close while low levels of transcription continue; and (vi) Active repression (230 fibers), promoter closes but U3/R remain accessible ( Fig. 4C , Fig. S4F,G). This trajectory aligns with known developmental timing, with the Z element initially activating low and moderate transcription, followed by UPS activation via autoregulation leading to sustained, bursty transcription ( 29 ). As observed for sna , hyperbursting is associated with widespread accessibility of the regulatory landscape, particularly the specialized bursting elements (e.g., U2 and U3). This extensive remodeling of the regulatory landscape and ftz gene body exposes repressor motifs (e.g., Hb, Slp1, Run, and Hairy) across the U2, U3 and R elements ( Fig. 4E , Fig. S4E) . This observation is consistent with a role for hyperbursting in sensitizing ftz to repressive regulatory inputs. Globally, the ratio of repressive and activating motifs within hyperburst-associated CREs is more balanced than the motif distributions seen for CREs associated with other states such as pausing and termination. These observations suggest that bifunctional activation and attenuation may be widespread ( Fig. 4E ; see below) . Discussion We have presented evidence that many or most developmental control genes display hyperbursting, a short-lived peak of transcriptional activity coinciding with extensive eviction of nucleosomes across gene bodies. A similar mass eviction was previously reported for hsp70 upon heat shock in Drosophila , where denuded chromatin fibers facilitate peak rates of transcription ( 34 ). The developmental control genes we characterized appear to achieve similar peaks of activity, consistent with previous live-imaging studies suggesting similar on-rates for Drosophila gap genes ( 35 ). This suggests that eviction of nucleosomes may be a general mechanism to achieve maximized transcription, and we expect that this state will also be observed in mammalian systems requiring extreme transcriptional outputs like erythropoiesis and somitogenesis. Our analyses connecting CRE accessibility to transcriptional states suggest that classically defined developmental enhancers are often composed of multiple elements with distinct functions. Many of these elements are only accessible on fewer than 10% of chromatin fibers, yet their functional importance is supported by strong associations with specific transcriptional states. This observation supports an outsize role for rapid chromatin remodeling in mediating developmental enhancer switching and fine-tuning gene expression. For example, the full-length distal sna (aka shadow) enhancer appears to contain three distinct cis-regulatory elements (D1, D2, and D3) with differing activities. D3 accessibility coincides with paused Pol II at the sna promoter, whereas D1 and D2 accessibility coincides with Pol II elongation and hyperbursts, respectively. These distinct activities could explain why the distal sna enhancer directs especially robust expression of reporter genes in transgenic embryos ( 36 ). An interesting implication of this analysis is the apparent necessity of long-range cooperation of numerous regulatory elements to achieve hyperbursting of developmental control genes. Hyperbursting of sna does not depend solely on the D1 and D2 elements within the distal enhancer but also appears to require maximal accessibility of the C2 and P2 elements located 3 kb and 6 kb downstream of D2, respectively. Similarly, hyperbusting of ftz coincides with peak accessibility of all three elements comprising the autoregulatory UPS enhancer as well as the proximal zebra enhancer, separated by over 4 kb. The extensive regulatory accessibility associated with hyperbursting suggests two mechanisms for triggering auto-attenuation. First, many hyperburst-associated regulatory modules comprising developmental enhancers appear to possess dual activator and repressor activities. For example, the U1 element of the ftz UPS autoregulatory enhancer achieves maximum accessibility, along with the proximal zebra enhancer, as transcription increases. U1 contains known activator binding motifs (e.g., Bicoid and Caudal), but these are not sufficient for high levels of transcription. Hyperbursts are observed only when the neighboring U2 and U3 modules also become accessible, possibly through spreading of U1 activation. U2 and U3 are enriched for both autoregulatory activator (e.g., Ftz) and repressor binding motifs (Hairy and Runt), supporting a function both in enhancing hyperbursts and subsequent attenuation ( Fig. 4D ). The second potential mechanism of attenuation is more direct. Namely, the ftz R element appears to become accessible only upon hyperbursting-driven nucleosome eviction in the 3’ UTR (Fig. S4E). The R element also appears to mediate repression at later stages of development when ftz is regulated by a new set of enhancers during neurogenesis (Fig. S4H) (Heffer et al., 2013). It is unclear whether antisense transcription arising from the R element directly contributes to the attenuation of ftz expression in interstripe regions. It is possible that antisense transcription precludes transcription of sense RNAs. Alternatively, antisense RNAs might work post-transcriptionally to degrade lingering ftz mRNAs produced in interstripe regions at earlier stages of development when ftz is preferentially transcribed in future interstripes ( Fig. S1I ) ( 5 ). Our Fiber-seq analysis provides a glimpse into the inter-connected activities of the complex regulatory landscapes controlling gene activity during development. This analysis suggests a level of regulatory integration and cooperativity not envisioned by the traditional modular view of the regulatory genome, whereby separate enhancers function in an autonomous manner to produce composite patterns of gene expression. Overall, enhancers and their constituent modules are rarely active as individual entities, but rather, many transcriptional states (most notably hyperbursts) appear to require the cooperation of multiple modules. Recent studies provide some evidence for regulatory integration, including examples of long-range enhancer-promoter compatibility across neighboring TADs ( 37 , 38 ). We expect that further investigation of regulatory networks at the single-fiber scale will reveal that this cooperativity is a general feature of the gene activities controlling dynamic developmental processes across the metazoa. Author contributions Conceptualization: TWT, ML Investigation: TWT, TB Validation: TWT, TB Formal analysis: TWT, TB Visualization: TWT, TB Funding acquisition: ML Supervision: ML Writing - original draft: TWT, ML Writing - review and editing: TWT, TB, ML Declaration of Interests TWT is an inventor on a provisional patent application related to Fiber-seq applied to human sperm chromatin. This application is unrelated to the current study. Download figure Open in new tab Figure S1: Quantification of transcriptional states A. (Left) example Fiber-seq reads at the gt locus showing promoter-proximal paused Pol II footprints. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). (Right) min-max normalized enrichment of Pol II footprints at the promoter proximal paused sites for genes binned into deciles of PRO-seq coverage within the promoter region. Error shading indicates 95% confidence interval from 1000x bootstrapped resamplings. B. (Left) example Fiber-seq reads at the gt locus showing PIC footprints. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). (Right) min-max normalized enrichment of PIC footprints a the TSS for genes binned into deciles of PRO-seq coverage within the gene body region. Error shading indicates 95% confidence interval from 1000x bootstrapped resamplings. C. (Left) example Fiber-seq reads at the gt locus showing terminating Pol II footprints. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). (Right) min-max normalized enrichment of Pol II footprints within the termination region (polyA site to polyA site +500bp) for genes binned into deciles of PRO-seq coverage within the termination region. Error shading indicates 95% confidence interval from 1000x bootstrapped resamplings. D. (Top) example in situ from the Berkeley Drosophila Genome Project showing spatial expression of cnc . (Bottom) corresponding Fiber-seq tracks showing PRO-seq, a metaprofile of nucleosome and subnucleosome footprints, and 100 randomly sampled Fiber-seq reads at the locus. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). E. (Top) example in situ from the Berkeley Drosophila Genome Project showing spatial expression of gt . (Bottom) corresponding Fiber-seq tracks showing PRO-seq, a metaprofile of nucleosome and subnucleosome footprints, and 100 randomly sampled Fiber-seq reads at the locus. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). F. (Top) example in situ from the Berkeley Drosophila Genome Project showing spatial expression of btd . (Bottom) corresponding Fiber-seq tracks showing PRO-seq, a metaprofile of nucleosome and subnucleosome footprints, and 100 randomly sampled Fiber-seq reads at the locus. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). G. Scatterplot showing correlation of fraction nuclei expressing a gene and fraction accessible chromatin fibers for a subset of 63 highly expressed developmental patterning genes. Histograms on each axis indicate distributions. Pearson correlation is indicated on plot. H. Scatterplot showing correlation of fraction nuclei expressing a gene and fraction chromatin fibers showing any evidence of transcription (paused Pol II, PIC, elongating Pol II) in Fiber-seq for a subset of 63 highly expressed developmental patterning genes. Histograms on each axis indicate distributions. Pearson correlation is indicated on plot. I. (Top) mean convoy count and (Bottom) mean convoy size for convoys called using a range of inter-polymerase distance cutofs. The chosen cutof (100bp) is indicated with a dashed line. J. Boxplot comparing distribution of inter-nucleosome and inter-polymerase distances. K. Min-max normalized enrichment of single or convoy elongating Pol II footprints within gene bodies for genes binned into deciles of PRO-seq coverage within the gene body region. L. Example Fiber-seq reads from human lymphoblastoid cell lines from the donor for the HG002 assembly within rDNA repeats. Footprints are labeled by predicted identity (green, nucleosome; purple to blue, putative Pol I of 40-80bp size range). M. Histogram of the ratio of single to convoy Pol II footprints across all genes showing expression in scRNA-seq in NC14. Download figure Open in new tab Figure S2: Extended quantification of hyperburst kinetics A. Histogram of spacing between elongating Pol II footprints split between single Pol II footprints and convoy Pol II footprints. N. Lineplot showing count of Pol II footprints within convoys and overall count of nucleosomes within the first 1.5kb of a gene. Error shading indicates 95% confidence interval from 1000x bootstrapped resamplings. B. Barplot showing fraction of polymerases in diferent convoy sizes, split into 10 bins of nucleosome eviction rate. C. Barplot showing the percentage of larger than 2 Pol II convoys in non- and hyperburst fibers. The percent diference is indicated on the plot. D. Barplot showing the mean count of Pol II convoys per read in non- and hyperburst fibers. The percent diference is indicated on the plot. E. Barplot showing the Pearson correlation of the 9kb upstream pattern of nucleosome occupancy for active hyperburst fibers to refractory or to non-hyperburst fibers. F. Scatterplots for the top twelve developmental patterning genes by hyperburst activity. Points represent the percentage of fibers with a Pol II footprints binned by centiles of percent nucleosome eviction. Download figure Open in new tab Figure S3: Prediction of regulatory element function via logistic regression A. Histogram of the count of transcriptional states predicted (ΔAUC > .05) by the combination of CRE activities observed at a given gene for regression models trained on those combinations on a per-gene basis. B. Simplex plot showing the balance of logistic regression coeficients for the prediction of five transcriptional states from CRE accessibility for all combinations of two CREs identified at genes with hyperbursting (n= 1250), calculated per-gene. Points are colored based on their strongest correlation coeficient. C. Simplex plot showing the balance of logistic regression coeficients for the prediction of five transcriptional states from CRE accessibility for all combinations of three CREs identified at genes with hyperbursting (n= 1250), calculated per-gene. Points are colored based on their strongest correlation coeficient. D. Fraction of genes where each individual transcriptional state was predictable (ΔAUC > .05) by the combination of CRE activities observed at a given gene for regression models trained on those combinations on a per-gene basis. E. (Top) Histogram of per-gene ΔAUC for promoter accessibility. Median is indicated on the plot. (Bottom) Scatter plot of per-CRE AUC and ΔAUC for promoter accessibility. Density is indicated via colormap. F. (Top) Histogram of per-gene ΔAUC for pausing. Median is indicated on the plot. (Bottom) Scatter plot of per-CRE AUC and ΔAUC for pausing. Density is indicated via colormap. G. (Top) Histogram of per-gene ΔAUC for non-hyperburst elongation. Median is indicated on the plot. (Bottom) Scatter plot of per-CRE AUC and ΔAUC for non-hyperburst elongation. Density is indicated via colormap. H. (Top) Histogram of per-gene ΔAUC for hyperburst. Median is indicated on the plot. (Bottom) Scatter plot of per-CRE AUC and ΔAUC for hyperburst. Density is indicated via colormap. I. (Top) Histogram of per-gene ΔAUC for termination. Median is indicated on the plot. (Bottom) Scatter plot of per-CRE AUC and ΔAUC for termination. Density is indicated via colormap. J. Percent gene-body accessibility for combinations of hyperburst-associated CREs, sorted identically to Fig. 3E . Combinations of elements are schematized below, with black boxes showing accessibility. Bars are colored by count of elements accessible. Errorbars represent the 95% confidence interval from 10,000x bootstrapped resamplings of the fibers. K. Barplot showing the percentage of genes with more than one CRE primarily predictive of hyperbursting, defined as their maximum regression coeficient being the hyperburst coeficient. Download figure Open in new tab Figure S4: Reconstruction of the regulatory trajectory of ftz from Fiber-seq reads A. Heatmap showing logistic regression coeficients for the five CREs identified by FIRE. Coeficients are listed within cells of the heatmap. B. 100 sampled footprinted Fiber-seq reads at the ftz locus where the R element is accessible. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). C. PRO-seq signal within the ftz gene. Positive track indicates sense transcription, negative track indicates antisense transcription. The region covered by the R element is schematized on the plot. D. Representative smFISH image of R-element asRNA localization. Arrows are added to highlight nuclei showing R-element asRNA activity. E. Enrichment of motifs for activators (top) and repressors (bottom) at the R element, with each submodule indicated via colored overlays. Motifs are filtered based on a permissive log-odds threshold (≥ 2.0), and support of overlapping subnucleosomal footprints in Fiber-seq in at least 5% of total fibers. F. Graph showing the top route (based on overall count of reads included) for trajectory of transcription activation and attenuation of the ftz gene based on Fiber-seq. Each node represents either opening (top 5 steps) or closing (bottom 5 steps) of the indicated element on the x-axis. Groups of transition steps are colored and labeled on the y-axis as described in the results. G. Graph showing all possible routes (15 total) for trajectory of transcription activation and attenuation of the ftz gene based on Fiber-seq. Each node represents either opening (top 5 steps) or closing (bottom 5 steps) of the indicated element on the x-axis. Groups of transition steps are colored and labeled on the y-axis as described in the results. Edges are thickened based on the count of reads across all trajectories supporting the step. H. Active footprinted Fiber-seq reads at the ftz locus during the 6-8hr period. Footprints are labeled by predicted identity. Footprints are labeled by predicted identity (green, nucleosome; pink, Pol II; grey, unlabeled). The neurogenic enhancer and R element are indicated. Reads are sorted to match the 2-4hr trajectory. View this table: View inline View popup Download powerpoint Table S1: Probe sequences used for smFISH Methods Fiber-seq 2-4hr and 6-8hr embryos were collected and dechorionated for 1 minute 45 seconds in 50% bleach. Embryos were then resuspended in homogenization bufer (250mM Sucrose,10mM Tris-CL, pH 8.0, 25mM KCl, 5mM MgCl 2 , 0.1 mM, 0.1% Triton X-100) and dounced 40 times using a tight dounce to isolate nuclei. Nuclei were filtered via a 40um filter and spun down at 350 × g for at 4°C for 10 minutes. Pellets were resuspended in 30 μL Bufer A (15 mM Tris-Cl pH 8.0, 15 mM NaCl, 60mM KCl, 1 mM EDTA, 0.5 mM EGTA, 0.5 mM Spermidine). 0.8 mM SAM and 200 U of Hia5 were then added, followed by incubation for 10 minutes at 25°C. The reaction was quenched with the addition of 1% final SDS and mixing using a wide-bore pipette tip. Samples for several collections were pooled, and their DNA was purified using the Promega Wizard HMW DNA Extraction Kit (Promega A2920) and submitted for library preparation and sequencing. A total of 24 collections for each timepoint were pooled, barcoded, and sequenced on the Pacific Biosciences Revio. 11 total sequencing runs were carried out to achieve around ∼4000x coverage of the genome. FiberHMM Footprints were called on Fiber-seq reads using FiberHMM. FiberHMM is based on a hidden Markov model (HMM) with two hidden states--accessible and inaccessible. To account for sequence-related biases from Hia5 or the methylation caller, at each position the model takes into account both the base and its surrounding 6 bp sequence context. The emission probabilities used in the model are the probabilities of methylation of a given base with its +/− 3 bp sequence context in an accessible or inaccessible state based on experimental-derived methylation rates from control datasets. The probability of methylation given an accessible state was based on the methylation frequency in a dataset generated from dechromatinized Drosophila S2 cell genomic DNA. The probability of methylation given an inaccessible state was based on the methylation frequency in a dataset generated from Drosophila S2 cells untreated with Hia5. Transition and starting probabilities for the HMM were carried over from a model trained on Drosophila S2 genomic DNA with initial probabilities picked from the Dirichlet distribution with all parameters set to 1 ( 7 ). Defining accessibility and footprints Nucleosome footprints were in general defined as footprints larger than 90bp. Pol II footprints were defined as footprints between 40 and 60 bp, overlapping signal from PRO-seq. PIC footprints were defined as either 20-40bp overlapping a TATA box (TBP), or 60-80 bp overlapping the TSS, as defined via CAGE-seq ( 40 ). Transcription factor footprints were defined as <90bp, not overlapping PRO-seq signal. To account for transcription factor footprints in close proximity leading to merged footprints, footprints larger than 90bp and less than 200bp were called as transcription factors if there was support from at least 50 reads with subnucleosomal footprints that perfectly aligned with both ends of the larger footprint. Accessible regions in general were defined as all gaps between nucleosomal footprints, ignoring called transcription factor and polymerase footprints. Comparing Fiber-seq to PRO-seq, scRNA-seq, and in situ hybridizations Processed PRO-seq and scRNA-seq datasets from nuclear cycle 14 were retrieved from GSE211220 and DVEX respectively. For PRO-seq, genes were binned into equal deciles of signal within the pause site (TSS to +200bp), body (+200bp to end), or termination region (end + 500bp). Mean fiber-seq footprint occupancy (pause, PIC, elongating, terminating) within the same regions was then calculated for each gene, with the pooled means plotted for each decile. For comparisons to scRNA-seq, per-cell matrices were accessed and the count of genes with at least one cell showing one transcript was calculated. This was then compared to the fraction of reads showing an accessible promoter or the fraction of reads showing any Pol II footprint within any of the above regions was plotted. Gene expression patterns for comparison to accessibility in Figure S1 were retrieved from the BDGP in situ hybridization database (Berkeley Drosophila Genome Project, https://insitu.fruitfly.org ; accessed 27 Oct 2025). Subsets of genes For analyses involving all genes, active genes from scRNA-seq as defined above (n=4329) were used in subsequent analyses. For analyses involving a subset of developmental genes, this curated list of 63 NC14 developmental patterning genes were used: bcd, cad, hb, Kr, kni, gt, tll, hkb, eve, run, ftz, odd, prd, slp1, slp2, opa, en, inv, wg, hh, ptc, smo, ci, dl, twi, sna, rho, sim, vnd, ind, msh, tor, trk, tsl, Egfr, spi, vn, dpp, scw, sog, tld, tsg, tkv, sax, put, Mad, Medea, brk, shn, pan, zen, lab, pb, Dfd, Scr, Antp, Ubx, abd-A, and abd-B. Defining convoys Convoys were defined as groups of two or more Pol II footprints within a gene body separated by 100bp or less (approximately two Pol II footprints maximum). This definition was picked based on calibration curves shown in Figure S1I, which found that a maximum 100bp separation captured 95% of grouped polymerases. Transcription state definitions Transcription states were defined as follows for regression analyses. An accessible promoter was defined as at least 100bp of accessibility on average surrounding the TSS of a gene on a read. A paused Pol II was defined as mean footprint size of 40-60 bp within the first 50bp downstream of the TSS. Elongating Pol II was defined as at least one called Pol II footprint found within the gene body. The hyperburst state was defined as at least 50% accessibility of the gene body. Defining on-rate, amplitude On-rate was estimated as the fraction of fibers with at least one convoy of Pol II. Relative burst amplitude was estimated as the mean count of polymerases found on a fiber. Defining hyperbursts, active and refractory Globally, hyperbursts were defined as fibers with less than 50% nucleosome occupancy within the first 1.5kb of the gene body (or less, if the gene was shorter, with all derivative calculations accounting for the shorter length). Refractory states were defined as fibers with that level of nucleosome eviction, but with no convoys of polymerase, while active states were defined as the population of fibers with convoys. Comparison of on-rates between Fiber-seq and live-imaging On-rate for live imaging data for ftz was calculated using 181 background-normalized fluorescence measurements for 1541 cells taken at 15s intervals ( 5 ). A gaussian mixture model was applied to identify a threshold for a 99% confident active state, and the total count of active and inactive states were counted. For Fiber-seq, the fraction of reads at the ftz locus showing convoy activity was calculated and compared to the live-imaging result. Fiber-seq Inference of Regulatory Elements (FIRE) FIRE was applied to the datasets using fibertools v0.7.0 ( 41 ). Aggregate traces of FIRE signal were then generated, with peaks called using scipy peak calling (minimum width = 20 bp, minimum height = 1/5 th of the maximum peak in the window). Peaks less than 100bp were expanded to 100bp, and overlapping peaks were merged. Accessibility of elements within these regions called based on footprinting as described above. Transcription state definitions Transcription states were defined as follows for regression analyses. An accessible promoter was defined as at least 100bp of accessibility on average surrounding the TSS of a gene on a read with no polymerase footprints overlapping the gene. A paused Pol II was defined as a mean footprint size of 40-60 bp within the first 50bp downstream of the TSS with no other polymerase footprints. Elongating Pol II was defined as at least one Pol II footprint found within the gene body. The hyperburst state was defined as at least 50% accessibility of the gene body. Termination was defined as at least one Pol II footprint found within the 500bp after the annotated gene end. Logistic Regression Modeling of Transcriptional States The probability that an individual chromatin fiber belonged to a specific transcriptional state—accessible promoter, paused, elongating, hyperbursting, or terminating—was modeled as a function of the accessibility of nearby cis-regulatory elements (CREs) identified by Fiber-seq Inference of Regulatory Elements (FIRE). For each gene, all putative CREs located within a 9 kb window upstream of the annotated transcription start site (TSS) were identified using FIRE. Accessibility was encoded on a per-read basis. The resulting binary feature matrix (reads × features) was constructed for each gene. Columns corresponding to transcriptional states (accessible, pause, elongating, hyperburst, terminating) served as response variables, and all remaining binary features were used as predictors. Each transcriptional state was modeled independently for every gene using penalized logistic regression of the form: where yᵢ denotes the presence of the transcriptional state on read i and xᵢⱼ denotes the accessibility of feature j on that read. Models were fit using the liblinear solver in scikit-learn (LogisticRegression) with an L1 penalty (C = 1.0, class_weight = ’balanced’, max_iter = 500). L1 regularization encouraged sparsity and yielded interpretable coeficients indicating the direction and magnitude of association between each CRE and the transcriptional state. In addition to these multi-feature models , parallel single-feature regressions were fit in which each CRE was considered individually against the same response variable. The single-CRE fits provided baseline predictive contributions and facilitated later analyses of synergy between elements. To prevent overfitting, stratified K-fold cross-validation (K = 5 or the maximum feasible value given class sizes) was performed within each gene. Predicted probabilities for held-out reads were collected to compute out-of-fold (OOF) receiver operating characteristic area under the curve (ROC-AUC) and average precision (AP) scores. Because class imbalance was common, AUC was used as the principal, threshold-independent performance metric. A within-gene permutation null was generated by shufling transcriptional-state labels across reads 50 times and recomputing OOF AUCs. The diference between the observed and permutation-mean AUC (ΔAUC) quantified predictive separation beyond random expectation. For each gene and transcriptional state, the fitted L1 coeficients were aligned to the full feature set to generate a feature × state matrix of weights. Positive coeficients indicated features whose accessibility increased the probability of the corresponding transcriptional state on the same molecule; negative coeficients indicated a negative association. Coeficient magnitude reflected predictive contribution but not causal efect. Features retained after regularization were used in subsequent analyses of feature importance and regulatory synergy by comparing model performance (AUC) after selective feature removal. CRE Dropout analyses Feature dropout analyses were performed to quantify the contribution of individual CREs to transcriptional-state prediction. This procedure measured the change in cross-validated model performance when one predictor features was omitted. The performance change for each feature was computed as: Positive ΔAUC values indicated that a feature contributed to predictive performance, with magnitude reflecting relative importance within that gene and state. Synergy between CREs To test for cooperative efects, models were refit after simultaneous removal of feature subsets S of size k (typically k = 2–4). Synergy was quantified as: Negative synergy values indicated that the joint performance drop from removing all features in S exceeded the sum of the individual drops, implying cooperative rather than additive efects in predicting the transcriptional state. All dropout and synergy AUCs were obtained using stratified K-fold cross-validation (K = 5 or smaller when limited by class counts) to ensure that ΔAUC values reflected out-of-fold performance diferences rather than overfit training improvements. Significance was assessed by repeating the analysis on label-permuted data (N = 50 permutations per gene/state) to generate an empirical null distribution. smFISH Probe Synthesis Probes were designed computationally with a 25-nucleotide hybridization region with a melting temperature of at least 45 C, and specific 20-nucleotide forward primer and 30-nucleotide reverse primer with an incorporated T7 promoter. Probe synthesis followed a modified version of the protocol that was outlined in Mateo et al. Nature Protocols 2021 ( 39 ). Briefly, probes were ordered using an oligo library from IDT at 10 pmol/oligo concentration. Probes were amplified with PCR and then treated with a HiScribe T7 in-vitro-transcription (IVT) utilizing the incorporated T7 promoter at 37 C overnight. A specific ATTO-conjugated primer (also from IDT) was added and incorporated into the probes following a reverse transcription reaction with the Maxima H Minus RT kit at 55 C for 5 hours. Resulting ssDNA probes were treated with RNAse A for 0.5 hours at 37 C to remove leftover ssRNA and stored at 4 C in the dark. Intermediary and final products were cleaned with a Zymo RCC-25 kit and reconstituted in nuclease-free H20 following each step. smFISH Single-molecule FISH (smFISH) was performed on embryos fixed in 4% formaldehyde and stored in methanol at – 20 C. We performed smFISH utilizing ATTO565 and ATTO633 conjugated probes. Embryos were washed triple-washed in 1x PBS + 0.1% Tween 20, and then incubated in 35% formamide, 0.1% Tween-20, and 4x SSC for 1 hour at room temperature. Embryos were then hybridized with 750 ng of fluorophore-conjugated probe diluted in 5X SSC Bufer, 50% formaldehyde, 0.1% Tween-20, 10% dextran sulfate, and 100 μg/ml herring sperm DNA at 37 C overnight. The following day, embryos were triple washed in 35% formamide, 0.1% Tween-20, and 4x SSC over a period of 2.5 hours, triple washed and incubated in 1x PBS + 0.1% Tween 20 over 30 minutes and mounted in Prolong Diamond (Invitrogen) to cure overnight. smFISH imaging The embryos were imaged using a Zeiss LSM 880 confocal microscope (Zen software 2.3 SP1) with a Plan-Apochromat ×40/1.3 N.A. oil-immersion objective. Lasers at 488 nm, 561 nm, and 633 nm were used to excite the fluorophores (Alexa 488, ATTO 565, ATTO633). Volumetric field of view was typically set at 354.2 µm x 354.2 µm x 7.5 µm, with a voxel size of 86 nm x 86 nm x 250 nm. Imaging was done with 4x averaging with a pixel dwell time of 8 µs. Graph-based regulatory trajectory Accessible elements surrounding the ftz locus, including the core promoter, were defined by FIRE, and accessibility on each read was assigned based on footprinting patterns. Each read spanning the −6 kb to +2 kb region relative to the ftz TSS was represented by a binary vector of element accessibility. Combinations of accessibility present in <1% of fibers (n < 15) were discarded. All possible routes were then enumerated from initial states (promoter open alone) to peak transcriptional activity (all elements open) and to terminal states (R open alone or U3 open alone). States were not permitted to repeat, and each transition involved the change of only one CRE at a time. Routes were ranked by the number of reads supporting each path, with all paths representing 847 of 916 total active reads and the top-ranked route containing 576 reads. Motif enrichment Transcription-factor motifs were identified within regulatory elements using position-frequency matrices for developmental regulators (e.g. for ftz : bcd, cad, ftz, run, odd, prd, h, vfl) were retrieved from the JASPAR 2024 CORE database via the jaspardb interface ( 42 ). Each motif’s counts matrix was normalized with a pseudocount of 0.8 and converted to a log-odds position weight matrix (PWM). Sequences were scanned base-by-base using a permissive log-odds threshold (≥ 2.0), and per-position scores were combined across all motif variants for each TF by taking the maximum value. Motifs were further filtered to retain only those overlapping a putative transcription-factor footprint observed on ≥5% of fibers (generally, n = 100), thereby enriching for actively engaged motifs. For each annotated CRE, motif intensities were summed across the interval, producing a per-CRE, per-TF motif-enrichment matrix. For genome-wide analyses, a restricted set of activator and repressor motifs was used to enable comparisons in motif balance across loci (activators: bcd, cad, dl, vfl; repressors: hb, Kr, run, odd). Motif balance was calculated as the ratio of minimum to maximum activator and repressor motif scores, with zero representing complete imbalance and 1 representing perfect balance. Acknowledgements This study was supported by a grant from the NIH (GM118147). We thank the members of the Levine lab for discussion and feedback on this work. PacBio sequencing was carried out by the University of Washington PacBio Sequencing Core. Funder Information Declared NIH Common Fund , GM118147 References 1. ↵ Lettice LA , Williamson I , Wiltshire JH , Peluso S , Devenney PS , Hill AE , Essafi A, Hagman J, Mort R, Grimes G, DeAngelis CL, Hill RE. Opposing functions of the ETS factor family define Shh spatial expression in limb buds and underlie polydactyly . Dev Cell . 2012 Feb 14; 22 ( 2 ): 459 – 67 . doi: 10.1016/j.devcel.2011.12.010 . PMID: 22340503; PMCID: PMC3314984. OpenUrl CrossRef PubMed Web of Science 2. ↵ Jeong Y , El-Jaick K , Roessler E , Muenke M , Epstein DJ . A functional screen for sonic hedgehog regulatory elements across a 1 Mb interval identifies long-range ventral forebrain enhancers . Development . 2006 Feb; 133 ( 4 ): 761 – 72 . doi: 10.1242/dev.02239 . Epub 2006 Jan 11. PMID: 16407397. OpenUrl Abstract / FREE Full Text 3. ↵ Rodríguez-Carballo E , Lopez-Delisle L , Willemin A , Beccari L , Gitto S , Mascrez B , Duboule D . Chromatin topology and the timing of enhancer function at the HoxD locus . Proc Natl Acad Sci U S A . 2020 Dec 8; 117 ( 49 ): 31231 – 31241 . doi: 10.1073/pnas.2015083117 . Epub 2020 Nov 23. PMID: 33229569; PMCID: PMC7733857. OpenUrl Abstract / FREE Full Text 4. ↵ Dunipace L , Saunders A , Ashe HL , Stathopoulos A . Autoregulatory feedback controls sequential action of cis-regulatory modules at the brinker locus . Dev Cell . 2013 Sep 16; 26 ( 5 ): 536 – 43 . doi: 10.1016/j.devcel.2013.08.010 . PMID: 24044892; PMCID: PMC3782659. OpenUrl CrossRef PubMed Web of Science 5. ↵ Lim B , Fukaya T , Heist T , Levine M . Temporal dynamics of pair-rule stripes in living Drosophila embryos . Proc Natl Acad Sci U S A . 2018 Aug 14; 115 ( 33 ): 8376 – 8381 . doi: 10.1073/pnas.1810430115 . Epub 2018 Jul 30. PMID: 30061421; PMCID: PMC6099890. OpenUrl Abstract / FREE Full Text 6. ↵ Stergachis AB , Debo BM , Haugen E , Churchman LS , Stamatoyannopoulos JA . Single-molecule regulatory architectures captured by chromatin fiber sequencing . Science . 2020 Jun 26; 368 ( 6498 ): 1449 – 1454 . doi: 10.1126/science.aaz1646 . PMID: 32587015. OpenUrl Abstract / FREE Full Text 7. ↵ Tullius TW , Isaac RS , Dubocanin D , Ranchalis J , Churchman LS , Stergachis AB . RNA polymerases reshape chromatin architecture and couple transcription on individual fibers . Mol Cell . 2024 Sep 5; 84 ( 17 ): 3209 – 3222 .e5. doi: 10.1016/j.molcel.2024.08.013 . Epub 2024 Aug 26. PMID: 39191261; PMCID: PMC11500009. OpenUrl CrossRef 8. ↵ Isomura A , Kageyama R. Progress in understanding the vertebrate segmentation clock . Nat Rev Genet. 2025 Jul; 26 ( 7 ): 479 – 493 . doi: 10.1038/s41576-025-00813-6 . Epub 2025 Mar 4. PMID: 40038453. OpenUrl CrossRef 9. ↵ Kassouf M , Ford S , Blayney J , Higgs D . Understanding fundamental principles of enhancer biology at a model locus: Analysing the structure and function of an enhancer cluster at the α-globin locus . Bioessays . 2023 Oct; 45 ( 10 ): e2300047 . doi: 10.1002/bies.202300047 . Epub 2023 Jul 5. PMID: 37404089; PMCID: PMC11414744. OpenUrl CrossRef PubMed 10. ↵ Levine M . Transcriptional enhancers in animal development and evolution . Curr Biol . 2010 Sep 14; 20 (17):R754-63. doi: 10.1016/j.cub.2010.06.070 . PMID: 20833320; PMCID: PMC4280268. OpenUrl CrossRef PubMed Web of Science 11. ↵ Yáñez-Cuna JO , Kvon EZ , Stark A . Deciphering the transcriptional cis-regulatory code . Trends Genet . 2013 Jan; 29 ( 1 ): 11 – 22 . doi: 10.1016/j.tig.2012.09.007 . Epub 2012 Oct 23. PMID: 23102583. OpenUrl CrossRef PubMed Web of Science 12. ↵ Keller SH , Deng H , Lim B . Regulation of the dynamic RNA Pol II elongation rate in Drosophila embryos . Cell Rep . 2023 Oct 31; 42 ( 10 ): 113225 . doi: 10.1016/j.celrep.2023.113225 . Epub 2023 Oct 12. PMID: 37837623; PMCID: PMC10842316. OpenUrl CrossRef PubMed 13. ↵ Hunt , G. , Vaid , R. , Pirogov, S. et al. et al. Tissue-specific RNA Polymerase II promoter-proximal pause release and burst kinetics in a Drosophila embryonic patterning network . Genome Biol 25 , 2 ( 2024 ). doi: 10.1186/s13059-023-03135-0 OpenUrl CrossRef PubMed 14. ↵ Kwak H , Fuda NJ , Core LJ , Lis JT . Precise maps of RNA polymerase reveal how promoters direct initiation and pausing . Science . 2013 Feb 22; 339 ( 6122 ): 950 – 3 . doi: 10.1126/science.1229386 . PMID: 23430654; PMCID: PMC3974810. OpenUrl Abstract / FREE Full Text 15. ↵ Krebs AR , Imanci D , Hoerner L , Gaidatzis D , Burger L , Schübeler D . Genome-wide Single-Molecule Footprinting Reveals High RNA Polymerase II Turnover at Paused Promoters . Mol Cell . 2017 Aug 3; 67 ( 3 ): 411 – 422 .e4. doi: 10.1016/j.molcel.2017.06.027 . Epub 2017 Jul 20. PMID: 28735898; PMCID: PMC5548954. OpenUrl CrossRef PubMed 16. ↵ Chatsirisupachai K , Moene CJI , Kleinendorst R , Kreibich E , Molina N , Krebs A . Mouse promoters are characterised by low occupancy and high turnover of RNA polymerase II . Mol Syst Biol . 2025 May; 21 ( 5 ): 447 – 471 . doi: 10.1038/s44320-025-00094-5 . Epub 2025 Mar 31. PMID: 40164797; PMCID: PMC12048509. OpenUrl CrossRef PubMed 17. ↵ Karaiskos N , Wahle P , Alles J , Boltengagen A , Ayoub S , Kipar C , Kocks C , Rajewsky N , Zinzen RP . The Drosophila embryo at single-cell transcriptome resolution . Science . 2017 Oct 13; 358 ( 6360 ): 194 – 199 . doi: 10.1126/science.aan3235 . Epub 2017 Aug 31. PMID: 28860209. OpenUrl Abstract / FREE Full Text 18. ↵ Jonkers I , Kwak H , Lis JT . Genome-wide dynamics of Pol II elongation and its interplay with promoter proximal pausing, chromatin, and exons . Elife . 2014 Apr 29; 3 : e02407 . doi: 10.7554/eLife.02407 . PMID: 24843027; PMCID: PMC4001325. OpenUrl CrossRef PubMed 19. ↵ Fukaya T , Lim B , Levine M . Rapid Rates of Pol II Elongation in the Drosophila Embryo . Curr Biol . 2017 May 8; 27 ( 9 ): 1387 – 1391 . doi: 10.1016/j.cub.2017.03.069 . Epub 2017 Apr 27. PMID: 28457866; PMCID: PMC5665007. OpenUrl CrossRef PubMed 20. McKnight SL , Miller OL Jr . Ultrastructural patterns of RNA synthesis during early embryogenesis of Drosophila melanogaster . Cell . 1976 ; 8 ( 3 ): 305 – 319 . doi: 10.1016/0092-8674(76)90014-3 . OpenUrl CrossRef PubMed Web of Science 21. ↵ Jonkers I , Lis JT . Getting up to speed with transcription elongation by RNA polymerase II . Nat Rev Mol Cell Biol . 2015 Mar; 16 ( 3 ): 167 – 77 . doi: 10.1038/nrm3953 . Epub 2015 Feb 18. PMID: 25693130; PMCID: PMC4782187. OpenUrl CrossRef PubMed 22. ↵ Kim JM , Larson DR . Timing is everything: transcription bursting in development . Genes Dev . 2025 Sep 2; 39 ( 17-18 ): 995 – 1011 . doi: 10.1101/gad.352465.124 . PMID: 40588418; PMCID: PMC12404199. OpenUrl Abstract / FREE Full Text 23. ↵ Vollger MR , Swanson EG , Neph SJ , Ranchalis J , Munson KM , Ho CH , Cheng YHH , Sedeño-Cortés AE , Fondrie WE , Bohaczuk SC , Dippel MA , Mao Y , Parmalee NL , Mallory BJ , Harvey WT , Kwon Y , Garcia GH , Hoekzema K , Meyer JG , Cicek M , Eichler EE , Noble WS , Witten DM , Bennett JT , Ray JP , Stergachis AB. A haplotype-resolved view of human gene regulation. bioRxiv [Preprint]. 2025 Jun 2:2024.06.14.599122. doi: 10.1101/2024.06.14.599122 . PMID: 40501892; PMCID: PMC12157683. OpenUrl Abstract / FREE Full Text 24. ↵ Leptin M , Grunewald B . Cell shape changes during gastrulation in Drosophila . Development . 1990 Sep; 110 ( 1 ): 73 – 84 . doi: 10.1242/dev.110.1.73 . PMID: 2081472. OpenUrl Abstract 25. ↵ Grande MT , Sánchez-Laorden B , López-Blau C , De Frutos CA , Boutet A , Arévalo M , Rowe RG , Weiss SJ , López-Novoa JM , Nieto MA . Snail1-induced partial epithelial-to-mesenchymal transition drives renal fibrosis in mice and can be targeted to reverse established disease . Nat Med . 2015 Sep; 21 ( 9 ): 989 – 97 . doi: 10.1038/nm.3901 . Epub 2015 Aug 3. Erratum in: Nat Med. 2016 Feb;22(2):217. doi: 10.1038/nm0216-217d. PMID: 26236989. OpenUrl CrossRef PubMed 26. ↵ Ip YT , Park RE , Kosman D , Yazdanbakhsh K , Levine M. dorsal-twist interactions establish snail expression in the presumptive mesoderm of the Drosophila embryo . Genes Dev . 1992 Aug; 6 ( 8 ): 1518 – 30 . doi: 10.1101/gad.6.8.1518 . PMID: 1644293. OpenUrl Abstract / FREE Full Text 27. ↵ Hong JW , Hendrix DA , Levine MS . Shadow enhancers as a source of evolutionary novelty . Science . 2008 Sep 5; 321 ( 5894 ):1314. doi: 10.1126/science.1160631 . PMID: 18772429; PMCID: PMC4257485. OpenUrl Abstract / FREE Full Text 28. ↵ Bothma JP , Garcia HG , Ng S , Perry MW , Gregor T , Levine M . Enhancer additivity and non-additivity are determined by enhancer strength in the Drosophila embryo . Elife . 2015 Aug 12; 4 : e07956 . doi: 10.7554/eLife.07956 . PMID: 26267217; PMCID: PMC4532966. OpenUrl CrossRef PubMed 29. ↵ Birnie A , Plat A , Korkmaz C , Bothma JP . Precisely timed regulation of enhancer activity defines the binary expression pattern of Fushi tarazu in the Drosophila embryo . Curr Biol . 2023 Jul 24; 33 ( 14 ): 2839 – 2850 .e7. doi: 10.1016/j.cub.2023.04.005 . Epub 2023 Apr 27. PMID: 37116484; PMCID: PMC10373528. OpenUrl CrossRef 30. ↵ Hiromi Y , Gehring WJ . Regulation and function of the Drosophila segmentation gene fushi tarazu . Cell . 1987 Sep 11; 50 ( 6 ): 963 – 74 . doi: 10.1016/0092-8674(87)90523-x . PMID: 2887293. OpenUrl CrossRef PubMed Web of Science 31. ↵ Schier AF , Gehring WJ . Direct homeodomain-DNA interaction in the autoregulation of the fushi tarazu gene . Nature . 1992 Apr 30; 356 ( 6372 ): 804 – 7 . doi: 10.1038/356804a0 . PMID: 1574120. OpenUrl CrossRef PubMed 32. ↵ Fischer MD , Graham P , Pick L . The ftz upstream element drives late ftz stripes but is not required for regulation of Ftz target genes . Dev Biol . 2024 Jan; 505 : 141 – 147 . doi: 10.1016/j.ydbio.2023.11.004 . Epub 2023 Nov 15. PMID: 37977522; PMCID: PMC10843599. OpenUrl CrossRef PubMed 33. ↵ Dearolf CR , Topol J , Parker CS . Transcriptional control of Drosophila fushi tarazu zebra stripe expression . Genes Dev . 1989 Mar; 3 ( 3 ): 384 – 98 . doi: 10.1101/gad.3.3.384 . PMID: 2498165. OpenUrl Abstract / FREE Full Text 34. ↵ Petesch SJ , Lis JT . Rapid, transcription-independent loss of nucleosomes over a large chromatin domain at Hsp70 loci . Cell . 2008 Jul 11; 134 ( 1 ): 74 – 84 . doi: 10.1016/j.cell.2008.05.029 . PMID: 18614012; PMCID: PMC2527511. OpenUrl CrossRef PubMed Web of Science 35. ↵ Chen PT , Levo M , Zoller B , Gregor T . A conserved coupling of transcriptional ON and OFF periods underlies bursting dynamics . Nat Struct Mol Biol . 2025 Oct; 32 ( 10 ): 1959 – 1971 . doi: 10.1038/s41594-025-01615-4 . Epub 2025 Jul 15. PMID: 40664961. OpenUrl CrossRef 36. ↵ Perry MW , Boettiger AN , Bothma JP , Levine M . Shadow enhancers foster robustness of Drosophila gastrulation . Curr Biol . 2010 Sep 14; 20 ( 17 ): 1562 – 7 . doi: 10.1016/j.cub.2010.07.043 . PMID: 20797865; PMCID: PMC4257487. OpenUrl CrossRef PubMed Web of Science 37. ↵ Batut PJ , Bing XY , Sisco Z , Raimundo J , Levo M , Levine MS . Genome organization controls transcriptional dynamics during development . Science . 2022 Feb 4; 375 ( 6580 ): 566 – 570 . doi: 10.1126/science.abi7178 . Epub 2022 Feb 3. PMID: 35113722; PMCID: PMC10368186. OpenUrl CrossRef 38. ↵ Balasubramanian D , Borges Pinto P , Grasso A , Vincent S , Tarayre H , Lajoignie D , Ghavi-Helm Y . Enhancer-promoter interactions can form independently of genomic distance and be functional across TAD boundaries . Nucleic Acids Res . 2024 Feb 28; 52 ( 4 ): 1702 – 1719 . doi: 10.1093/nar/gkad1183 . PMID: 38084924; PMCID: PMC10899756. OpenUrl CrossRef 39. ↵ Mateo , Leslie J. , Nasa Sinnott-Armstrong , and Alistair N. Boettiger . “ Tracing DNA paths and RNA profiles in cultured cells and tissues with ORCA .” Nature protocols 16 . 3 ( 2021 ): 1647 – 1713 . OpenUrl PubMed 40. ↵ Schor IE , Bussotti G , Maleš M , Forneris M , Viales RR , Enright AJ , Furlong EEM . Non-coding RNA Expression, Function, and Variation during Drosophila Embryogenesis . Curr Biol . 2018 Nov 19; 28 ( 22 ): 3547 – 3561 .e9. doi: 10.1016/j.cub.2018.09.026 . Epub 2018 Nov 1. PMID: 30393032; PMCID: PMC6264527. OpenUrl CrossRef PubMed 41. ↵ Jha A , Bohaczuk SC , Mao Y , Ranchalis J , Mallory BJ , Min AT , Hamm MO , Swanson E , Dubocanin D , Finkbeiner C , Li T , Whittington D , Noble WS , Stergachis AB , Vollger MR . DNA-m6A calling and integrated long-read epigenetic and genetic analysis with fibertools . Genome Res . 2024 Nov 20; 34 ( 11 ): 1976 – 1986 . doi: 10.1101/gr.279095.124 . PMID: 38849157; PMCID: PMC11610455. OpenUrl Abstract / FREE Full Text 42. ↵ Rauluseviciute I , Riudavets-Puig R , Blanc-Mathieu R , Castro-Mondragon JA , Ferenc K , Kumar V , Lemma RB , Lucas J , Chèneby J , Baranasic D , Khan A , Fornes O , Gundersen S , Johansen M , Hovig E , Lenhard B , Sandelin A , Wasserman WW , Parcy F , Mathelier A JASPAR 2024: 20th anniversary of the open-access database of transcription factor binding profiles Nucleic Acids Res. 2024 Jan 5; 52 ( D1 ): D174 – D182 .; doi: 10.1093/nar/gkad1059 OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted October 29, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Comprehensive Single Molecule View of Transcriptional Dynamics in Development Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Comprehensive Single Molecule View of Transcriptional Dynamics in Development Thomas W Tullius , Tohn Borjigin , Michael Levine bioRxiv 2025.10.28.685222; doi: https://doi.org/10.1101/2025.10.28.685222 Share This Article: Copy Citation Tools Comprehensive Single Molecule View of Transcriptional Dynamics in Development Thomas W Tullius , Tohn Borjigin , Michael Levine bioRxiv 2025.10.28.685222; doi: https://doi.org/10.1101/2025.10.28.685222 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Developmental Biology Subject Areas All Articles Animal Behavior and Cognition (7633) Biochemistry (17685) Bioengineering (13890) Bioinformatics (41930) Biophysics (21449) Cancer Biology (18587) Cell Biology (25496) Clinical Trials (138) Developmental Biology (13375) Ecology (19897) Epidemiology (2067) Evolutionary Biology (24309) Genetics (15608) Genomics (22499) Immunology (17736) Microbiology (40390) Molecular Biology (17177) Neuroscience (88586) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4823) Physiology (7641) Plant Biology (15150) Scientific Communication and Education (2045) Synthetic Biology (4293) Systems Biology (9823) Zoology (2271)