Epigenetic regulation underlying Plasmodium berghei gene expression during its developmental transition from host to vector

Epigenetic regulation of gene expression is an important attribute in the survival and adaptation of the malaria parasite Plasmodium in its human host. Our understanding of epigenetic regulation of gene expression in Plasmodium developmental stages beyond asexual replication in the mammalian host is sparse. We used chromatin immune-precipitation (ChIP) and RNA sequencing to create an epigenetic and transcriptomic map of the murine parasite Plasmodium berghei development from asexual blood stages to male and female gametocytes, and finally, to ookinetes. We show that heterochromatin 1 (HP1) almost exclusively associates with variantly expressed gene families at subtelomeric regions and remains stable across stages and various parasite lines. Variant expression based on heterochromatic silencing is observed only in very few genes. In contrast, the active histone mark histone 3 Lysine 9 acetylation (H3K9ac) is found between heterochromatin boundaries and occurs as a sharp peak around the start codon for ribosomal protein genes. H3K9ac occupancy positively correlates with gene transcripts in asexual blood stages, male gametocytes and ookinetes. Interestingly, H3K9ac occupancy does not correlate with transcript abundance in female gametocytes. Finally, we identify novel DNA motifs upstream of ookinete-specific genes thought to be involved in transcriptional activation upon fertilization.


INTRODUCTION 32
Malaria is caused by apicomplexan parasites of the genus Plasmodium and is transmitted to humans 33 through bites of anopheline mosquitoes. Malaria clinical cases and deaths decreased significantly 34 over the past decade but began to stabilize since 2015 indicating that current measures have now 35 reached their maximum capacity and that new measures are urgently needed (1). Parasite 36 transmission through the mosquito vector is a natural bottleneck in its development and a favourable 37 stage for interventions aiming at malaria control and elimination. Therefore, research towards 38 understanding parasite development in the mosquito has been intensified in recent years. 39 Haploid parasites cyclically infect and asexually replicate in the infected red blood cells (iRBCs) of 40 the host causing disease pathology. Eventually, some parasites escape this cycle and form sexual 41 forms called gametocytes, the stage infective to mosquitoes. In the mosquito gut lumen, ingested 42 gametocytes are activated to form gametes. Female gametocytes exit iRBCs, and messenger RNAs 43 (mRNAs), stored in a messenger ribonucleoprotein (mRNP) complex, become available for 44 translation (2, 3). Activation of male gametocytes involves three rapid rounds of endomitosis leading 45 to eight flagellated microgametes, a process called exflagellation (4). After fertilization, the zygote 46 embarks on a meiotic endoreplication cycle while traversing the mosquito midgut epithelium in the 47 form of an ookinete that upon arrival at the midgut basal side transforms into an oocyst (5). Over the 48 next days, endomitotic replication in the oocyst produces hundreds of sporozoites that, upon oocyst 49 rupture, travel to the mosquito salivary gland for inoculation into a vertebrate host with the next 50 mosquito bite. 51 Epigenetic regulation is important for parasite survival within the human host (6). Genes involved in 52 host-parasite interactions or coding for virulence factors or ligands involved in RBC invasion are 53 epigenetically regulated (7, 8), while some genes involved in drug resistance are epigenetically 54 switched on or off in an environment-dependent manner (9). Importantly, expression of parasite 55 antigens in the mouse host is reset after its passage through the mosquito, suggesting that 56 epigenetic imprinting may be erased during the mosquito stage (10). 57 P. falciparum HP1 (PfHP1) is shown to bind exclusively to histone 3 lysine 9 tri-methylation 58 (H3K9me3) and be the hallmark of transcriptionally silent heterochromatin (7, 11). Heterochromatic 59 loci are largely confined to telomeric and subtelomeric regions as well as chromosome-central 60 islands and almost invariably associated with variantly expressed multigene protein families (7, 11-61 15). In contrast, the universal histone mark linked to euchromatin H3K9ac is shown to be related to 62 active transcription in P. falciparum asexual blood stages (12,14). 63 Asexual blood stage parasites were harvested via heart puncture and passed through a Plasmodipur 97 filter (Europroxima) to remove leucocytes, resuspended in RPMI-1640 medium (Sigma-Aldrich) and 98 crosslinked with 1% formaldehyde in PBS for 10min at 37⁰C. Crosslinking was quenched adding 99 glycine to an end concentration of 0.125M. RBC were then lysed with 0.15% saponin (in PBS) for 5-100 10 minutes on ice. To obtain nuclei the resulting parasite pellet was lysed with cell lysis buffer (20 101 mM Hepes, 10mM KCl,1mM EDTA, 1mM EGTA, 0.65% NP-40, 1mM DTT, 1x protease inhibitor 102 (Roche)). The pellet was resuspended in sonication buffer (1% SDS, 50mM Tris pH8, 10 mM EDTA, 103 1x protease inhibitor (cOmplete™, Mini, EDTA-free, Roche)) and sheared for 25 minutes (30sec ON, 104 30 sec OFF; settings HIGH) using a Bioruptor® Plus sonication device (Diagenode) to obtain DNA 105 fragments of around 100-300bp. 106 The method for the purification of gametocytes was modified from (27). Briefly, mice were pretreated 107 by intra peritoneal injection of 0.2 ml phenylhydrazine (6 mg/ml in PBS) to stimulate reticulocyte 108 formation two days prior to infection with parasites. Gametocyte-enriched blood was harvested via 109 heart puncture and blood was immediately resuspended in 4⁰C coelenterazine loading buffer ( H3K9ac and 6x PbHP1); ookinetes (2x H3K9ac and 5x PbHP1)). 138 For each stage, the following numbers of biological replicates were pooled to obtain enough material: 139 ABS (1); FG (5), MG(5), OOK (2). 140

ChIP sequencing 141
For each sequencing library up to 10 ng of ChIP or input DNA were end-repaired, extended with 30 142 A-overhangs and ligated to barcoded NextFlex adapters (Bio Scientific) as described previously (28 The bam files containing mapped reads from each ChIP was normalised against its input using 158 bamcompare in DeepTools2 (32). For PbHP1, the default settings were used with the following 159 changes: 100bp bin size, 0.01 pseudocount and 1000bp smoothing. 160 The mean of the log2 PbHP1/input ratio was then extracted for each gene, making each gene fit 161 500bp using one 500bp bin using computematrix in DeepTools2 (32) with default settings with the 162 following changes: missing values = 0 and skip 0. The resulting bigwig files were hierarchically 163 clustered with average linkage and Euclidean distance as similarity metric in Cluster3.0 (33). The 164 resulting heatmap and tree was inspected in Treeview (34). 165 The multigene family list is based on with (****) P ≤ 0.0001. 180

Motif identification 181
5' UTRs were extracted from PbANKA genome version 3. Sequences between annotated genes 182 were attributed to the nearest gene as follows. For head-to-head genes, the 5' UTRs were split 1:1. 183 For tail-to-head genes, the sequence was split 1:2 (1/3 was assigned as 3' downstream region to 184 the gene ending and 2/3 was assigned as 5' UTR to the gene starting). To find motifs enriched in 185 ookinete-expressed genes, we used DREME (v 5.0.0) (38) with default settings locally on the MEME 186 suite (39). 187 We used the build-in Gene Ontology tool of PlasmoDB with default settings to identify enriched GO 188 terms (40). 189

RNA preparation 190
Asexual blood stage parasites were harvested via heart puncture and resuspended in 15ml  HEPES and passed through a Plasmodipur filter (Europroxima) to remove white blood cells. Sheep Anti-Mouse IgG) and resuspened in 500ul TRIzol (Invitrogen) and stored at -80⁰C. 213 Genomic DNA was used as a control for the library preparation protocol. For this, asexual blood 214 stage parasites were harvested via heart puncture and resuspended in 15ml RPMI-HEPES and 215 passed through a Plasmodipur filter (Europroxima) to remove white blood cells. Parasite/RBC pellet 216 was lysed with 10ml RBC-lysis buffer (150mM NH4Cl, 10mM KHCO3, 1mM EDTA) for 20min on ice. 217 Parasite-pellets were washed once in ice-cold PBS and lysed in buffer A (500mM NaAc, 100mM 218 NaCl, 1mM EDTA, pH5.2) and 3% SDS. DNA was extracted using phenol:chloroform and sheared 219 to 150bp fragments on a Bioruptor® Plus sonication device (Diagenode). 220

RNA sequencing 221
RNA was extracted using the Direct-zol™ RNA MiniPrep kit (Zymo). Residual gDNA was digested 222 with the TurboDNA-free kit (Ambion). Stranded RNA sequencing libraries were prepared using the 223 RNA HyperPrep Kit (KAPA) following manufacturer's protocol with the exception of the amplification 224 step set to 60⁰C. The RNA library was sequenced in 43bp paired-end reads on a Illumina NextSeq 225 500 sequencer and each sample was split over four non-independent lanes. 226 RNAseq was performed using biological triplicates for each condition. 227 We used sheared genomic DNA (gDNA) as a technical control and did not observe a bias towards 228 GC-rich(er) sequences ( Figure S2A), a known problem in Plasmodium falciparum NGS (14). 229

RNA sequencing analysis 230
Quality of the reads were checked by eye using FASTQC. Reads were aligned separately for each 231 lane to PbANKA genome v3 using HiSat2 with default settings (v 2.0.5.2; --max intron size 5000; --232 fr) (41). Mapped reads were pooled into one bam file per condition and replicate. Importantly, to 233 adjust for potential GC-bias gDNA was processed alongside RNA samples as a control for the RNA 234 library preparation. GC-bias was calculated according to (42) using the default settings for 235 computeGCBias in DeepTools2 suite (32). 236 The PCA plot was calculated using default settings in the DeepTools2 suite (32). 237 We calculated the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value for 238 each gene using featurecounts with default settings (43). 239 For comparative transcriptomics we used Deseq2 (v2.11.39) with default settings (44). We 240 considered significantly differently regulated genes as having a p value of lower than 0.001 (adjusted 241 for multiple testing with the Benjamini-Hochberg procedure which controls false discovery rate 242 (FDR)). 243 All of the analysis has been performed using usegalaxy.org (an open source, web-based platform 244 for data intensive biomedical research) (45) unless stated otherwise. Additionally, usegalaxy.eu was 245 used for some calculations. 246

Generation of epigenetic and transcriptional profiles 249
ChIP was carried out on P. berghei mixed ABS, FG, MG and OOK using antibodies against P. from each other, respectively, suggesting that euchromatin and heterochromatin occupancies differ 257 from each other but are similar between the three stages and parasite lines ( Figure 1A). 258 To investigate the link between epigenetic profiles and gene expression, we complemented the 259 epigenetic profiles with transcriptomics data from the same stages and parasite lines by RNAseq. 260 PCA analysis showed that transcription profiles of ABS and MG are more similar to each other than 261 to the other two stages, and that transcription profiles of FG are more related to OOK 262 (Supplementary Figure 1). 263

Mutually exclusive profiles of heterochromatin and euchromatin 264
Bar plots of ChIPped vs input chromatin showed that heterochromatin is confined to telomeric and 265 subtelomeric regions of all chromosomes, except from the right arm of chromosome 4, in contrast to 266 H3K9ac that is detected in all other chromosomal regions ( Figure 1B). Ten subtelomeric regions 267 could only be partially mapped due to their repetitive nature. As shown in the example of 268 chromosome 8 (Figure 1C) The second chromosome-central heterochromatic region corresponds to cap380 and appears to be 284 epigenetically silenced in ABS, FG and MG but not in OOK ( Figure 1E). Indeed, OOK display 285 increased levels of relative cap380 transcripts, which coincide with a sharp H3K9ac peak within the 286 cap380 5'UTR. Cap380 transcription is controlled by AP2-O (23, 35) and the protein is expressed in 287 early-stage oocysts and localizes to the oocyst capsule (49, 50). Our data provide evidence that 288 heterochromatic silencing is an additional regulatory level of cap380 expression, presumably 289 preventing its premature transcription. 290 A third region that showed clear heterochromatic marks by visual inspection but could not be 291 classified as heterochromatic by our analysis algorithm encompasses the gene encoding AP2-G 292 (PbANKA_1437500; Figure 1F). Expression of P. falciparum AP2-G is epigenetically controlled by 293 PfHP1 (51). The weak levels of PbHP1 marking of this region in our analysis could be explained by 294 our choice of the ANKA 2.33 parasite line. This line carries a mutation in the ap2-g gene resulting in 295 expression of a truncated, non-functional protein unable to induce gametocytogenesis (21), thereby 296 making its epigenetic silencing redundant. 297

Epigenetic silencing of subtelomeric multigene families 298
Our analysis revealed that as many as 221 genes (including pseudogenes) located in subtelomeric 299 regions are significantly enriched in PbHP1 binding ( Figure 1G,  is present in all Plasmodium species, the fam-b family is specific to rodent malaria parasites. 315 Eleven of the heterochromatic, subtelomeric genes belong to the putative reticulocyte binding protein 316 family (aka Pb235) (55), and three of them encode unknown Plasmodium exported proteins. Finally, 317 erythrocyte membrane antigen 1, erythrocyte membrane associated protein 1 and a tryptophan-rich 318 protein pseudogene are also heterochromatic. 319 Only seven subtelomeric heterochromatic genes do not belong to RMP multigene families ( Figure  320 1H).  heterochromatic dynein gene in this study but not orthologous to neither of two heterochromatic 330 dynein heavy chain-encoding genes found in P. falciparum blood stages (Flueck et al., 2009). Our 331 data show that the gene is transcribed in MG, and its promoter is rich in H3K9ac binding in 332 gametocytes ( Figure S2). Its knockout is associated with slow ABS growth (57). It is interesting that 333 of six dynein genes located close to telomeres, only PBANKA_0601200 is found to be 334 heterochromatic in our study. 335 In summary, our data show that similar to other Plasmodium species (15), subtelomeric multigene 336 families are epigenetically silenced via PbHP1 throughout P. berghei development. Interestingly, and 337 in contrast to P. falciparum (15), we did not observe any expansion of heterochromatic boundaries 338 throughout the P. berghei lifecycle. 339

Variant developmental expression of heterochromatic genes 340
Stochastic changes in heterochromatin distribution result in clonal variant gene expression, which 341 forms the basis of P. falciparum antigenic variation (6). Since different parasite lines were used in 342 our study (except FG and MG that derived from the same line), we could not directly determine 343 whether differences in heterochromatin distribution are due to the line epigenetic background or 344 developmental stage, nonetheless, comparative heterochromatic profiling would capture the sum of 345 both. We compared the heterochromatic profiles of genes across development, and found a 346 surprisingly small number of 16 genes with differential heterochromatin occupancy between stages 347 (Figure 2A and Supplementary Figure S3). To identify if variegated clonal expression is present in P. berghei stages, we arbitrarily selected all 356 heterochromatic genes with more than 10 FPKM (mean of three replicates) in at least one 357 developmental stage. Using this approach, we identified an additional 35 genes in ABS, 17 genes in 358 FG, 37 genes in MG and 8 genes in OOK, respectively, which are likely variantly expressed ( Figure  359 2B). Most of these heterochromatic genes display additional euchromatic H3K9ac marks, mainly in 360 their 5'UTR and some of them show very high transcript levels ( Figure 2B). For example, despite 361 the PbANKA_0300600 ORF being heterochromatic, its 5'UTR is enriched in H3K9ac and the gene 362 is highly expressed in all stages ( Figure 2C). This suggests that PbANKA_0300600 is expressed in 363 only a subset of cells. Similarly, the epsilon-tubulin ORF is heterochromatic in all stages but shows 364 euchromatic traits and high transcript levels in both MG and FG ( Figure 2D). The function of epsilon-365 tubulin in Plasmodium is unknown, but the protein is known to mark the older of the two human 366 centrioles upon centrosome duplication (59) and be an essential part of the basal bodies in 367 Tetrahymena (60). 368 Another example of clonally variegated expression includes multigene family members on the right 369 arm of chromosome 7, displaying both heterochromatin and euchromatin occupancies and high 370 transcript levels ( Figure 2E). This is consistent with the finding that some heterochromatic genes 371 are transcribed by a subset of cells, effectively displaying variegated expression as shown before for 372 members of putative exported protein families in P. berghei (Fougère et  Genes displaying hallmarks of clonally variant expression are all located at the heterochromatin-374 euchromatin boundaries, which seems to facilitate clonal variation. 375

Distinct H3K9ac distribution in ribosomal protein genes 376
H3K9ac is a universal histone mark associated with active promoters across the animal kingdom, 377 including P. falciparum ABS, oocysts and sporozoites (12, 61, 62). We investigated the relationship 378 between H3K9ac distribution and gene transcription across the three P. berghei developmental 379 stages. H3K9ac enrichment in the gene ORF, 1 kb 5'UTR and 500 bp 3'UTR was examined against 380 transcript abundance for each stage ( Figure 3A). In consistence with previous findings from P. 381 falciparum, a positive correlation was detected between the gene 5'UTR H3K9ac enrichment and 382 transcriptional levels in ABS (12, 13). We also found that H3K9ac enrichment in the 5'UTR positively 383 correlates with transcription levels in MG and OOK but less so in FG ( Figure 3B). 384 A shift of H3K9ac occupancy toward the gene ORF was observed for highly expressed genes, mainly 385 in ABS and OOK (see arrowheads in Figure 3A). Many of the genes with high expression encode 386 ribosomal proteins. Closer investigation indicated that H3K9ac occupancy sharply peaks around the 387 start codon of ribosomal protein genes in ABS, FG and OOK ( Figure 3C). Intriguingly, H3K9ac 388 occupancy was extended further into the 5'UTR of MG resulting in no clear peak being detected 389 around the start codon ( Figure 3D). Further analysis confirmed that the mean H3K9ac enrichment 390 in the 500 bp 5'UTR of ribosomal genes was significantly different between MG and all other stages 391 (Wilcoxon signed rank test, p≤0.0001), but no difference was detected between any of the other 392 stages ( Figure 3E). As a control, the same analysis using the 500 bp 3'UTR of these genes did not 393 show any significant differences between any of the stages (Wilcoxon signed rank test, p>0.3) 394 ( Figure 3F). Interestingly, relative transcripts of ribosomal protein genes differ between all 395 developmental stages, with highest transcript numbers found in ABS ( Figure 3G). We did not detect 396 a correlation between transcriptional abundance and H3K9ac occupancy. These data clearly show 397 that genes encoding ribosomal proteins exhibit a different H3K9ac pattern from other genes, 398 highlighting interesting differences in epigenetic makeup for a specific gene group. It has been 399 previously suggested that a shift in peak shape can indicate a different function of a gene (63), and 400 highly dense and narrow distributions of H3K9ac near transcriptional start sites have been 401 associated with constitutive expression of genes involved in translation in plants (64). 402

H3K9ac enrichment does not correlate with transcript levels in FG 403
We further wanted to examine the relationship between H3K9ac occupancy and transcript levels in 404 P. berghei development. For this, we identified genes that are differentially expressed between two 405 related developmental stages and examined their H3K9ac occupancy for each of the two stages 406 (Figure 4A, B, C). Genes upregulated in ABS (1426), MG (1523) or OOK (1641) compared to FG 407 exhibit greater H3K9ac enrichment in their 5'UTR (1000bp upstream of the start codon) in each 408 respective stage than genes downregulated in any of these three stages (1429, 1518 and 1534,  409 respectively) compared to FG (Figure 4D, E, F). However, when the same groups of genes were 410 examined for H3K9ac enrichment in their 5'UTR in FG, upregulated genes in FG showed less (FG 411 vs. ABS and FG vs. MG) or the same (FG vs. OOK) H3K9ac enrichment than downregulated genes 412 ( Figure 4G, H, I). These results indicate that H3K9ac occupancy in the gene 5'UTR is a good 413 predictor for relative transcript levels in ABS, MG and OOK. 414 We further investigated the absence of positive correlation between H3K9ac and transcript levels by 415 generating two high-confidence lists of genes that are either upregulated (591 genes) or 416 downregulated (271 genes) in FG compared to all other stages, respectively ( Figure 5A and  417 Supplementary Table 3). Gene Ontology (GO) analysis using the GO slim version showed that the 418 former list includes genes involved in cell adhesion, reproduction, motility, differentiation, cell cycle 419 and locomotion, all of which are characteristic of cells preparing for fertilization and development into 420 the motile ookinete stage (Figure 5B and Supplementary Table 2). In contrast, downregulated 421 genes in FG are involved in ribosome biogenesis and translation, suggestive of a translationally less 422 active FG stage compared to the other three stages. FG are known to produce and store mRNAs 423 that are translationally repressed (2). Similar to before, these downregulated FG genes show 424 significantly higher H3K9ac enrichment in their 5'UTRs than genes that are upregulated in FG 425 (Figure 5C). 426 These data together strongly indicate that H3K9ac, a histone mark otherwise universally associated 427 with active promoters, does not directly correlate with transcript levels in FG. Instead, it appears that 428 H3K9ac marks the promoters of most genes in FG irrespective of their transcript levels. process. It is important to note that while additional OOK transcription factors have been described 461 (16), none of the published DNA-motifs matches those found here.

DISCUSSION 463
Our study provides new insights into the epigenetic regulation of Plasmodium gene expression 464 during its developmental transition from asexual to sexual and zygotic stages, which occurs as the 465 parasite moves from its vertebrate host to the insect vector and is essential for parasite transmission. 466 Epigenetic silencing has been associated with key strategies of parasite development and 467 environmental adaptations. Our data reveal that HP1-mediated epigenetic silencing remains largely 468 unchanged during this transition and confined to chromosomal subtelomeric regions. Mapping 469 Plasmodium subtelomeric genes is very challenging as most genes belong to multigene families with 470 high levels of sequence identity between them. Indeed, as many as 51 heterochromatic genes could 471 not be mapped to P. berghei chromosomes in our analysis : 14 fam-a, 8 fam-b, 1 fam-c, 27 pir, and  472 one encoding a tryptophan-rich antigen. 473 Comparison of our findings with those of Fraschka et al. (15) that analyzed PbHP1 occupancy in 474 ABS of the parental P. berghei ANKA strain revealed that 177 heterochromatic genes are shared 475 between the two studies, 46 are specific to our study and 14 are specific to the other study 476 (Supplementary Figure S4A). Whereas the vast majority of differences pertain to subtelomeric 477 regions, where our analysis appears to be more sensitive, the combination of the two studies 478 identified four non-subtelomeric heterochromatic genes, of which only cap380 encoding an oocyst 479 capsule protein is detected in both studies. On the one hand, PBANKA_0934600 that encodes a 480 protein of unknown function was only detected in our analysis but a closer visual inspection of the 481 profiles obtained by Fraschka et al. (15) shows extensive occupancy of the gene by PbHP1 482 (Supplementary Figure S4B). On the other hand, ap2-g and ap2-sp3/ap2-tel are not detected as 483 heterochromatic in our analysis but visual inspection indicates increased PbHP1 occupancy 484 (Supplementary Figure S4C). These data suggest that automated detection of moderate 485 heterochromatin enrichment remains a challenge and that conclusions about absence of epigenetic 486 silencing in such regions should be treated with caution. Nonetheless, the two studies together reveal 487 that P. berghei heterochromatin formation, maintenance and inheritance are largely hardwired 488 throughout development. 489 It has been previously proposed that epigenetic reprogramming or resetting of the expression of 490 genes involved in parasite virulence occurs during parasite passage through the vector (10, 66). 491 However, no such evidence for a reset of virulence is detected pre-and post-meiotically in our study, 492 which could be explained by two scenarios. Firstly, such reset is likely to start immediately after 493 nuclear fusion as in higher eukaryotes (67) and thus may have been missed by our study design. 494 Secondly, any such presumed reset may involve different epigenetic marks than those examined 495 here. Finally, the virulence-reset hypothesis is based on parasites serially maintained in rodents, 496 which are known to become more virulent. Therefore, it is probable that the observed phenotype is 497 based on ill-managed heterochromatin maintenance undetectable by our study. 498 Gametocytes are poised cells that are rapidly activated and transform into respective gametes upon 499 ingestion by a mosquito. In P. falciparum, heterochromatin boundaries are shown to expand outside 500 telomeres as the parasite lifecycle progresses from ABS to gametocytes (15). However, our data do 501 not show any major differences in heterochromatin distribution between the various P. berghei 502 lifecycle stages examined. This result can be explained by the fact that genes that become 503 heterochromatic in P. falciparum gametocytes and are involved iRBC remodeling (such as knob 504 formation) have no clear orthologues in P. berghei and highlights differences of heterochromatin 505 maintenance between the two parasites. 506 Male and female gametocytes are characterized by specific transcriptomes and proteomes (22, 68-507 70). However, in our analysis, differential heterochromatin occupancy between these two stages is 508 limited to 12 genes, suggesting that epigenetic gene silencing is not a key regulatory mechanism of 509 differential gene expression between male and female gametocytes. In addition, all these genes are 510 located in subtelomeric regions and 11 of them belong to a multigene family, and thus some of the 511 detected differences may be due to the method's sensitivity. It remains to be seen whether any of 512 the 12 genes show true expression differences between male and female gametocytes and relate to 513 sex-specific functions. 514 A striking difference between male and female gametocytes is the level of H3K9ac occupancy in the 515 5'UTR of genes, which correlates with transcript abundance in male but not in female gametocytes. 516 Thus it appears that gene transcription in female gametocytes is independent of H3K9ac levels. We 517 propose that the female gametocyte is a stage of epigenetic remodeling that involves universal 518 marking of gene promoters with H3K9ac. At the same time, genes that are involved in ribosome 519 biogenesis and translation are downregulated, a phenomenon that together with translational mRNA 520 repression by the DOZI-complex suggests that a multilayered mechanism regulating zygotic 521 development operates in Plasmodium. 522 Mature oocytes in mammals, flies and worms are transcriptionally silent when passing through 523 meiosis I preceding fertilization (71). Additionally, early embryonic development depends solely on 524 maternally deposited proteins and RNAs and coincides with low or undetectable transcription (72). 525 Zygotic transcription often resumes after several hours or days post fertilization, a process called 526 maternal-to-zygotic transition accompanied by zygotic genome activation, and in flies, for example, 527 is controlled by a master transcription factor that is not necessarily associated with H3K9ac (73). 528 Additionally, epigenetic reprogramming occurs before and after fertilization, ensuring the transition 529 from a highly specific cell type back to a totipotency state able to form a new organism. However, 530 direct comparison of malaria parasite female gametocytes to oocytes of multicellular organisms is 531 complicated by the fact that malaria parasites undergo meiosis directly after fertilization (74)