An hourglass pattern of inter-embryo gene expression variability and of histone regulation in fly embryogenesis

Gene expression is in part a stochastic process 1–3, which is at odds with the need for precise expression levels in biological processes such as embryonic development. To study the tension between precision and stochasticity of gene expression during development, we quantified expression variability of single embryo transcriptomes throughout fly Drosophila melanogaster embryogenesis, from 150 samples over eight developmental stages. We found that expression variability follows an hourglass pattern, with lower variability at extended germband, the phylotypic stage. We explain this pattern by stronger histone modification mediated transcriptional noise control at this stage. In addition, we propose that histone modifications contribute to mutational robustness in regulatory elements, and thus to conserved expression levels. These results provide insight into the role of robustness in the phenotypic and genetic patterns of evolutionary conservation in animal development.

p-values (Benjamini-Hochberg method) are shown in Table S2; they are all < 10 −5 except 95 E2 vs. E4, for which p-values = 0.24. 96 97 The variation in expression variability could either be due to changes in the set of active 98 genes, with genes differing in their intrinsic variability levels, or to genome-wide changes 99 in the regulation of variability. To test this, we first reproduced our results restricted to the 100 subset of genes which are expressed at all stages. Under the first explanation, we would 101 expect to lose the hourglass variability pattern, but the pattern is maintained ( Figure 2B). 102 We performed additional tests: restricting to genes with constant expression level over 103 development ( Figure  promoters seem to be more robust to mutations, which might also translate into robustness 106 to noise. Despite a loss of power with fewer genes, there remains an hourglass pattern of 107 expression variability in all cases. Interestingly, the precise promoter genes have higher 108 variability than the dispersed promoter genes except at E3, thus a strongest hourglass pattern. 109 Overall, these results suggest that the lower variability at E3 is due to genome-wide 110 regulation mechanisms more than to changes in the gene set. For each gene, we calculated the mean modification signal (background-subtracted tag 118 density) separately for proximal promoters and for gene bodies 23 . Higher modification signal 119 genes tend to have lower variability for all histone modifications ( Figure 3A). This supports a 120 role in minimizing transcriptional noise, and is consistent with previous studies in yeast and 121 mammals 22,23 .

H3K4Me3
Zscore (promoter) we normalized gene and promoter signal by that on intergenic regions (Methods), which are 152 not expected to change histone modification signal between stages. All histone marks present 153 an hourglass-like pattern, with the highest signal at 8-12h (except for H3K4Me1 on gene body,  154 where it is a local but not global maximum), corresponding to E3 and E4, i.e. the lowest 155 expression variability, for both promoters and gene body (Figures 3B-C, S11). Moreover, for 156 all histone marks on gene body, as well as H3K4Me1 on promoters, there is another local 157 maximum at 16-20h, corresponding to E6. Generally, these results support changes in histone 158 modification signal over development, with a correspondence between stronger histone 159 modification signal and lower expression variability. 160 161 Several studies have suggested that mechanisms which confer robustness to stochastic 162 variation can also buffer the effects of genetic variation 7,27,28 . If histone modifications can 163 buffer the effect of genetic variation on gene expression, we should observe that genes with 164 higher histone modification signal are less sensitive to mutations in their regulatory regions, 165 and are thus less conserved. Indeed, genes with higher histone modification signal tend to 166 have less conserved core promoter sequences 29 (49 bp upstream TSS and 10 bp downstream 167 from the TSS) between species (phastCons score; Figure 4A). They are also less conserved 168 within a population (promoter nucleotide diversity π; Figure S12). The phastCons pattern 169 remains using 200 bp or 400 bp regions, but disappears using 1 kb regions ( Figure S13   results were based solely on promoter analysis, and it remains to be seen how much these 217 observations extend to other regulatory elements. Second, we only considered nucleotide 218 substitution changes, but not insertion, deletion and turnover. For these large effect 219 mutations, purifying selection might be more efficient than mutational robustness to keep 220 conserved expression level. Thus it is possible that both purifying selection and mutational 221 robustness together shape the lower expression divergence in the fly phylotypic stage. We performed pairwise Wilcoxon test between any two stages to test the significance. The 287 multiple test corrected p-values (Benjamini-Hochberg method) are shown in Table S9. The legend is the same as for Figure 2. We performed pairwise Wilcoxon test between any 294 two stages to test the significance. The multiple test corrected p-values (Benjamini-Hochberg 295 method) are shown in Tables S10 and S11. 296 For each stage, the first and the second box represents dispersed promoter genes and precise 300 promoter genes respectively. The legend is the same as for Figure 2. We performed pairwise 301 Wilcoxon test between any two stages to test the significance separately for dispersed promoter 302 genes and for precise promoter genes. The multiple test corrected p-values (Benjamini-

Figure S11: Histone modification signal across development 306
The legend is the same as for Figure 3B and 3C. The median signal value in each 307 development stage is indicated above each box. We performed pairwise Wilcoxon test 308 between any two stages to test the significance. The multiple test corrected p-values 309 (Benjamini-Hochberg method) for H3K4Me1, H3K27Ac and H3K9Ac are shown in Tables 310 S14-S19. 311 312 Figure S12: Spearman's correlation coefficient between histone modification signal 313 and promoter nucleotide diversity (π). 314 The legend is the same as for Figure 4A. The figure legend is the same as in Figure 4A. The legend is the same as for Figure 4C. We calculated the ratio of mean expression between genes from the X chromosome and from 334 the autosomes for each sample. Red represents high ratio, blue represents low ratio. For 335 Drosophila, dosage compensation is achieved by increasing expression of X chromosome 336 genes in males. Since the dosage compensation is still incomplete during development, 337 females should have a higher ratio of mean expression between genes from the X chromosome and from the autosomes. Here, we found both high ratio samples and low ratio samples are 339 well mixed in both the cluster and large clusters. Thus, we reject the hypothesis that the two 340 different clusters are due to sex. (400 g/L) in 95% ethanol, 1 L Water. 100 to 150 flies were transferred to cages, which were 495 sealed to a grape agar plate (1:1 mixture of 6% agar and grape juice). We used 4 separate cages 496 to collect the embryos. The adults were kept overnight on this plate before being transferred to 497 a new plate supplemented with yeast paste. Synchronization of eggs on this plate lasted for 2 498 hours before being swapped with a new plate supplemented with yeast paste. We let the adults 499 lay eggs for 30 min before removing the plate and letting the eggs incubate for the desired time. 500 Eggs were harvested using the following protocol. First a 1:1 bleach (Reactolab 99412) 1x PBS 501 mix was poured on the plate and incubated for 2 min. During this incubation, we used a brush 502 to lightly scrape the surface to mobilize the embryos. We then poured the PBS-bleach mixture 503 through a sieve, washed the plate with 1x PBS, and poured the wash on the same sieve. We 504 washed the sieve several time with 1x PBS until the smell of bleach disappeared. Single 505 embryos were then manually transferred to Eppendorf containing 50 µL beads and 350 µL 506 Trizol (lifetechnologies AM9738). The tubes were homogenized in a Precellys 24 Tissue 507 Homogenizer at 6000 rpm for 30 seconds. Samples were then transferred to liquid nitrogen for 508 flash freezing and stored at -80°C. For RNA extraction, tubes were thawed on ice, 509 supplemented with 350 µL of 100% Ethanol before homogenizing again with the same 510 parameters. We then used the Direct-zol™ RNA Miniprep R2056 Kit, with the following 511 modifications: we did not perform DNAse I treatment, we added another 2 min centrifugation 512 into an empty column after the RNA Wash step, finally elution was performed by adding 8 µL 513 of RNAse-free water to the column, incubation at room temperature for 2 min and then 514 centrifugation for 2 min. RNA was transferred to a low-binding 96 well plate and stored at -515 80°C. 516 517 Bulk RNA Barcoding and sequencing (BRB-seq) 518 The BRB-seq is a technique for multiplexed RNA-seq 16 which is able to provide high-quality 519 3' transcriptomic data at a low cost (e.g. 10-fold lower than Illumina Truseq Stranded mRNA-520 seq). The data (fastq files) generated from BRB-seq are multiplexed and asymmetrical paired 521 reads. Read R1 contains a 6 bp sample barcode, while read R2 contains the fragment sequence 522 to align to the reference genome. 523 1. Library preparation 524 RNA quantity was assessed using picogreen (Invitrogen P11496). Samples were then grouped 525 according to their concentration in 96-well plates and diluted to a final concentration 526 determined by the lowest sample concentration on the plate. RNA was then used for gene 527 expression profiling using BRB-seq. In short, the BRB-seq protocol starts with oligo-dT 528 barcoding, without TSO for the first-strand synthesis (reverse transcription), performed on each 529 sample separately. Then all samples are pooled together, after which the second-strand is 530 synthesized using DNA PolII Nick translation. The sequencing library is then prepared using 531 cDNA augmented by an in-house produced Tn5 transposase preloaded with the same adapters 532 (Tn5-B/B), and further enriched by limited-cycle PCR with Illumina compatible adapters. 533 Libraries are then size-selected (200 -1000 bp), profiled using High Sensitivity NGS Fragment 534 Analysis Kit (Advanced Analytical, #DNF-474), and measured using Qubit dsDNA HS Assay 535 Kit (Invitrogen, #Q32851). In total, we generated four libraries. For details of library 536 information, please check Table S20. 537

Sequencing 538
Libraries were mixed in equimolar quantities and were then sequenced on an Illumina Hi-Seq polyA sequences of the demultiplexed files by using the "Trim" tool. Next, the STAR aligner 547 31 was used to map the trimmed reads to the reference genome of fly Drosophila melanogaster 548 (BDGP6, Ensembl release 91 32 ). Finally, the read count of each gene was obtained with HTSeq 549 33 . 550

Filtering samples and genes 551
Low-quality samples need to be filtered out, since they might bias results of downstream 552 analyses. In order to assess sample quality, we calculated the number of uniquely mapped reads 553 and of expressed genes for each sample 34 . We removed samples with <0.3 million uniquely 554 mapped reads or with <4500 expressed genes ( Figure S1). We confirmed that these filtered 555 samples are indeed outliers in a multidimensional scaling analysis (MDS) ( Figure S15). Since 556 lowly expressed genes have larger technical error, to minimize the technical noise, we need 557 to remove lowly expressed genes as well. We first calculated counts per million (cpm) with 558 the edgeR package 35 for each gene. Then we removed genes with mean cpm across samples 559

Multidimensional scaling analysis (MDS) 573
A number of factors could be invoked to explain the two groups observed in our 574 multidimensional scaling analysis (MDS) ( Figure 1B), but they should also explain that only 575 one group is structured according to developmental time. The obvious hypothesis that they 576 correspond to male and female embryos does not explain that structure, and is also not 577 supported by examining X/autosome gene expression ratios ( Figure S16). An alternative 578 hypothesis is that samples in the small cluster are unfertilized eggs. If an egg is not fertilized, 579 after completion of meiosis, development will be arrested 39 , but they are visually 580 indistinguishable. This hypothesis is confirmed by at least two lines of evidence, in addition 581 to the lack of developmental time structure. First, for expression correlation, all samples in 582 the small cluster are highly correlated with unfertilized egg, while the correlations in the 583 other samples gradually decrease with development ( Figure S3A). Second, all the samples 584 from the small cluster are enriched with meiosis related genes ( Figure S3B). Thus we 585 excluded the small cluster for downstream analyses, i.e. we used 150 embryos with an 586 average of 18 individuals per developmental stage. 587 588

Metrics of expression variability 589
Expression variability is generally measured by the coefficient of variation (CV) 40 . However, 590 a gene's CV is strongly dependent on its RNA abundance ( Figure S4). While this is an inherent 591 property of time-interval counting processes (such as a Poisson process), it makes the 592 comparison of variability between different conditions difficult 38,41 . Distance to median (DM,593 the distance between the squared CV of a gene and its running median) has been proposed as 594 a variability metric that is independent of expression level 38,41,42 . However, the DM is still 595 strongly negatively correlated with the mean expression level in our data ( Figure S5). To avoid 596 this dependency, we developed another variability measure, the adjusted standard deviation 597