RNA-seq identifies rapid changes in gene expression during the schizont-to-ring transition

While gene expression in the main phases of the IDC has been studied in detail, limited insight into gene expression in merozoites and how it relates to the other stages is available [7, 20]. Given that the transition from one erythrocyte to the next requires numerous changes in parasite biology, we were interested in characterizing transcriptional changes that occur in parasites during merozoite formation. To this end, we performed RNA-seq using three biological replicates of synchronized early schizonts (40 h post infection or hpi) and merozoites collected as outlined in Additional file 1: Fig. S1A. The 40 hpi schizont time point was selected because parasites collected at a later time point (44–48 hpi) may already include fully segmented schizonts, due to the challenges of fully synchronizing P. falciparum cultures. Merozoites were prepared by treating early schizonts (collected at 40 hpi) with the egress inhibitor E64. E64-treated parasites progress through schizogony normally, but are blocked at the final step of egress, just prior to rupturing of the erythrocyte membrane. Fully segmented schizonts were then mechanically ruptured to release merozoites.

To interrogate differences in gene expression between early schizonts (40 hpi) and merozoites, we performed differential gene expression analysis with these two parasite populations. We observed two sets of differentially expressed genes between these populations, demonstrating that as parasites transition from schizonts to merozoites, the transcriptional landscape is altered (Fig. 1A, Additional file 2: Tables S1, S2). Notably, a group of genes appear to be turned on as parasites transition to merozoites, countering the idea that this stage is transcriptionally silent. Unexpectedly, we observed prominent heterogeneity in our merozoite samples. To probe whether this was the result of technical issues in our sample preparation and sequencing or reflective of true biological differences between replicates, we took advantage of data available from a recent study by Wichers et al. [33]. This study generated RNA-seq data at several time points across the IDC, including early schizonts (40 hpi) and merozoites. For the genes that were upregulated in merozoites as compared to schizonts in our data set, the merozoite samples from Wichers et al. also displayed prominent variability (Fig. 1A, Additional file 2: Table S2). Importantly, the trends in expression observed in our data, including upregulation of a subset of genes, were also present in the Wichers et al. data set, demonstrating that despite heterogeneity within timepoints, both data sets captured similar transcriptional changes between early schizonts and merozoites.

Fig. 1

Gene expression changes rapidly during transitions between schizont, merozoite, and ring stage parasites. A Heatmap depicting differentially expressed genes (log2(fold change) > 1.5 or < -1.5 and adjusted p-value < 0.1) between early schizont samples (40 hpi, n = 3) and merozoites (MZ, n = 3), all collected by us. Expression of these genes in early schizont (40 hpi, n = 3) and merozoite (MZ, n = 3) samples generated by Wichers et al. are shown on the right. B Principal component analysis of all RNA-seq samples. Merozoite subpopulations defined in panel C are indicated by blue numbers. Amount of variance between samples accounted for by each component is shown on the x- and y- axes. C Heatmap depicting differentially expressed genes (log2(fold change) > 1.5 or < -1.5 and adjusted p-value < 0.1) between early schizont samples (40 hpi, n = 6) and merozoites (n = 6). Genes are divided into 4 clusters based on expression pattern by k-means clustering as indicated on the left of the heatmap. Merozoite subpopulations are indicated by the labels on the bottom of the heatmap. Enriched gene ontology (GO) terms in each cluster are displayed to the right of the heatmap. D Boxplots of average gene expression for each group as shown in panel C. Gene expression differed between groups (Kruskal–Wallis test). P values indicated in the graph are from Dunn’s post hoc tests. *, P < 0.05; **, P < 0.01. For all boxplots, box edges indicate the first and third quartiles, mid-line indicates the median, and whiskers indicate the minimum and maximum data values. E H3K9ac and H3K27ac expression shown as log2(ChIP/input) around the translation start site (ATG) for genes from each group as indicated in panel C in early schizonts (40 hpi), merozoites, and early rings (4 hpi). Graphs were generated using computeMatrix from the deepTools package

Given the heterogeneity in merozoite samples and the similarities between the two data sets, we chose to combine our data with that from Wichers et al. for further analysis. In addition to the early schizont and merozoite samples, Wichers et al. also performed RNA-seq on mid schizonts (44 hpi, n = 3), late schizonts (48 hpi, n = 3) and mid rings (8 hpi, n = 3). By including these time points in our analysis, we were able to interrogate changes in transcription at higher resolution across the schizont-to-merozoite transition as well as characterize alterations in gene expression between merozoites and mid rings.

We first performed principal component analysis (PCA) to understand global differences between parasite stages. This analysis reaffirmed the heterogeneity we initially observed in merozoites (Fig. 1B). The first component of the PCA mainly captured differences between the stages, while the second component captured variability between the two data sets (Additional file 1: Fig. S2). Generally, schizont stages (early schizont, 40 hpi; mid schizont, 44 hpi; late schizont, 48 hpi) clustered together and mid ring stage parasites (8 hpi) clustered separately. The merozoite population had the most spread, with some samples showing similarity to early schizonts (subpopulation 1), some split between late schizonts and mid rings (subpopulation 2), and the final set located near the mid rings (subpopulation 3) (Fig. 1B). This suggests that within the six merozoite samples, three subpopulations were captured.

To characterize the differences in gene expression across the schizont-to-ring transition in more detail, we performed differential gene expression analysis on the combined data set. Since we observed significant differences in expression between early schizonts and merozoites in the initial analysis of only our samples (Fig. 1A) we sought to identify genes that were differentially expressed between early schizonts (40 hpi, n = 6) and merozoites (n = 6). Based on a magnitude log2(fold change) > 1.5 and an adjusted p-value < 0.1, we identified 936 differentially expressed genes (Additional file 2: Table S3). We then assessed the expression of these genes in early schizonts (40 hpi), mid schizonts (44 hpi), late schizonts (48 hpi), merozoites, and mid rings (8 hpi) (Fig. 1C, Additional file 2: Table S4). Genes were grouped through k-means clustering based on expression across the stages into 4 clusters (Fig. 1C). Genes expressed from early to late schizonts (n = 110) in our differential gene expression analysis were captured in cluster 1 (Fig. 1C, D). As parasites progressed through the developmental cycle from early schizonts to mid rings, expression of these genes decreased. This cluster was enriched for genes related to the apicoplast as determined by gene ontology (GO) analysis (Additional file 2: Table S5). Genes with increased expression in late schizonts and merozoites (cluster 2, n = 205) were enriched for GO terms relating to actin filament organization and cytoskeleton organization (Additional file 2: Table S5), which are critical processes required for successful invasion (reviewed in [34]). These transcripts are likely rapidly downregulated following invasion, as expression levels in mid rings were comparable to those seen in schizonts (Fig. 1D).

Cluster 3 (n = 450) was comprised of genes expressed at elevated levels in mid rings and enriched for GO terms related to ribosome biogenesis and rRNA processing (Fig. 1D, Additional file 2: Table S5). In merozoites, expression of these genes was heterogeneous. In two merozoite samples (subpopulation 1 in Fig. 1B, C), we observed low levels of expression comparable to that seen in early, mid, and late schizonts. Two merozoite samples (subpopulation 2 in Fig. 1B, C). had expression levels intermediate between that seen in schizonts and that seen in mid rings. In the final two merozoite samples (subpopulation 3 in Fig. 1B, C), we observed high levels of these transcripts, comparable to or higher than the levels observed in mid rings. We hypothesize that heterogeneity in the merozoite samples reflects three different stages along the merozoite developmental trajectory, with expression of genes related to ribosome activity and RNA processing increasing as merozoites mature. This gene expression program appears to reach its highest level in late stage merozoites (subpopulation 3) and continue during ring development. Interestingly, differential gene expression analysis between merozoites and mid rings (8 hpi) revealed that expression of genes involved in DNA replication and mitosis gradually decreased in merozoites (Additional file 1: Fig. S3A, Additional file 2: Tables S3, S6, S7, cluster A, n = 1,253). This analysis also identified a second set of genes enriched for GO terms related to translation and ribosome biogenesis that was upregulated in mid rings (8 hpi) as compared to all other stages (Additional file 1: Fig. S3A, Additional file 2: Tables S3, S6, S7, cluster C, n = 1,057). Combined, our differential gene expression analyses suggests that, prior to invasion, parasites downregulate genes involved in schizogony while upregulating genes required for translation and RNA processing. After invasion, another wave of transcription results in the expression of additional genes involved in translation and metabolism that are likely necessary to further facilitate ring stage development in a new erythrocyte.

Genes in the final cluster of differentially expressed genes between early schizonts and merozoites (cluster 4, n = 171) were elevated in merozoites and mid rings; however, expression of this cluster was more variable than what we observed in the other three clusters (Fig. 1C, D). This cluster was enriched for genes with roles in lipase activity (Additional file 2: Table S7). Interestingly, it appears these genes were consistently elevated in samples from Wichers et al. compared to our samples, suggesting variation in this cluster may be due to stochastic differences in parasite populations. Consistent with this possibility, genes known for variable expression within parasite populations, such as rifins, were present within this cluster.

Together, these data demonstrate that parasites undergo rapid changes in gene expression during or following segmentation into merozoites, resulting in heterogeneity among merozoite populations.

Changes in the distribution of H3K9ac and H3K27ac accompany gene expression changes in merozoites

Given the changes in gene expression as parasites transition from early schizonts to mid rings, we were interested in investigating potential mechanisms behind this shift. Epigenetic regulation of gene expression by the distribution of histone PTMs is known to influence stage-specific expression of genes during the IDC [25]. To uncover differences in the histone PTM landscape that could be related to the differences in gene expression we observed in our RNA-seq analysis, we performed ChIP-seq for the histone PTMs known to be associated with active gene expression (H3K9ac, H3K27ac, and H3K4me3) on early schizonts (40 hpi), merozoites, and early rings (4 hpi). Samples for ChIP-seq and RNA-seq were collected separately as outlined in Additional file 1: Fig. S1A, B. Multiple samples from each stage were pooled to control for heterogeneity; however, this means our ChIP-seq data is reflective of the average histone PTM profile across each stage and does not capture differences in histone PTMs that may exist between the three subpopulations of merozoites identified in our RNA-seq analysis.

First, we assessed the similarity of global histone PTM distribution between marks and stages. Generally, within and between stages, the distributions of H3K4me3, H3K9ac, and H3K27ac showed moderate to high correlations (Additional file 1: Fig. S4, blue box). This indicates that within each stage these marks likely colocalize and that, generally, these marks have similar distribution patterns between stages. However, for each mark, the correlation between early schizont and early ring stages was higher than between either of the two stages and merozoites. Additionally, the correlation between H3K9ac and H3K27ac was stronger in merozoites (r = 0.92) than the other two stages (early schizonts, r = 0.56; early rings, r = 0.43). These results suggest that the histone PTM landscape is dynamically remodeled during the schizont-to-ring transition, with a distinct distribution in merozoites.

Given the higher correlation between H3K9ac and H3K27ac in merozoites (compared to the other stages) and their well-known association with active gene expression, we investigated how the distribution of these marks changed during the schizont-to-ring transition. Globally, H3K9ac and H3K27ac followed similar distribution patterns in early schizonts, merozoites, and early rings (Additional file 1: Fig. S5). To probe how differences in gene expression relate to differences in histone PTM distribution, we looked at the enrichment of H3K9ac and H3K27ac for the four clusters of differentially expressed genes identified in our RNA-seq analysis (Fig. 1E). For all groups, both marks had similar distribution patterns with sharp increases at the gene start site and enrichment into the gene body. As expected, genes with higher expression levels in schizonts (cluster 1) had increased H3K9ac in early schizonts. H3K9ac and H3K27ac decreased for this cluster of genes in merozoites, which corresponded with a decrease in transcript levels at this stage. The promoter regions of genes in cluster 3, that were upregulated in merozoites, were enriched in H3K9ac and H3K27ac in merozoites as compared to these marks in the other stages. In rings, the levels of H3K9ac and H3K27ac in these genes reverted back to the levels seen in schizonts, even though the corresponding transcripts were still detected in the mid ring stage. This could be reflective of maintenance of these transcripts when the corresponding genes have already been turned off. Collectively, these results show changes in the distribution of H3K9ac and H3K27ac that are closely related to differences in gene expression as P. falciparum transitions from the schizont to ring stage.

Unique histone PTM profile defined by depletion of H3K4me3 is associated with gene expression in P. falciparum merozoites

In addition to the increase in H3K9ac and H3K27ac in merozoites, we observed differences in H3K4me3 between the stages, primarily in the coding regions (Additional file 1: Fig. S5). While all stages had sharp enrichment of this mark at the gene start and end, H3K4me3 was only enriched in the gene body in merozoites. Surprisingly, we observed similar levels of H3K4me3 in schizonts and rings (Additional file 1: Figs. S5 and S6). Previously, this mark was shown to be present at low levels in ring stages, with peak enrichment during the trophozoite stage [26]. Our results suggest that H3K4me3 is relatively constant between early schizonts, merozoites, and early rings with the most prominent difference being enrichment of H3K4me3 in the gene bodies in merozoites (Additional file 1: Fig. S5). In keeping with this previous study, we did not observe changes in H3K4me3 in relation to gene expression (Additional file 1: Fig. S6). Our results suggest that H3K4me3 is present in coding regions throughout the IDC and is uncoupled from gene expression.

Since H3K9ac and H3K27ac are associated with transcription and H3K4me3 is associated with permissive chromatin, we sought to assess the combinatorial distribution of these marks across the genome. To accomplish this, we performed ChromHMM analysis, a machine learning approach to identifying combinatorial chromatin states [35]. This tool utilizes histone PTM distribution information to develop a model of chromatin states based on the frequency of each mark across the genome in 200 bp segments. Using the ChIP-seq data we generated for H3K9ac, H3K27ac, and H3K4me3 in early schizonts (40 hpi), merozoites, and early rings (4 hpi), we identified 5 unique chromatin states (Fig. 2A, Additional file 2: Tables S8, S9) and annotated the genome with these states to understand the relationship between each chromatin state and specific genomic features (intergenic regions, coding sequences, gene start sites, etc.).

Fig. 2

Unique histone remodeling is associated with a subset of merozoite-expressed genes. A State emissions profile generated by ChromHMM analysis. B Heatmap depicting the frequency of each ChromHMM state in the promoters (1,000 bp region upstream of the gene ATG) of all genes (n = 5,602) in early schizonts (SZ, 40 hpi), merozoites (MZ), and early rings (R, 4 hpi). C H3K27ac, H3K9ac, and H3K4me3 ChIP tracks depicting a region of H3K4me3 depletion upstream of AP2-EXP (PF3D7_1466400) and corresponding ChromHMM states. Black rectangle indicates the region where an H3K4me3 depletion is flanked by enrichment of H3K9ac and H3K27ac in merozoites. At the bottom, the AP2-EXP gene body is indicated by a green box, while the region with the 1 – 5 – 1 ChromHMM state pattern is marked by the blue box (coordinates on chromosome 14: 2,719,800–2,722,200 bp). Data range of ChIP tracks is indicated in brackets on the right. Scale bar depicting 1,000 bp is shown. D Average ChIP/input of H3K27ac, H3K9ac, and H3K4me3 across 66 regions with H3K4me3 depletion flanked by enrichment of H3K9ac and H3K27ac. Regions were scaled to the same size by computeMatrix, and the enrichment 1.0 kb upstream and downstream of the regions was also graphed. E Heatmap depicting the frequency of each ChromHMM state in the promoter region (1.000 bp upstream of the gene ATG) of schizont-expressed genes, ring-expressed genes, and genes with H3K4me3 depletion in the upstream intergenic region in early schizonts (S, 40 hpi), merozoites (M), and early rings (R, 4 hpi). F Functions associated with genes with H3K4me3 depletion. Functions were assigned based on information available in PlasmoDB and a review of the literature and can be found in Additional file 2: Table S12. Genes annotated as ‘conserved with unknown function’ (n = 9) in PlasmoDB are not shown in the graph but were included in all analyses. G Heatmap depicting normalized gene expression in protein-coding genes with H3K4me3 depletion (n = 70). Four non-coding genes with H3K4me3 depletion are not plotted because these genes were excluded from the RNA-seq differential gene expression analysis. H Boxplots depicting average gene expression of the genes with H3K4me3 depletion (panel) and a set of randomly selected genes (n = 99). Gene expression differed between time points for the panel of genes with H3K4me3 depletion (Kruskal–Wallis test). P values indicated in the graph are from Dunn’s post hoc tests. **, P < 0.01; ***, P < 0.001; ****, P < 0.0001. I) DNA motif in the H3K4me3 depletion regions. J) DNA motif in the H3K4me3 depletion regions of genes associated with PTEX

We were interested in identifying differences in histone PTM landscape between parasite stages that could be related to the regulation of gene expression, so we focused our analysis on the intergenic regions upstream of gene start sites, which are likely to have promoter functions. We determined the abundance of each of the states across all promoter regions (1,000 bp upstream of the gene start site) in early schizont, merozoite, and early ring genomes. As parasites matured from early schizonts to merozoites, promoters already marked during the early schizont stage with H3K9ac and H3K4me3 appeared to gain H3K27ac as demonstrated by a reduction in ChromHMM state 4 and an increase in ChromHMM state 1 in merozoites (Fig. 2B, Additional file 2: Table S10). This may reflect a transition to a state of active gene expression in merozoites as H3K27ac is known to mark active promoter and enhancer regions [36, 37]. In keeping with this theory, genes with the highest expression in merozoites also had the highest enrichment of H3K27ac (Additional file 1: Fig. S6A, B). Following re-invasion, H3K9ac and H3K27ac levels appeared to decrease in gene promoters. These regions became mainly characterized by H3K4me3 and low levels of H3K9ac as seen in the increase in ChromHMM states 3 and 5 and decrease in ChromHMM states 1 and 4 in early rings compared to merozoites (Fig. 2B). This may reflect that the histone PTM landscape is altered following invasion to allow for the stage-specific expression of genes.

While visualizing our ChIP-seq data in IGV, we noticed that some intergenic regions contained a strong depletion of H3K4me3 flanked on either side by enrichment of H3K9ac and H3K27ac in merozoites (representative gene shown in Fig. 2C and average enrichment patterns shown in Fig. 2D). To confirm this finding, we generated a biological replicate for H3K4me3 in merozoites. The two H3K4me3 replicates were strongly correlated (Pearson correlation coefficient, r = 0.91). Importantly, the depletions we originally observed in intergenic regions were also present in the biological replicate (representative genes shown in Additional file 2: Fig. S7A–C). This pattern of H3K4me3 depletion flanked by enrichment of H3K9ac and H3K27ac corresponded with the ChromHMM state pattern 1–5–1 (merozoites, Fig. 2C). Those same regions in early schizonts and early rings showed a slight depletion of H3K4me3 that was not flanked by H3K9ac and H3K27ac and were mainly defined by ChromHMM state 2 (Fig. 2C). Using a custom Perl script, we identified all regions of the genome with H3K4me3 depletion marked by the 1–5–1 ChromHMM state pattern. After filtering for regions located in intergenic areas, we identified 66 regions with the H3K4me3 depletion pattern (Fig. 2D, Additional file 2: Table S11). Across these regions, the distribution of H3K9ac, H3K27ac, and H3K4me3 in merozoites differed dramatically from that in early schizonts and early rings where there was little fluctuation in all three marks.

The 66 areas with the H3K4me3 depletion pattern (i.e., ChromHMM 1–5–1 pattern) were located in intergenic regions upstream of a panel of 74 genes (Additional file 2: Table S11). These regions were found at varying distances from the genes and were not restricted to the 1000 bp upstream of gene start sites that we defined as promoter regions in our ChromHMM analysis (Fig. 2B). However, promoter regions and upstream regulatory elements in Plasmodium are poorly defined [28, 38, 39], meaning that while outside of the region we defined as promoters, these H3K4me3 depletions may still fall in regulatory regions.

To gain a better understanding of how the chromatin landscape of this gene panel compared to that of genes with known stage-specific expression patterns, we analyzed the frequency of ChromHMM states present in the promoter regions (defined as the 1,000 bp preceding the gene start site) of the H3K4me3 depletion genes (Additional file 2: Table S12) as well as a set of schizont-expressed genes (905 genes, Additional file 3: Table S13) and a set of ring-expressed genes (139 genes, Additional file 2: Table S14). As expected, promoters of schizont-expressed genes were enriched for high H3K9ac, high H3K27ac, and medium H3K4me3 in early schizonts (ChromHMM state 1, Fig. 2E, Additional file 2: Table S15). As early schizonts transitioned to merozoites and then early rings, schizont-expressed gene promoters lost this chromatin state, likely reflecting silencing of these genes following the schizont stage. Ring-expressed genes followed the opposite pattern, with promoters defined by ChromHMM state 1 increasing in frequency as parasites transitioned from early schizonts to early rings (Fig. 2E). This likely reflects the activation of ring-expressed gene promoters in the early ring stage. Interestingly, the promoters of the genes with the H3K4me3 depletion (i.e., ChomHMM 1–5–1 pattern) followed a similar pattern as the ring-expressed genes with the frequency of active promoters (characterized by ChromHMM state 1) increasing as parasites transitioned from early schizonts to early rings (Fig. 2E). Additionally, genes with H3K4me3 depletion were enriched for state 5 in merozoites, which was not observed for ring genes. These differences in the enrichment of state 5 between the promoters of the H3K4me3 depletion genes and ring-expressed genes demonstrates that the promoters of the genes with H3K4me3 depletion are subject to epigenetic remodeling not seen in the ring-expressed gene promoters.

Genes with H3K4me3 depletion (i.e., ChromHMM 1–5–1 pattern) in the upstream intergenic region had roles in numerous cellular processes, including egress and invasion, RNA-binding, and metabolism (based on a review of the literature, Fig. 2F, Additional file 2: Table S12). Exported proteins (n = 10) and proteins involved in export pathways (n = 12) were the two largest groups of genes. Protein export in early ring parasites is crucial to facilitate host cell remodeling following invasion [40,41,42]. Together, this suggests that the depletion of H3K4me3 in putative regulatory regions in merozoites is associated with a select group of genes that play important roles in establishing a productive infection in a newly invaded erythrocyte.

Given that the chromatin state of the H3K4me3 depletion gene promoters bore similarity to the promoters of ring-expressed genes, we next investigated the timing of expression of these genes. Using our RNA-seq data supplemented with the data from Wichers et al., we assessed the expression of the genes with H3K4me3 depletion in the upstream intergenic regions in early schizonts (40 hpi), mid schizonts (44 hpi), late schizonts (48 hpi), merozoites, and mid rings (8 hpi). The majority of these genes were found cluster 3 from our differential gene expression analysis between early schizonts and merozoites (Fig. 1C). Generally, expression was highest in more mature merozoites (subpopulations 2 and 3, Fig. 2G, H, Additional file 2: Table S16, S17) and mid rings, consistent with roles for these genes in ring stage development following invasion. In keeping with this observation, nearly half of the genes with H3K4me3 depletion in the upstream intergenic region (n = 36) reached peak transcription at 1 hpi based on the data available from Painter et al. [20], which is a significant enrichment as compared to all genes (675/5373, P < 0.0001, Chi-squared test). Combined, these results indicate that genes with H3K4me3 depletion in the upstream intergenic region undergo active transcription in merozoites or immediately after invasion.

The uniform pattern of expression we observed for genes with H3K4me3 depletion (i.e., ChromHMM 1–5–1 pattern) in the upstream intergenic region suggests they are subjected to common regulatory processes. A protein binding to the DNA would result in depletion of H3K4me3 and lack of enrichment for H3K9ac and H3K27ac at the binding location in our ChIP samples much like the pattern we observed. To probe the possibility that a protein bound the region of H3K4me3 depletion, we performed a motif search in these regions (indicated by ChromHMM state 5, Fig. 2C) using MEME [43]. In 70 out of the 74 genes assessed, we identified a common motif (Fig. 2I). This motif bore similarity to the core DNA motif of AP2-EXP (CATGC) as well as the binding motif of AP2-11a [22, 44, 45]. Using Simple Enrichment Analysis [46], we observed that this motif was enriched in the H3K4me3 depletion regions (P = 0.0004) but not in randomly selected intergenic regions (P = 0.8), suggesting that these genes may be subjected to regulation by an AP2 transcription factor. Based on these results, we hypothesize that the observed histone PTM pattern in upstream intergenic regions is the result of a protein binding to this region and shielding the histone PTMs from antibody recognition during ChIP.

H3K4me3 depletion in promoter regions is associated with PTEX and PTEX-related genes

Genes encoding proteins involved in export pathways (n = 12) was the largest group of genes within the set of genes with H3K4me3 depletion in upstream intergenic regions (Fig. 2F). Out of the 12 genes involved in protein export, 7 were part of or known to interact with the Plasmodium translocon of exported proteins (PTEX), a protein complex important for export of proteins from the parasitophorous vacuole into the host cell [41]. The promoters of three components of PTEX (PTEX150, HSP101, and PTEX88) had a similar histone PTM profile, with depletion of H3K4me3 flanked by H3K9ac and H3K27ac in the upstream intergenic regions. However, these regions were annotated primarily with a 1 – 5 – 2 – 5 – 1 ChromHMM state pattern and so they were not identified in our original analysis of 1 – 5 – 1 ChromHMM regions (Additional file 1: Fig. S8A–C). Like the other genes with H3K4me3 depletion in the upstream intergenic regions, these three genes were also expressed at higher levels in merozoites and early rings (Additional file 1: Fig. S8D–F, Additional file 2: Table S3, S4). Additionally, 9 out of the 10 PTEX-associated genes (including PTEX150, PTEX88, and HSP101) contained a common motif in the H3K4me3 depletion region (Fig. 2J) suggesting potential regulation by a common transcription factor; however, this motif does not match any known Plasmodium transcription factor motifs.

Stage-specific depletion of H3K4me3 in intergenic regions is seen in both hepatic and erythrocytic merozoites

Liver stage Plasmodium development culminates in the formation of hepatic merozoites that will invade erythrocytes to initiate a blood stage infection [5]. The processes of egress and invasion that occur during the parasite transition from the liver to the blood have similarities to those that occur as parasites transition between erythrocytes, thus gene expression in hepatic merozoites and erythrocytic merozoites may be subjected to similar regulatory processes [47, 48]. To probe this possibility, we performed ChIP-seq for H3K4me3 and H3K9ac in Plasmodium berghei hepatic merozoites (68–75 hpi). Due to the technical challenges of producing hepatic P. falciparum merozoites, we utilized P. berghei for this analysis. Merozoites were prepared by collecting merosomes and releasing hepatic merozoites by saponin treatment. Due to limited input material, we only assessed H3K4me3 and H3K9ac.

In P. berghei hepatic merozoites, H3K4me3 and H3K9ac were enriched in intergenic regions and depleted in gene bodies (Fig. 3A). Similar to our observations in P. falciparum, H3K4me3 in hepatic merozoites also had modest peaks associated with the gene start and end sites (Fig. 3A). However, in contrast to what we saw in P. falciparum, we also observed enrichment of these marks in the intergenic regions immediately flanking coding regions.

Fig. 3

H3K4me3 depletion region upstream of genes in P. berghei hepatic merozoites. A Distribution of H3K4me3 and H3K9ac in hepatic merozoites (68–75 hpi) shown as log2(ChIP/input) across all coding sequences (CDS) in the genome, 1.0 kb upstream, and 0.5 kb downstream of the genes. All genes (n = 5,042) were scaled to the same size by computeMatrix. Enrichment of each mark for individual genes is depicted in the heatmap. Genes are sorted by expression level in hepatic merozoites. B Gene expression (log10FPKM) plotted against H3K4me3 or H3K9ac enrichment (log2(ChIP/input)) in the promoter region (1,000 bp upstream of the gene ATG). Line of best fit as determined by linear regression analysis is depicted by red line. The 95% confidence interval of the line of best fit is indicated by dashed red lines. C H3K9ac and H3K4me3 ChIP tracks depicting a region of H3K4me3 depletion upstream of PBANKA_0605700. Gene coding sequence is indicated by the green box. The black rectangle indicates the region of H3K4me3 depletion. Scale bar depicting 500 bp is shown. D Average enrichment of H3K9ac and H3K4me3 in hepatic merozoites depicted as log2(ChIP/input) across 58 regions with H3K4me3 depletion. These regions were scaled to the same size by computeMatrix, and the enrichment 1.0 kb upstream and downstream of the regions is displayed. E Functions associated with genes belonging to the H3K4me3 depletion panel. Functions were assigned based on a literature review and can be found in Additional file 2: Table S19. Genes annotated as ‘conserved with unknown function’ (n = 10) in PlasmoDB are not shown but are included in all analyses. F Heatmap depicting DESeq2 normalized gene expression of H3K4me3 depletion panel genes (n = 64) in 48 hpi hepatic schizonts (48 hpi SZ, n = 2), 54 hpi hepatic schizonts (54 hpi SZ; n = 2), 60 hpi hepatic schizonts (60 hpi SZ; n = 2), hepatic merozoites/DCs (MZ; 69 hpi, n = 2), and erythrocytic early rings (4 hpi, n = 2). G Boxplots depicting the average gene expression of genes with H3K4me3 depletion (Panel, n = 64) and a set of randomly selected genes (Random, n = 147). Gene expression did not differ between time points for the panel of genes with H3K4me3 depletion (Kruskal–Wallis test). H) DNA motif in the H3K4me3 depletion regions

Given the differences in distribution between these histone PTMs in hepatic and erythrocytic stages, we sought to characterize the relationship between these marks and gene expression in hepatic merozoites. A recent study by Caldelari et al. analyzed RNA-seq data from hepatic schizonts (48, 54, and 60 hpi), merosomes/detached cells (DCs) containing hepatic merozoites (69 hpi), and erythrocytic early rings (4 hpi) from P. berghei [49]. Using this data set, we first compared the level of H3K4me3 or H3K9ac to the level of gene expression in hepatic merozoites from the 69 hpi merosomes/DCs. In keeping with observations in P. falciparum, H3K4me3 in the promoter region did not appear to be