Our ChIP-seq results show that unique patterns of H3K36me2&3 occupancy occur at each of the gametocyte stages and that both hPTMs are differentially enriched in the stage II gametocytes (Fig. 1c). The presence of H3K36me2&3 is dynamic with patterns of “low-high-low” hPTM abundance associated with pre-stage II, stage II and post-stage II gametocyte populations, respectively. These stage-specific occupancy profiles were confirmed independently by ChIP-qPCR (Additional file 1: Figure S2). Consistent with the stage-specific enrichment dynamics reported previously for P. falciparum H3K36me2&3, these results validate the experimental strategy used to detect changes in the hPTM landscape during gametocyte development. As highlighted in a central region of chromosome 12 (1250–1290 kb), both H3K36me2&3 were detected at low levels in pre-stage II gametocytes. However, this contrasts sharply with an abundance of both hPTMs in the stage II gametocytes in which the enrichment was preferentially intergenic and upstream of the coding regions of genes. In post-stage II gametocytes, the levels of the hPTMs dissipate, as evidenced by residual H3K36me2&3 occupancy at sites that were highly enriched in the preceding stage II gametocyte samples (Fig. 1c). Given the unique abundance of H3K36me2&3 in the stage II gametocytes and the congruency of this with our previous proteomics data, our downstream analyses focussed on the preferential enrichment of these hPTMs upstream of genes in the stage II gametocytes.
Stage II gametocytes have a unique pattern of H3K36me2&3 enrichment upstream of genes
To determine the genome-wide positioning of H3K36me2&3, we interrogated the hPTM occupancy spanning each of P. falciparum genes for which sequencing data were obtained (5602 genes, Additional file 2: Table S2, Fig. 2a). The occupancy profiles of H3K36me2&3 are distinctly uniform across the intergenic non-coding regions in stage II gametocytes, contrasting with more variable patterns in the pre- and post-stage II gametocytes (Fig. 2a). In stage II gametocytes, both H3K36me2&3 are concentrated 750 bp upstream of gene start sites (SS) with average occupancies (log2ChIP/input) over this region of 0.2 and 0.15 (peaking at 0.24 and 0.2) for H3K36me2&3, respectively, with a contrasting lack of enrichment (log2ChIP/input < 0.2 and < 0.15 for H3K36me2&3, respectively) or complete depletion (log2ChIP/input ≤ 0) in coding regions (Fig. 2a). In the pre- and post-stage II gametocytes, both hPTMs are on average depleted 750 bp upstream of the gene SS (e.g. average H3K36me2 occupancy of -0.14 and − 0.001 over this region in pre- and post-stage II gametocytes, respectively), emphasising the prominent abundance of H3K36me2&3 upstream of genes that is exclusive to the stage II gametocytes.
Next, we quantified the number of genes for which H3K36me2/3 occupancy was detected and stratified this according to location (Fig. 2b). In stage II gametocytes, H3K36me2&3 were detected upstream of 89% (4969) and 83% (4677) of genes occupied (log2ChIP/input > 0) by either of these modifications, respectively (Fig. 2b). Of these, 3942 and 3776 genes had occupancy exclusive to the upstream regions while only a minor proportion (85 and 75) of the genes had H3K36me2&3 exclusively present in the coding regions, respectively (Additional file 1: Figure S3A). This contrasts sharply with the distribution of the hPTMs in pre- and post-stage II gametocytes, where coding region occupancy is evident (Fig. 2b, Additional file 1: Figure S3B). In post-stage II gametocytes, both H3K36me2&3 were distributed relatively equally across the upstream and coding regions of genes (Fig. 2b). Unexpectedly, low level H3K36me2 occupancy (average log2ChIP/input = 0.07) was pervasive in the coding regions of pre-stage II gametocytes (4292 genes, Additional file 1: Figure S3A). Since the pre-stage II gametocyte samples used for ChIP-seq still contained a relatively large subpopulation of asexual parasites, this occupancy likely results from the detection of H3K36me2 in these stages where coding region occupancy has been reported to be linked to gene repression (18, 22).
The wide-spread occupancy patterns of H3K36me2&3 in stage II gametocytes each translate to a large number of genes where the hPTMs are enriched (Fig. 2c). For H3K36me2, 2480 genes were enriched with the hPTM 750 bp upstream of the SS and a further 226 genes were robustly enriched (log2ChIP/input ≥ 0.5). Although also detected in the coding regions of 23% of genes (599) with upstream enrichment, only four of these had H3K36me2 enrichment extending to their coding region (Additional file 1: Figure S3C). Similarly, H3K36me3 was enriched upstream of 2691 genes in stage II gametocytes with 85 of these exhibiting robust enrichment (Fig. 2c). Approximately 25% of these enriched genes also had H3K36me3 associated with the coding regions however, only 15 were enriched for the hPTM in this region (Additional file 1: Figure S3C). Taken together, these results indicate that H3K36me2&3 occupancy occurs in a pattern unique to and characteristic of stage II gametocytes. Furthermore, a concentration of the hPTMs upstream of gene SS, where they would be ideally situated to influence expression by virtue of their proximity to promoters, suggests functional roles for H3K36me2&3, exclusive to stage II gametocytes.
Enrichment of H3K36me2&3 upstream of genes is associated with reduced transcript levels in stage II gametocytes
To explore the functional roles of H3K36me2&3 in stage II gametocytes, we compared the levels of hTPM occupancy upstream of gene SS with expression levels from a publicly available gene expression time course data set (34). Corresponding expression profiles were obtained for days 2–6 of gametocyte development (days on which stage II gametocytes were present) for 96% and 95% (2594 and 2545 genes) of the H3K36me2&3-enriched genes, respectively (Additional file 2: Table S3). The global analysis of the hPTMs upstream of genes and the transcript levels in stage II gametocytes (average log2Cy5/Cy3 on days 4 and 5 of development associated with stage II gametocytes) indicated that the degree of H3K36me2&3 occupancy is not correlated with the expression (Pearson correlations, r2= -0.12 and r2= -0.1 for H3K36me2&3, respectively, Fig. 3). The contrasting depletion of H3K36me2&3 upstream of genes in the other stages, particularly in pre-stage II gametocytes (Additional file 1: Figure S4A), implies that in the stage II gametocytes, these hPTMs are associated with an active repression of genes that are not required for gametocyte development. Indeed, 63% of these H3K36me2/3-enriched and repressed genes (Additional file 2: Table S3) have been shown to be indispensable for asexual proliferation previously (39). This is furthermore supported by the increasing anti-correlation (r2= -0.026 and − 0.02 on day 2 and r2= -0.13 and − 0.14 on day 6 for H3K36me2&3, respectively, Additional file 1: Figure S4A) between hPTM abundance and transcript levels for the enriched genes in stage II gametocytes. Additionally, for genes without H3K36me2/3 enrichment, this relationship was weakened or reversed (corresponding r2 = 0.02 and 0.03 on day 2 and r2= -0.03 and − 0.01 on day 6). The H3K36me2-enriched genes clustered into three groups based on expression profiles (Fig. 3a) with ~ 98% (2548 genes) of the genes showing reduced transcript abundance. A similar pattern was present for H3K36me3 with ~ 98% (2545 genes) of the genes enriched for this hPTM repressed in the stage II gametocytes (Fig. 3b). The H3K36me2&3-enriched repressed gene sets are significantly (P-value < 0.05) associated with shared and unique biological processes (Fig. 3), reflecting the large proportion of genes in the stage II gametocytes that were enriched for both H3K36me2&3 (2012 genes, i.e. 80%), Additional file 1: Figure S4B, Additional file 2: Table S3) and 20% of each set uniquely enriched (695 and 680 genes for H3K36me2&3, respectively). This overlap suggests that one of the H3K36 methylation states is a transient intermediate in the generation of the other (22).
The biological processes associated with genes enriched for H3K36me2&3 include those that are well established to be essential in the proliferative parasite stages (e.g. signal transduction, antigenic variation, gene expression and DNA replication) and include genes coding for Maurer’s cleft proteins (e.g. rex1, PF3D7_0935900 and rex4, PF3D7_0936400) members belonging to the major invasion families (e.g. rhs, ebas, vars, rifs and msps) and components of the cytoadherence complex (e.g. kahsp40, PF3D7_0201800 and kahrp, PF3D7_0202000, Fig. 3, Additional file 2: Table S4) (34, 39, 40). Similar to the role previously ascribed to H3K36me3 as a repressor of invasion gene families when associated with regions upstream of genes in asexual parasites (18, 22), we find the majority of the erythrocyte membrane proteins (emp1s, 75%), stevors (94%) and rifins (70%) are enriched for H3K36me3 in stage II gametocytes, suggesting the control of variant gene expression by H3K36me3 is not limited to asexual parasites but is conserved across multiple life cycle stages. Additionally, most of these genes are also enriched for H3K36me2 as an intermediate, as in the asexual stages (22). Interestingly, within the process of chromatin remodelling, all four P. falciparum histone acetyltransferase (HAT) -encoding genes (14, 41) that generate transcriptionally active chromatin structures (17, 19) were enriched and repressed by H3K36me2&3. Various process that have been documented to be functionally important during sexual commitment and early gametocyte differentiation including chromatin remodelling (P = 0.04), protein translocation (P = 0.01) and phosphatidylcholine (PC) biosynthesis (P = 0.003) (27, 29, 42–46) are significantly associated with the H3K36me2&3-enriched gene sets (Fig. 3). Interestingly, the only processes where H3K36me2&3 enrichment was not associated with a definitive reduction in transcript abundance were actin filament capping, acetyl-CoA biosynthesis and intracellular signalling. Overall, these results indicate that H3K36me2&3 are associated with an active repression of genes that encode protein products involved in proliferation and sexual commitment when they become irrelevant in the developing gametocyte.
In stage II gametocytes, H3K36me2&3 are associated with transcriptional repression of genes that are activated during sexual commitment.
Next, we sought to further investigate the relationship between H3K36me2&3 enrichment in stage II gametocytes and transcriptional repression of commitment- and differentiation-related genes. Of a panel of 1500 sexual development-related genes compiled from previous studies (25, 27, 29, 34, 42–46), we find 64% are enriched for H3K36me2&3 in stage II gametocytes (Additional file 2: Table S5). Importantly, all the enriched genes have dramatically reduced transcript abundance levels during early gametocyte development compared to those without H3K36me2&3 (average of 6-fold and ≤ 12-fold lower on days 3–7, Fig. 4a). The H3K36me2&3-enriched set includes genes encoding the phosophoethanolamine methyltransferase (PMT) and ethanolamine kinase (EK) that are involved in the synthesis of PC from an alternative substrate under lysophosphatidylcholine (lysoPC) limited conditions (42, 46–49). Additionally, H3K36me2&3 enrichment is evident for genes encoding the chromatin modifying enzymes SNF2L, ISWI, the histone deacetylases HDA1 and HDA2 and the NAD+-dependent deacetylase, SIR2A, the latter of which is known to be repressed during early development and thereby blocks proliferative processes (50). Interestingly, we identify H3K36me2&3 enrichment for two histone methyltransferase (HMT) -encoding genes, set3 (PF3D7_0827800) and set9 (PF3D7_0508100). Although the targets of SET9-mediated histone methylation remain unclear, SET3 is known to di- and tri-methylate H3K9 (23, 51); H3K9me3 is well established as a repressive hPTM (16, 25, 27). Moreover, H3K36me2&3 are enriched upstream of the gene encoding Jumonji-C domain-containing protein 2 (JMJC2), a putative P. falciparum HDM containing the Jumonji-C domain (52) that is present in histone lysine demethylase 4 (KDM-4)/JMJC2 HDMs that demethylate di- and tri-methylated H3K9 and H3K36 in other organisms (53–55).
In addition to the above, the ap2-g (PF3D7_1222600) locus was enriched for H3K36me2&3 in stage II gametocytes. AP2-G, a member of the ApiAP2 DNA-binding protein family, has been shown to be a master regulator of sexual commitment in P. falciparum, P. yoelii and P. berghei (28, 30, 39, 56, 57). Here we find some differentiation in the distribution of H3K36me2 vs. -me3 upstream of ap2-g with the former somewhat less defined in this region (Fig. 4b). Furthermore, similar H3K36me2/3 enrichment was present upstream of 67% of the genes (average upstream occupancy of 0.3 and 0.24 respectively, Fig. 4b, Additional file 1: Figure S5A) directly targeted by the AP2-G transcription factor during commitment and subsequent development (29, 58). As before, the H3K36me2&3 enrichment at the AP2-G target genes is associated with their transcriptional repression, contrasting with the expression patterns of those without enrichment (Fig. 4b). In the pre-stage II gametocytes, these AP2-G binding sites are generally depleted of H3K36me2&3 (average upstream occupancy of -0.18 and − 0.3, respectively, Additional file 1: Figure S5A), including the ap2-g locus itself with < 1% of the target genes associated with H3K36me2/3 enrichment in this stage (Additional file 1: Figure S5B). This shows that H3K36me2&3 are involved in the active repression of genes once the protein products are no longer required following the establishment of gametocytogenesis. Similarly, in the post-stage II gametocytes, these AP2-G binding sites are generally associated with an absence of H3K36me2&3 (average occupancy of 4.2e− 4 and 1.8e− 3, respectively, Additional file 1: Figure S5A) with < 7% of these genes exhibiting residual H3K36me2/3 enrichment from the preceding stage II gametocytes. This is associated with a gradual increase in transcript abundance, indicating that these genes are no longer actively repressed once depleted of H3K36me2&3 enrichment in the post-stage II gametocytes (Additional file 1: Figure S5C). Although the targets of AP2-G vary between gametocytes arising from the same-cycle conversion (SCC) and next-cycle conversion (NCC) pathways (29), we find no considerable differences in the proportion of SCC and NCC AP2-G target genes that are enriched for H3K36me2&3 in stage II gametocytes (Additional file 1: Figure S5D), suggesting that the H3K36me2&3-associated transcriptional repression occurs post-commitment irrespective of the route taken during sexual differentiation.
We also assessed the H3K36me2/3 enrichment in our gametocyte populations at sites occupied by a second ApiAP2 transcription factor, AP2-G2 (PF3D7_1408200), during intermediate gametocyte development (59). AP2-G2 is not essential for asexual parasite proliferation or sexual commitment however, the genetic disruption of ap2-g2 stalls normal gametocyte development beyond stage III (39, 59). In pre-stage II gametocytes, 3% of the upstream AP2-G2 binding sites were enriched with H3K36me2/3. This increases drastically in stage II gametocytes where 92% of these sites are now enriched with H3K36me2/3 (Fig. 4c). This suggests an interaction between these two regulatory mechanisms and supports previous descriptions of AP2-G2 as a transcriptional repressor (59). The number of H3K36me2&3-enriched sites is somewhat reduced in the post-stage II gametocytes (54%) as expected given the residual hPTM occupancy patterns in these later gametocyte stages that are reminiscent of H3K36me2/3 enrichment in the stage II gametocytes. For AP2-G2 binding sites present in coding regions, 0.3% have H3K36me2/3 enrichment in the stage II gametocytes (Fig. 4c) and in the pre- and post-stage II gametocytes, this this increases slightly to 7.5% and 11%, respectively. This implies the potential interaction of H3K36me2&3 with AP2-G2 occurs predominantly in regions upstream of genes with the enrichment of these hPTMs at these locations suggesting a stage-specific role for this interaction.
Lastly, we compared H3K36me2&3 with other epigenetic mechanisms that function in P. falciparum sexual commitment and development or have been implicated in gene repression and chromatin remodelling (16, 25, 27, 60, 61). The interaction between H3K9me3 and histone protein 1 (HP1) mediates global heterochromatin formation during gametocyte development in addition to its role in regulating commitment (27, 33). We find that of the 89 genes preferentially occupied by HP1 in stage II/III gametocytes relative to schizonts (27), 88% are enriched for H3K36me2&3 in stage II gametocytes, resulting in diminished transcript abundance (average of 1.5-fold and ≥ 2.6-fold lower on day 3–7 compared to those without H3K36me2&3, Fig. 4d). This implies some degree of co-operativity exists between the repressive actions of H3K36me2&3 and H3K9me3/HP1 binding. Interestingly, each of the 15 genes with reduced HP1 occupancy (27) were also enriched for H3K36me2&3 in stage II gametocytes (Fig. 4d) with the presence of the two hPTMs overlapping with residual HP1 occupancy upstream of gene SS (Additional file 1: Figure S5E) These genes, including ap2-g and the early gametocyte markers PF3D7_1476500, PF3D7_1476600, PFD7_1477300 (Pfg14_744), PF3D7_1477400, PF3D7_1477700 (Pfg14_748) and PF3D7_1478000 (gexp17) (43, 44, 62, 63), exhibit peak transcript abundance during commitment and early differentiation (days 1 and 3) followed by a steady decline in expression levels associated with the enrichment of these genes with H3K36me2&3 in the stage II gametocytes (Fig. 4d). This finding, together with the documented absence of H3K9me3 at these sites in late (stage IV/V) gametocytes (33), suggests that although the identities of the H3K36me2&3 reader protein/s remain unclear, these hPTMs are also involved in the repression of transcription independent of the H3K9me3/HP1 mechanism during gametocyte development. Interestingly, for genes that do gain H3K9me3 occupancy in late-stage gametocytes, 82% were enriched for H3K36me2&3 in stage II gametocytes (Fig. 4e). Once again, the H3K36me2&3-enriched gene set had substantially lower transcript levels compared to those without such enrichment, supporting the idea of co-operative functioning of H3K9 and H3K36 tri-methylation in P. falciparum gametocytes. Together, these results indicate that H3K36me2&3 are involved in generating a post-commitment transcriptional shift to allow differentiated gametocytes to continue through development.
Inhibition of H3K36 demethylation by JIB-04 is associated with altered patterns of transcription in P. falciparum gametocytes
Lastly, we investigated the potential enzymes that generate the unique stage II gametocyte-specific enrichment of H3K36me2&3 described here and previously (11). In asexual parasites, SET2 (PF3D7_1322100) methylates H3K36 (52, 64). The transcript levels of this HMT increase on day 1 of gametocytogenesis and subsequently decline thereafter (Fig. 5a), suggesting that the H3K36-specific methyltransferase activity of SET2 is extended to the gametocyte stages. The transcript levels of the three Jumonji KDMs of P. falciparum (jmjc1 (PF3D7_0809900), jmjc2 (PF3D7_0602800) and jmj3 (PF3D7_1122200)) spanning stage II gametocyte development strongly suggests that at least one of these enzymes is likely involved in H3K36 demethylation (Fig. 5a). To investigate the functional relevance of these enzymes for H3K36me2&3, we chemically interrogated HDM activity and determined the effects of this inhibition on H3K36 methylation levels on days associated with stage II gametocyte development. Jumonji HDM inhibitors with activity against gametocytes (38, 65, 66). JIB-04 (pan-selective) and ML324 (targeting KDM4, PfJMJ3) resulted in the hypermethylation of both H3K36me2&3 relative to the untreated controls (36- to 89-fold, Fig. 5b, Figure S6A), similar to hypermethylation of H3K4me3 and H3K9me3 induced by JIB-04 and ML324 in asexual parasites and gametocytes, respectively (38, 66). This hypermethylation was not observed for GSK-J4 (KDM6, H3K27me3 selective, (67) or PCPA-2 (targeting LSD1). These results confirm that at least one of the P. falciparum Jumonji HDMs demethylate H3K36me2&3 after they peak in abundance in stage II gametocytes (11). Furthermore, mechanistic studies in asexual parasites have demonstrated that the inhibition of demethylase activity by JIB-04 arises from the direct targeting of PfJMJ3 (38). In the asexual parasites, the inhibition of PfJMJ3 by JIB-04 results in transcriptional de-regulation (38). Since this inhibitor is significantly more potent towards the sexual stages (38, 66), we were particularly interested in the transcriptional effects of JIB-04 in gametocytes since this inhibitor is significantly more active against the sexual stages. As such, we performed genome-wide transcriptional profiling of JIB-04 activity on day 3 and 4 gametocytes each following 24 h treatment with JIB-04 (Additional file 3: Table S6). JIB-04-treated gametocytes displayed differential expression of ~ 13% of the genome (711 and 696 genes with decreased and increased transcript abundance, respectively, Fig. 5c, Additional file 1: Figure S6B), similar to the restricted effect induced by this inhibitor in asexual parasites (38). Differentially expressed genes with decreased transcript abundance are associated with various processes including chromatin organisation (3%), gene expression (12%) and transport (6%), each of which have confirmed roles in gametocyte differentiation and development (27, 34, 43, 45, 46, 63, 68–70) (Fig. 5c). Of these genes, 63% were enriched for H3K36me2&3 (93% with occupancy) in stage II gametocytes including the invasion-related gene families (Fig. 5c) whose post-commitment transcriptional repression we show to be associated with these hPTMs. A large proportion (93%) of the differentially expressed genes with increased abundance also had some degree of H3K36me2/3 occupancy with 56% enriched for the hPTMs in the stage II gametocytes. Genes functionally associated with cellular proliferation that were enriched for H3K36me2&3 in stage II gametocytes and repressed as a result, also had increased transcript abundance in response to JIB-04 treatment (Fig. 5c). Interestingly, we also find an increase in the transcript levels of jmjc2 and the HAT-encoding genes, gcn5 (PF3D7_0823300) and myst (PF3D7_1118600) (Fig. 5c) (14, 41), all of which may reflect an attempt to counteract the abnormal H3K36me2&3-associated repression induced by JIB-04. Combined with the depletion of H3K36me2&3 beyond stage II gametocyte development, these results suggest that aberrant H3K36 hypermethylation poses a barrier to further development and highlight the importance of histone methylation for transcriptional reprogramming during gametocyte differentiation and development. Shared disruption of certain biological processes, including chromatin organisation, cell motility and kinase/phosphatase activity (38) was present (Fig. 5d). The JIB-04 treated gametocytes do not share common differentially expressed genes with the G9a methyltransferase specific inhibitors, BIX-01294 (71) but some overlap is present with HDAC inhibition by Trichostatin A (72), in line with other studies which showed similarities in the gene expression signatures of JIB-04 and TSA in cancer cell lines (73).