Stress‐sensitive dynamics of miRNAs and Elba1 in Drosophila embryogenesis

To identify sensitive periods during Drosophila development, where stress exposure can induce long‐term memory, we used the position‐effect variegation strain w^m4h (Fig 1A). This strain has an inversion on the X chromosome, positioning the white gene close to the pericentric heterochromatin. The expression of white, which is needed for eye pigmentation, is therefore controlled by the centromeric chromatin state and the adult eye color can be used as a reporter for heterochromatin at this locus (Elgin & Reuter, 2013). Variegation of w^m4h is controlled by the methyltransferases Su(var) 3–9, Su(var) 4–20, E(z), and HP1 (Phalke et al, 2009). In addition, we have previously showed that this reporter is sensitive to paternal diets (Öst et al, 2014).

To enable high‐resolution mapping of sensitive periods for environmentally induced changes in variegation, we used a heat shock intervention. In contrast to other environmental challenges such as suboptimal nutrition and exposure to toxins, heat shock allows for a sharp and distinct intervention time. We performed a 1‐h heat shock during different time points of Drosophila development, throughout embryogenesis (Fig 1B and C), as well as during the larva stages (Appendix Fig S1). After assessing the eye color in 5‐day‐old male adults, we found that the first 2 h in embryogenesis is the only sensitive period for heat shock induction of long‐term effects on heterochromatin (Fig 1C; Appendix Fig S1). This finding was confirmed by repeating the experiment with an even shorter heat shock exposure (30 min; Fig 1D). Our finding is consistent with previous work identifying the first 0–3 h after fertilization as a time in which the epigenome is sensitive to heat shock stress (Hartmann‐Goldstein, 1967; Lu et al, 1998; Seong et al, 2011; Bughio et al, 2019). Importantly, we performed heat shock in more narrow intervals than previously reported (Seong et al, 2011), and while we did notice effects in 0–1 and 1–2 h, we did not detect any significant long‐lasting effects on white expression in 2‐ to 3‐h‐old embryos, indicating that in order for long‐term effects to occur, the exposure of a stressor must happen before the MBT.

The first 2 h of Drosophila embryogenesis, entailing stages 1–3 (at 22°C), is characterized by rapid mitotic cycles dependent on maternally loaded mRNAs and proteins (Fig 2A; Bushati et al, 2008; Tadros & Lipshitz, 2009; Vastenhouw et al, 2019). In concordance with the rapid cell divisions, there is no clear higher‐order chromatin architecture in this period (Li et al, 2014; Hug et al, 2017; Ogiyama et al, 2018). Chromatin states are fully established at MBT, around stage 5, when the tempo of mitosis subsides and zygotic transcription is activated (Rudolph et al, 2007; Zenk et al, 2021). As sncRNAs are highly present already in stages 1–3 and are known to modulate higher‐order chromatin structure, we hypothesized that heat shock‐induced changes of sncRNA might precede and guide the de novo heterochromatin formation. As there are no high‐resolution timelines of changes in sncRNA during these early stages of Drosophila embryogenesis, we began by performing sncRNA sequencing on single embryos, embracing stages 1–5 (Fig 2A). We found that the relative proportions (Fig 2B) and size distributions (Fig 2C) of different sncRNA classes were highly dynamic during these stages.

As expected during the maternal‐to‐zygotic transition, stage 3 is characterized by an increased proportion of degradation products (“other,” Fig 2B), which then decrease at the end of MZT. We also noted that there is a larger proportion of rRNA fragments (rsRNAs) in stage 1 than in the other stages, while tRNA fragments (tsRNAs) are present at low levels at all stages (Fig 2B and C). By contrast, piRNAs are highly expressed throughout all stages. The most striking change, however, was the proportion of miRNA. It increased from approximately 9% in stage 1 embryos to over 50% in stage 5 (Fig 2B).

While multiple miRNAs were reduced during this developmental time window, suggestive of a maternal origin, the zygotic mir‐309 cluster and a few other miRNAs showed a pronounced upregulation (Fig 2D). Controlled by Zelda, a maternally provided pioneering transcription factor (Liang et al, 2008; Fu et al, 2014), the mir‐309 cluster plays an important role in the zygotic‐driven pathway that degrades maternal transcripts (Bushati et al, 2008). Sequential comparison of stage 2 against 1, 3 against 2, 4 against 3, and 5 against 4 revealed that even though many transcripts from the mir‐309 cluster had a sharp increase in stages 4 and 5 (Appendix Fig S2C and D), there were several members of this cluster that showed a significant increase already between stage 1 and 2 (Fig 2D; Appendix Fig S2A). To test that this early activation was not an artifact driven by a few outlier miRNA sequences, we compared unique miRNA sequences per sample and stage (Appendix Figs S2 and S3; Dataset EV1). This revealed an upregulation of several unique miRNA sequences of the mir‐309 cluster and strengthened the notion that there is an upregulation of this cluster between stages 1 and 2. Previous findings have shown that members of this cluster are expressed at low levels in 0‐ to 1‐h‐old embryos and are strongly induced 2–3 h after egg laying (Aravin et al, 2003; Ruby et al, 2007; Bushati et al, 2008; Fu et al, 2014; Ninova et al, 2014). Furthermore, low levels of Zelda have been detected in the nucleus already at nuclear cycle 2 (Nien et al, 2011). Our results align with these findings and support a scenario where some members of the mir‐309 cluster are starting to be transcribed at low levels already between embryonic stages 1 and 2.

To investigate whether a heat shock in the identified sensitive period results in changes to the sncRNA‐profile at cellularization and de novo heterochromatin formation, we next heat‐shocked 0‐ to 0.5‐h‐old embryos for 30 min and then aged them to stage 5 (Figs 3A and EV1). Precisely hand‐staged single embryos—with completed cellularization typical for the later stage 5 (Bownes, 1975)—were selected for sequencing (Fig 3A). As the gene expression is very dynamic during embryonic cellularization and gene activation, we used two strategies to ensure that minor mistakes during staging would not influence the data. First, we ensured a good sample size by selecting 24 embryos for each condition. Second, in parallel to the sncRNA‐seq, we performed rRNA‐depleted sequencing of long RNA from the same single embryo. Since there is published gene expression data of stage 5 (nuclear cycle 14) that is chronologically divided into four parts (Lott et al, 2011), we reasoned that it could be used to control our staging. By looking at the linear regression of the four parts of stage 5 gene expression, we classified genes as early (slope ≤ −1), late (slope ≥ 1), or stably expressed (slope between −1 and 1; Fig EV1B and C). We could not detect any temporal bias between our groups using this classification (Fig EV1D). Neither could we detect any temporal bias between our groups when we compared it against the expression of maternally provided (Fig EV1E; Lott et al, 2011) nor early zygotic genes (transcribed within 1–2 h of embryogenesis; Fig EV1F; De Renzis et al, 2007). Thus, we concluded that there was no temporal bias between control and heat‐shocked embryos.

Analysis of the sncRNA‐seq data showed, as earlier (Fig 2B and C), that rRNA fragments and miRNAs (Fig 3B) dominate the sncRNA profiles of stage 5 embryos. We did not detect any changes in size distribution between the conditions (Fig 3B), but we found a significant increase of miRNAs in heat‐shocked samples (Fig 3C), a result of several upregulated unique miRNA sequences (Figs 3D and EV2A; Datasets EV2 and EV3). In addition, differential expression analysis of all unique reads revealed upregulation of several piRNAs, lincRNAs, and snoRNAs (Fig EV2B, E, and F), while tsRNA and rsRNA showed diverse responses (Fig EV2C and D). As expected, we found that heat shock induced more tRNA halves originating from the 5′ terminal of mature Gly‐GCC (Fig EV2G). This tsRNA, which is also called tRNA‐derived stress‐induced small RNA, is cleaved at the anticodon loop by ribonucleases like angiogenin and is generated in response to different kinds of cellular stress leading to the formation of stress granules (Emara et al, 2010). In mammals, such tRNA 5′ halves are known to be changed in sperm in response to diet and to modulate early embryonic processes (Chen et al, 2016; Nätt & Öst, 2020), but their role in Drosophila embryogenesis is unknown. Nonetheless, the most prominent effect of early embryonic heat shock on stage 5 embryos was the significant increase of 184 unique miRNA sequences, indicative of a miRNA‐dependent stress response preceding MBT and zygotic gene activation (ZGA). Importantly, we did not detect any difference in the expression of the miR‐309 cluster, accounting for 69–75% of total miRNA at stage 5 (Fig 3E).

The heat shock‐induced upregulation of miRNA could be due to either increased zygotic transcription or decreased degradation of maternal transcripts. These two scenarios were discriminated by comparing the miRNA profiles from stages 1 to 5 embryos (Figs 2D and 3F and G, Appendix Fig S2 and S3). Since we found a clear distinction between the miRNA species found at the first stage compared with stages 4–5 (Fig 2D; Appendix Fig S3), we classified miRNAs having their peak expression during stage 1 as maternally loaded and miRNAs that had their peak in later stages as zygotic. To confirm that this classification was correct, we compared our results with the miRNA expression profiles of 0‐ to 2‐ and 2‐ to 4‐h‐old embryos made by Zhou et al (2018) and found that our classifications largely align with theirs. Using this division, we found that only a few zygotic miRNAs increased in response to heat shock, whereas there was a distinct increase in maternally provided miRNAs (Fig 3F–H).

Since the w^m4h locus is known to be controlled by classical H3K9me3‐dependent mechanisms, we hypothesized that the upregulation of miRNA would result in the downregulation of factors controlling the epigenetic state at this locus. Differential expression analysis showed that heat shock significantly altered the expression of several mRNAs (Fig 4A; Dataset EV4; 1,136 up vs. 459 down). We could, however, not detect any significant changes in the expression of H3K9me3‐related epigenetic enzymes such as Su(var)3–9, HP1, Eggless (SETDB1), or ATF‐2 (Fig 4B). When comparing all differentially expressed genes with Drosophila genes annotated with the GO‐term “heterochromatin formation” (GO:00315007), we found only two downregulated genes, pgc and Elba1 (Fig 4C). In line with this, a GO‐term enrichment analysis of up‐ or downregulated genes did not show an enrichment of epigenetic factors (Appendix Fig S4A and B).

Unsupervised correlation and clustering analysis between the upregulated miRNA and downregulated long RNA revealed a cluster of strong inverse correlation (Fig 4D, cluster 1; Dataset EV5). This cluster consisted of 11 miRNA, including mir‐190, mir‐2a and b, and bantam, and 13 mRNA. Using modENCODEs temporal expression data for all genes in each RNA cluster (Fig 4E–I) we found that our identified gene cluster displays a distinct temporal expression (Fig 4E). To get further insight into this cluster, we compared it with carefully categorized pre‐MBT genes from Chen et al (2013) and found that 9 out of 13 genes in this cluster overlapped with nonpaused pre‐midblastula transition (pre‐MBT) genes (Fig 4J). It is interesting to note that genes from the identified pre‐MBT gene cluster have their dominant expression precisely overlapping the end of the stress‐sensitive period (Figs 1C and D, and 4E). Stress‐induced downregulation of such temporally expressed genes could explain why stress at this time, but not after, would modulate variegation of the w^m4h locus. We therefore investigated this cluster with 13 genes in more detail. Functionally, we found no GO‐term enrichment for this cluster, and only one of them, Elba1, has been reported to be involved in chromatin silencing (Appendix Fig S4C). More specifically, Elba1 has been shown to be a transcriptional repressor and insulator‐binding protein that works in concert with Elba2, Elba3, and Insv to ensure that there is no leakage between transcriptional units (Aoki et al, 2012; Ueberschär et al, 2019). To verify our findings from the RNA‐seq data, we repeated the heat shock intervention during the sensitive period in embryos expressing Elba1‐GFP. This time, we collected heat‐shocked and control embryos 3–3.5 h after egg laying for fixation and staining with GFP and HP1a. Late stage 5 embryos were carefully staged using a confocal microscope with HP1a as a guide (Fig 5A). We next quantified the Elba1‐GFP expression in 10 peripheral cells per embryo and, in alignment with the RNA‐seq data, we found that a pre‐MBT heat shock significantly reduced Elba1 at embryonic stage 5 (Fig 5B).

To find potential miRNA seed sequences complementary to the Elba1 transcript, we conducted an in silico analysis using TargetScan Fly's script (Agarwal et al, 2018). We found three types of target sites (7mer‐m8, 7mer‐A1, and 6mer) and identified miR‐283‐3p as having a candidate targeting 6mer to the 3′UTR (Fig 5C). Aside from miR‐283‐3p, we identified several miRNAs with candidate targeting sites in the CDS. Similar to mammals, miRNA targeting to the 3′UTR is most effective, but targeting sites within the CDS also have some silencing effects (Agarwal et al, 2018). While multiple miRNA candidates presumably have some silencing effect on Elba1, miR‐283‐3p remains the primary candidate. This is supported by its highest negative correlation to the Elba1 transcript in our earlier correlation analysis.

To experimentally test the causal relationship between heat shock‐induced miRNA and downregulation of Elba1 (Fig 4D; Appendix Fig S5), we used qPCR to compare the Elba1 expression in stage 5 embryos of control and heat‐shocked w^m4h and Dicer‐1 (Dcr‐1) mutants. Dcr‐1 is essential and specific for the miRNA synthesis, with minimal impact on the synthesis of other sncRNA (Lee et al, 2004). As before, a short 30 min heat shock reduced the Elba1 expression in w^m4h flies (Fig 5D). In Dcr‐1‐mutant embryos, however, this downregulation was not detectable (Fig 5D). Thus, we concluded that the heat shock‐induced downregulation of Elba1 relies on a functional miRNA‐processing pathway. To directly test whether the Ago1‐miRNA RISC complex binds Elba1, we performed an Ago1 immunoprecipitation followed by qPCR using an IgG antibody as negative control (Fig 5E and F; Appendix Fig S6). In line with our previous results, we detected an enrichment of Elba1 in the Ago1 IP fraction compared with input, the Ago1 unbound fraction (UB), and the IgG IP control (Fig 5F).

Considering that the Elba complex is important to minimize transcriptional leakage (Aoki et al, 2012; Ueberschär et al, 2019), we reasoned that the heat shock‐induced upregulation of genes (Fig 4A) might be a consequence of reduced levels of Elba1. Unsupervised cluster analysis of Pearson's r scores between insulator‐binding factors and heat shock‐induced genes revealed four gene clusters of which cluster 2, the largest cluster, had a strong negative correlation to Elba1, Elba2, Insv, and CP190 (Pearson's r: Elba1 mean = −0.66, SD = 0.10; Elba2 mean = −0.76, SD = 0.11; Insv mean = −0.68, SD = 0.10 and CP190 mean = −0.63, SD = 0.11; Fig 6A, cluster 2; Dataset EV6). Functional analysis of clusters showed an over‐representation of genes involved in developmental and morphological progress in cluster 2 (Fig 6B, top right). As developmental and morphological associated genes were identified by Ueberschär et al (2019) to be controlled by the Elba complex, we compared our data with their ChIP‐data for Elba1‐3 and Insv (Appendix Fig S7A–E). In agreement with the loss of Elba1‐restricted transcription, we found an overlap between genes in cluster 2 and genes associated with Elba1‐3 ChIP‐peaks.

To test whether an early heat shock decreases the Elba1 binding to these genes, we performed CUT&RUN using GFP antibody on five sets of 20 stage 5 Elba1‐GFP embryos exposed to either 30 min heat shock during the sensitive period or kept as controls. As we had little starting material, we merged the 5 replicates within each experimental group (individual datasets are available under Data availability). Looking at the peak scores associated with these gene regions, we detected increased Elba1 binding at the TSS region of cluster 2 genes, compared with genes from the other clusters (Figs 6C and EV3A and B). Heat‐shocked embryos showed reduced peak scores at the TSS region of cluster 2 genes, despite having similar profiles (Figs 6C and EV3A and B). CUT&RUN tracks over representative genomic loci are found in Fig EV3C–E. We further aligned the consensus peaks and compared them to all upregulated genes. As before, the greatest overlap was detected between Elba1‐GFP peaks and cluster 2 (Fig 6D).

Since not only Elba1 but also Elba2 and Elba3 showed a similar clustering (Fig 6A), we next looked more closely at insulator‐binding factors with different temporal expression profiles during embryogenesis (Fig EV4). We analyzed their expression between the heat‐shocked and control embryos and found that several insulating binding factors were statistically significantly downregulated following heat shock (Fig EV4). Most interestingly, we found that the expression of all members of the Elba complex and Insv was reduced, although not enough (with the exception of Elba1) to reach the Log₂ fold change threshold of < −1 or the FDR used when analyzing the whole dataset.

To test whether the reduced expression of the Elba complex during MBT would result in long‐term epigenetic effects in the adult fly, we next crossed female virgins from the w^m4h PEV‐strain with homozygous Elba1‐3 mutant males (Fig 6E). Intriguingly, eye pigmentation of 5‐day‐old males showed loss of white silencing in all Elba heterozygous mutants, thus mimicking the effect of an early embryonic heat shock (Fig 6F and G). In all, our study points to a central role for miRNA‐Elba dependent fine‐tuning of the emerging chromatin landscape, that—if disrupted—may have life‐long consequences on the phenotype.

The ln(1)w^m4h Drosophila strain (Muller, 1930) was kindly provided from Gunter Reuter's lab and has been maintained in a climate‐controlled 22°C incubator and kept on standardized food. Flies used for experiments were inbred for > 10 generations and flies with complete loss of PEV were not used for crossing to enable capture of differences of variegation. Dicer‐1^Q1147X mutants (Lee et al, 2004; #RRID:BDSC_32066) containing a nonsense codon at the PAZ domain were kept at 22°C for > 6 generations on standard food before egg collection. Elba1 sk6, Elba2 sk2, and Elba3 sk5 mutant flies homozygously carrying the respective frame‐shift mutation were kindly provided from Qi Dai's lab (fly strains are described in Ueberschär et al, 2019) and kept in room temperature (at approximately 22°C) on standardized food. Elba1‐3^−/− − ln(1)w^m4h crossings were kept in a climate‐controlled 22°C incubator on standardized food. The Elba1‐GFP fly strain (RRID:BDSC_83657) and w¹¹¹⁸ were kept in a 26°C incubator on standardized brown food.

For screening of sensitive periods: eggs were collected on juice agar plates in tight intervals (30 min—1 h) and exposed to one heat shock session in a 37°C incubator, or were kept as controls. The selection was random. For PEV expression in Elba1 mutants: virgin ln(1)w^m4h were crossed with either Elba1‐3 mutant males or ln(1)w^m4h males and left to mate and lay eggs. Five or six different vials per crossing were set up and all were flipped 3 times. All experiments: flies were left to develop in a climate‐controlled 22°C incubator. Males were decapitated 5 days after eclosure and their heads, collected in groups of 3–10, were first frozen in liquid nitrogen and then homogenized with a 5 mm ∅ metal bead (Qiagen) for 2 min at 40 Hz using TissueLyser LT (Qiagen). 500 μl PBS‐tween (0.01%) was added and samples were shaken, kept at room temperature for 1 h, and centrifuged. Absorbance at A480 was measured on supernatant in technical doublets using VersaMax (Molecular Devices) microplate reader. ln(1)w^m4h heat shock experiments: At least two biological replicates were collected per heat shock time and experiment, and heat shock experiments were performed 7 times. Elba1‐3 mutants—ln(1)w^m4h experiments: 6–10 heads were analyzed per sample. n = 58 w^m4h × w^m4h‐, n = 51 Elba1^SK6 × w^m4h‐, n = 19 Elba2^SK2 × w^m4h‐, and n = 13 Elba3^SK5 × w^m4h crossings were collected and analyzed.

Eggs were collected on juice agar plates for 30 min and were immediately dechorionated. The staging was performed under SMZ 745 (Nikon) microscope using the criteria for Bownes' stages 1–5 (Bownes, 1975). Single embryos were collected in 2 μl RNase‐free water with Recombinant RNase inhibitor (TAKARA) and ruptured with an RNase‐free needle. One 5 mm ∅ metal bead (Qiagen) and 500 μl Qiazol (Qiagen) were added per sample and the samples were homogenized for 2 min at 40 Hz using TissueLyser LT (Qiagen). n = 5 of stage 1–3 and n = 4 of stage 4 and 5.

Eggs were collected on juice agar plates in 30 min intervals and immediately exposed to one session of heat shock at 37°C for 30 min or kept as controls. This selection was random. Embryos were thereafter kept in a climate‐controlled 22°C incubator for approximately 2 h, dechorionated in 3.5% bleach, and staged under SMZ 745 (Nikon) bright‐field microscope using the criteria for Bownes' stage 5 (Bownes, 1975), including formed cells at egg surface and round pole cells at the posterior axis. Single embryos were collected in 2 μl RNase‐free water with an RNase inhibitor and ruptured with an RNase‐free needle. One 5 mm ∅ metal bead (Qiagen) and 500 μl Qiazol (Qiagen) were added per sample and the samples were homogenized for 2 min at 40 Hz using TissueLyser LT (Qiagen). n = 24 single embryos per condition.

RNAs were extracted using miRNeasy Micro Kit (Qiagen) according to manufacture protocol. Quality was confirmed using Agilent RNA 6000 Nano kit (Agilent) on the 2100 Bioanalyzer Instrument (Agilent) prior to storage at −70°C. NEBNext Small RNA Library Prep Set for Illumina (New England Biolabs) was used for library preparation according to the manufacturer's protocol with some changes. We downscaled all samples to half volume and added 2S rRNA block oligo (5′‐TAC AAC CCT CAA CCA TAT GTA GTC CAA GCA‐SpcC3 3′; 10 μM; Wickersheim & Blumenstiel, 2013) to a final concentration of 2.5 μM together with SR‐RT primer (from the kit). Primers and adaptors from the kit were diluted 1:4 until PCR amplification, according to starting RNA concentration. PCR amplification was run for 15 cycles and NEBNext Index1‐24 primers for Illumina were used (New England Biolabs).

Libraries were cleaned using Agencourt AMPure XP (Beckman Coulter) and run on precasted 6% polyacrylamide Novex TBE gel (Invitrogen). Bands of sizes 140–170 bp were selected. Gel extraction was made by centrifugation at 15,000 g using gel breaker tubes (IST Engineering Inc) in DNA Gel Elution Buffer provided in the NEBNext kit. Samples were incubated at 37°C for 1 h on a shaker, frozen at −70°C for 15 min, and incubated at 37°C on a shaker for 1 h once more. Gel debris was removed by Spin‐X 0.45 μm tube. Libraries were precipitated overnight at −70°C in 1 μl GlycoBlue (Invitrogen), 0.1 × volume of 3 M acetate (pH 5.5), and 3 × volume of 100% ethanol. Library sizes were measured on 2100 Bioanalyzer instrument (Agilent) using the Agilent High Sensitivity DNA kit (Agilent) and concentration was determined using QuantiFluor ONE ds DNAsystem on Quantus fluorometer (Promega). Equal concentrations of libraries were pooled and sequenced on NextSeq 500 sequencer using NextSeq 500/550 High Output Kit v2 with 75 cycles (Illumina). Unique sample IDs are summarized in Dataset EV8.

We used Cutadapt version 1.18 (Martin, 2011) to trim the adaptor sequence (AGATCGGAAGAGCACACGTCTGAACTCCAGTCACAT) from sncRNA reads and FastQC v.0.11.5 (Andrews, 2015) for quality filtering. Reads between 14‐ and 80 nucleotides containing the adaptor and with more than 80% of the bases having a phred quality score (Q‐score) > 20 were retained. Mean sequence depth was 18.06 M reads (min = 12.01 M, max = 38.36 M) for the stage 1–5 data and 15.98 M reads (min = 12.15 M, max = 20.21 M) for the heat shock experiment data.

Trimmed reads were further mapped using SPORTS pipeline version 1.0.5 (Shi et al, 2018) with standard settings except following modifications; we replaced Rfam with repeatmasker (Dataset EV9), which was placed at the bottom of the hierarchy. Within this pipeline, Bowtie version 1.1.2 (Langmead et al, 2009) was used with the following input; −M 1 ‐‐strata ‐‐best ‐v 1, returning one single read allowing one mismatch. Alignment was performed to the dm6 reference genome and then to sncRNA‐specific references using the following hierarchy; miRNA, tRNA, rRNA, piRNA, other ncRNA, and repeats. For details and annotation sources, see Dataset EV9. The number of reads, read‐length, and annotation hits were retained per unique sequence.

A list containing experimental metadata, annotation information, and count table was retained and filtered; minimum of 20 reads per sequence in 17% of samples for stage 1–5 timeline experiment and minimum of 20 reads per sequence in 50% of samples for heat shock experiment. An additional filter removing sequences with less than 0.01 rpm per sequence (in 17% of samples for the timeline experiment or 100% of samples for heat shock experiment) was applied and sequences were assigned to a sncRNA class according to regular expressions retained from the annotation hits. Differential expression analysis was performed in DEseq2 (version 1.24.00) with the design ~ Intervention + flow cell run. We used Euclidean clustering within the pheatmap package (version 1.0.12) and a cutoff at Log2 ± 1 based on stage 5 vs. stage 1 differential expression to determine the maternal or zygotic origin of miRNAs.

DNA was digested from aliquots from the same RNA extracted as described above (RNA extraction and smallRNA library preparation) with RNase‐Free DNase Set (Qiagen) according to kit protocol and concentrated using Oligo Clean & Concentrator (Zymo Research) according to kit protocol but adjusted for sample volumes. RNA quality was determined on the 2100 Bioanalyzer Instrument (Agilent) using RNA 6000 Nano kit (Agilent).

cDNA was synthesized using Ovation RNA‐Seq Systems 1–16 for model organisms (NuGEN) according to the kit protocol. Samples were sonicated 6 times in 15 s on‐ 15 s off‐intervals using the Bioruptor Pico sonication device (diagenode). Library construction was done using the Ovation RNA‐Seq Systems 1–16 for model organisms (NuGEN) according to protocol, and cycles for library amplification were determined using 7900HT Fast Real‐Time PCR System (Applied Biosystems™). The amplification buffer and enzyme mixes provided in the library kit were used for the master mix together with EvaGreen for qPCR (Biotium). Libraries were amplified according to the mean cycle for exponential PCR amplification per experiment (16 cycles) and purified according to protocol. Library sizes were measured on the 2100 Bioanalyzer Instrument (Agilent) using the High Sensitivity DNA chip (Agilent) and concentrations were determined using QuantiFluor ONE ds DNAsystem on Quantus fluorometer (Promega). Equal concentrations of libraries were pooled and sequenced on the NextSeq 500 sequencer using NextSeq 500/550 High Output Kit v2 with 75 cycles (Illumina). Unique sample IDs are summarized in Dataset EV8.

We used Cutadapt version 1.18 (Martin, 2011) to trim the adaptor sequence (AGATCGGAAGAGCACACGTC) from long RNA reads and FastQC v.0.11.5 (Andrews, 2015) for quality filtering. Reads over 14 nucleotides and with more than 80% of the bases having a phred quality score (Q‐score) > 20 were retained. Depth per library was 21.54 M reads (min = 20.26 M, max = 31.66 M reads). STAR genome index files were generated using Drosophila_melanogaster.BDGP6.28.dna.toplevel fasta and Drosophila_melanogaster.BDGP6.28.101.gtf (Ensembl). These genome index files were then used on trimmed reads using STAR (v.2.5.0a; Dobin et al, 2013) with standard settings and indexed using SAMtools (v.1.3.1; Li et al, 2009a). Standard featureCounts (v.1.5.0‐p1; Liao et al, 2014) settings were used for assigning reads to genomic features, with minimal overlap of 15 bases. A list containing experimental metadata, annotation information, and count table was retained and filtered (minimum 10 counts per read in 50% of samples).

To compare maternally loaded and early zygotic genes between control and heat‐shocked embryos, we extracted the genes in our data that matched the classifications made by Lott et al (dataset S1 in Lott et al, 2011; maternal) or by De Renzis et al (table S8 in De Renzis et al, 2007; early zygotic). We further used the expression profile of sub‐stages of nuclear cycle 14 from dataset S1 in (Lott et al, 2011) to differentiate nuclear cycle 14 expressed genes. The web‐based gene set analysis toolkit (WebGestalt; Wang et al, 2017) was used for functional analysis using the over‐representation analysis of biological processes with the BH FDR method. To generate matrixes of correlation, we used the rcorr function with the Pearson's option within the Hmisc package (version 4.2‐0) and Euclidean clustering within the pheatmap package (version 1.0.12). For analysis of overlaps with pre‐MBT classifications, we used classifications made by Chen et al (table 1 in Chen et al, 2013), to compare against genes from indicated gene clusters.

Eggs from w^m4h, Dcr‐1‐ and Elba1 mutants were collected on juice agar plates for 1 h and immediately heat‐shocked for 30 min (as earlier described) or kept as control. All samples were dechorionated in 3.5% bleach and staged under an SMZ 745 (Nikon) microscope as described above. 2–5 stage 5 embryos were collected per sample in 4 μl RNase‐free water with RNase inhibitor and ruptured with an RNase‐free needle. One 5 mm ∅ metal bead (Qiagen) and 500 μl Qiazol (Qiagen) were added per sample and shaken for 2 min at 40 Hz using TissueLyser LT (Qiagen). RNAs were extracted using miRNeasy Micro Kit (Qiagen) according to manufacture protocol and good quality was confirmed using Agilent RNA 6000 Nano kit (Agilent) on the 2100 Bioanalyzer Instrument (Agilent). iScript cDNA Synthesis Kit (BIO‐RAD) was used for cDNA synthesis and a master mix of iTaq Universal SYBR Green Supermix (BIO‐RAD), RNase‐free H₂O, and primers (Merch) was prepared according to manufacturer's protocols. Samples and master mix were loaded in triplicates onto 96 well plates and read on a 7500 Fast Real‐Time PCR system (Thermo Fisher Scientific). Primers used: Elba1 forward: TGTCCTTAGCAGCTTCTCAG, reverse: CGCATTCAAGATGCAAATGAG. Dicer‐1 forward: AGGAGACAAAGCGGGCAAAG, reverse: TATGCGGTACAGGATGCAGG. RpL32 (for normalization) forward: CCGCTTCAAGGGACAGTATC, reverse: ACGTTGTGCACCAGGAACTT. Relative expression toward w^m4h control embryos was analyzed using the ΔΔCT method.

W¹¹¹⁸ Drosophila embryos were collected in 1 h intervals, aged for 1–1.5 h, dechorionated in 3.5% bleach, and washed with RNase‐free H₂0. Embryos were put in RNase‐free PBS and transferred in batches of 100–500 to Eppendorf tubes where the PBS was discarded. PBS with cOmplete Mini EDTA‐free protease inhibitor cocktail (Sigma) and RNase inhibitor was added to each batch and the samples were snap frozen. 1,000 embryos were pooled into one tube (n = 8), the PBS was removed and lysis buffer (50 nM Tris (pH 7.5), 100 mM KCl, 12 mM MgCl₂ 1% Nonidet P‐40, 1 mM DTT, 100 μg/ml Cyclohexamide, 1× cOmplete Mini EDTA‐free protease inhibitor cocktail and 1 μl/ml RNase inhibitor) was added. The samples were homogenized using the “tight” Dounce Tissue Grinder and centrifuged at 4°C for 10 min at 10,000 g. The supernatants were collected and 10% of each sample was collected as input. Of each, one half was used for RNA extraction and 400 μl Qiazol was added and the samples snap frozen on dry ice. For the other half, 1× NuPAGE LDS Sample Buffer (Invitrogen) and 0.5 μl 2‐mercaptoethanol were added. Protein samples were boiled at 80°C for 10 min and put on dry ice.

For immunoprecipitation, the embryo lysates were initially precleared using Dynabeads™ Protein G (Invitrogen) for 30 min at 4°C. The lysates were divided in two and incubated with rabbit anti‐Ago1 polyclonal (Abcam, ab5070) or rabbit IgG polyclonal (Abcam, 171870) on rotation at 4°C overnight followed by incubation with Protein G‐dyna beads for 2–4 h at 4°C on rotation. Unbound sample (UB) was collected and prepared for RNA or protein extraction (as described above). The beads were washed in high salt buffers (50 mM Tris (pH 7.5), 300 mM KCl, 12 mM MgCl_2, 1% Nonidet P‐40, 1:250 DTT, 0.5 μl/ml RNase inhibitor, 100 μg/ml Cyclohexamide and 1× cOmplete Mini EDTA‐free protease inhibitor cocktail) and extra high salt buffer (high salt buffer +300 nM NaCl). We diluted an aliquot of beads in lysis buffer and used those for western blot (as described below). The rest of the beads were resuspended in 400 μl Qiazol and processed for RNA extraction, quality assessment, and qPCR, as described earlier (see Quantitative Real‐Time PCR). Equal loading of tot RNA was assured by QuantiFluor ONE ds DNAsystem on Quantus fluorometer (Promega). Primers used: Elba1 forward: TGTCCTTAGCAGCTTCTCAG, reverse: CGCATTCAAGATGCAAATGAG. RpL32 (reference gene) forward: CCGCTTCAAGGGACAGTATC, reverse: ACGTTGTGCACCAGGAACTT. Enrichment was calculated by 2^{−(mean Ct (sample) – mean Ct(sample input))}.

For IP assay control, SDS–PAGE, protein transfer onto PVDF membrane, and western blotting were all performed using standard procedures on the prepared protein extracts. Primary rabbit anti‐Ago1 polyclonal (Abcam, ab5070) antibody in the dilution of 1:1,000 was used, and the proteins detected using the secondary antibody from Li‐COR diluted 1:15,000 (IRdye donkey anti‐rabbit 800CW). Nonspecific proteins were blocked using 4% BSA in PBS with 0.1% tween and all antibody dilutions were made in the same block solution. PBS with 0.1% tween was used as a wash buffer. Membranes were imaged using Odyssey^® DLx Imaging System (Li‐COR).

The Elba1 transcript was downloaded from FlyBase with transcript ID FBtr0077423. The seed sequences of all upregulated miRNA families with a strong negative correlation to Elba1 (Pearson's r < −0.5) were obtained from TargetScan Fly release 7.2 (Agarwal et al, 2018). The miRNA seed target analysis using the full Elba1 transcript was done using TargetScan Fly's script, targetscan_70.pl.

Eggs from Elba1‐GFP flies were collected on juice agar plates for 1 h and heat‐shocked at 37°C for 30 min (as described above) or collected for 2 h and kept at 26°C as controls. Eggs were detached from plates, rinsed, and dechorionated in 3.5% bleach for 3–3.5 h after respective cage setup. Eggs were rinsed in water and fixated for 30 min in a 1:1 solution of 4% PFA and n‐heptane. The PFA layer was removed and an equal volume of 99.9% methanol was added. For removal of the vitelline membrane, the vials were hand‐shaken for ~60 s and moved to new tubes. Eggs were washed several times with 99.9% methanol and stored at −20°C. The samples were stepwise rehydrated in 80/20, 60/40, 40/60, and 20/80 mixture of 99.9% methanol and 0.2% PBT (PBS + Triton X‐100) for 5 min each and then blocked in PBTN (0.2% PBT with 4% horse serum) for 1–3 h RT. The samples were incubated overnight at 4°C with primary antibodies (Rabbit anti‐GFP (Torrey Pines Biolabs # TP401) and mouse anti‐C1A9‐s (1ea; Developmental Studies Hybridoma Bank)) diluted 1:500 with PBTN. Samples were then washed x 4 in PBT and incubated with secondary antibodies (Alexa fluor 488 donkey anti‐rabbit; Life Technologies # A21206) and Rhodamine Red‐X‐conjugated donkey anti‐mouse (Jackson immunofluorescence laboratories # 715‐295‐150) diluted 1:500 with PBTN for 2 h RT. The samples were washed × 4 in PBT and finally mounted on slides in Vectashield Vibrance Antifade Mounting Medium (Vector Laboratories, Inc. Ref: H‐1700). We conducted a pilot trial mounting embryos in Vectashield containing DAPI to facilitate staging. However, the emission spectra from DAPI interfered with the detection of the Alexa fluor 488, and we therefore decided to exclude DAPI from further experiments.

For reliable quantification, we captured one middle stack intersection of control and heat‐shocked stage 5 embryos at 20x on an LSM700 upright confocal microscope (ZEISS microscopy) using the same gain and acquisition settings between all images. Fijis ImageJ2 (version 1.53f55, java version: 1.8.0_172) rolling ball radius algorithm was used for background subtraction and kept constant between samples. The Elba1‐GFP expression was quantified from 10 adjacent nuclei per embryo, using the freehand line option where the line width was set according to nucleus size. The integrated density was used to calculate the relative expression between the two conditions. For visualization purposes only, the brightness and the contrast for each channel were similarly modified for both samples using Fijis ImageJ2 brightness/contrast function.

Elba1‐GFP and w¹¹¹⁸ eggs were collected on juice agar plates in 45 min intervals and immediately exposed to one session of heat shock at 37°C for 30 min or kept as controls (as described above). Embryos were thereafter kept in a climate‐controlled 22°C incubator for approximately 2 h, dechorionated in 3.5% bleach, and staged under SMZ 745 (Nikon) bright‐field microscope using the criteria for Bownes' stage 5 (Bownes, 1975). 20 embryos were collected in 140 μl nuclear extraction buffer (described in Zambanini et al, 2022) and ruptured with an RNase‐free needle. n = 5 samples of 20 embryos per condition (Elba1‐GFP), n = 2 w¹¹¹⁸ (no GFP‐control) samples. The samples were centrifuged for 10 min 700 g at 4°C and the supernatant was discarded. Pellet was resuspended in 100 μl nuclear extraction buffer. Bead‐, primary antibody, and pAG‐MNase binding, digestion, fragment release, beads clean up, library preparation, and gel extraction was made exactly similar to the CUT&RUN low volume‐Urea protocol described in (Zambanini et al, 2022) with the addition of adding 0.1 ng/ml CUTANA™ Ecoli spike‐in to the stop buffer mix, and using the rabbit‐anti‐GFP (Abcam, ab290) 1:200. Concentrations were determined using QuantiFluor ONE ds DNAsystem on Quantus fluorometer (Promega) and equal library concentrations were pooled and sequenced (paired‐end) on the NextSeq 500 sequencer using NextSeq 500/550 High Output Kit v2.5 with 75 cycles (Illumina).

Quality control was made using FastQC (v.0.11.5; Andrews, 2015) and adaptor trimming using BBDuk from the BBTools suite (v. 39.01; Bushnell et al, 2017). Post trim QC, aligning target and spike‐in, spike‐in and w¹¹¹⁸ GFP‐control normalization, peak calling, and consensus peak reporting were made using nf‐core's CUT&RUN v3 pipeline (Ewels et al, 2020), using standard settings, setting the iGenome reference to dm6 (spike‐in: K12‐MG1655), and using the de‐duplication of target. We used macs2 as the primary peak caller and included peaks found in ≥ 2 samples as the threshold for consensus peaks.

As the starting material was low, we merged the normalized bigWig files per experimental condition. These files were used to compute Matrix and plot heatmaps and profiles using deepTools (v. 3.5.1; Ramírez et al, 2016), as well as to demonstrate representative genomic areas using IGV (v. 2.14.1). Unique sample files can be found under BioProject: PRJNA729249 (see Data availability). For overlaps between heat shock‐induced upregulated genes, CUT&RUN peaks and Insv and Elba factor‐binding sites, genes notated with corresponding binding sites in Ueberschär et al (supplementary data 4 in Ueberschär et al, 2019) was extracted. The consensus peaks (Dataset EV7) from the CUT&RUN experiment were aligned to the closest gene using ChIPseeker's (v.1.28.3; Yu et al, 2015) annotatePeak function using default values with TxDb.Dmelanogaster.UCSC.dm6.ensGene as reference. UpSetR (v. 1.4.0; Conway et al, 2017) was used to illustrate the intersections to the indicated clusters.

All statistical analysis was done in R 3.6.0, R 4.1.0, or GraphPad Prism v.8.4.3 and v.9.1.2. For eye pigment statistical analysis, ordinary one‐way ANOVA with the Dunnett's multiple comparison or two‐tailed t‐test was used as indicated. Four outliers (1 from dataset used in Fig 1C, at 12 h, and 3 from dataset in Fig 6F, one from w^m4h, Elba2^−/+, Elba3^−/+ each) was removed using the ROUT method (Q = 0.1%). For all qPCR measurements and fluorescence quantification, we used the unpaired one‐ or two‐tailed student t‐tests or Mann–Whitney test depending on the normal distribution as measured with the D'Agostino–Pearson normality test or Shapiro–Wilk normality test. As indicated, we used either rpm (Dataset EV1–EV4) or variance stabilizing transformation (vst) from DEseq2 (version 1.24.00) for normalization of sncRNA and long RNA‐sequencing results. For statistical analysis of sncRNA and long RNA‐seq data, the DEseq2's build‐in Wald test after negative binominal fitting was used. The unpaired one‐ or two‐tailed t‐test or Mann–Whitney test was used to test expression rates of specific targets where indicated.