Genome-wide localization of histone variants in Toxoplasma gondii implicates variant exchange in stage-specific gene expression

Toxoplasma gondii is a protozoan parasite that differentiates from acute tachyzoite stages to latent bradyzoite forms in response to environmental cues that modify the epigenome. We studied the distribution of the histone variants CenH3, H3.3, H2A.X, H2A.Z and H2B.Z, by genome-wide chromatin immunoprecipitation to understand the role of variant histones in developmental transitions of T. gondii parasites. H3.3 and H2A.X were detected in telomere and telomere associated sequences, whereas H3.3, H2A.X and CenH3 were enriched in centromeres. Histones H2A.Z and H2B.Z colocalize with the transcriptional activation mark H3K4me3 in promoter regions surrounding the nucleosome-free region upstream of the transcription start site. The H2B.Z/H2A.Z histone pair also localizes to the gene bodies of genes that are silent but poised for activation, including bradyzoite stage-specific genes. The majority of H2A.X and H2A.Z/H2B.Z loci do not overlap, consistent with variant histones demarcating specific functional regions of chromatin. The extent of enrichment of H2A.Z/H2B.Z (and H3.3 and H2A.X) within the entire gene (5’UTR and gene body) reflects the timing of gene expression during the cell cycle, suggesting that dynamic turnover of H2B.Z/H2A.Z occurs during the tachyzoite cell cycle. Thus, the distribution of the variant histone H2A.Z/H2B.Z dimer defines active and developmentally silenced regions of the T. gondii epigenome including genes that are poised for expression. Histone variants mark functional regions of parasite genomes with the dynamic placement of the H2A.Z/H2B.Z dimer implicated as an evolutionarily conserved regulator of parasite and eukaryotic differentiation.

transcription and DNA repair. Eukaryotic cells have evolved three main strategies to regulate nucleosomes. First, ATP-dependent chromatin remodeling complexes use the energy from ATP hydrolysis to remodel chromatin by replacing or sliding nucleosomes, thus affecting their interaction. Second, histones are extensively modified with posttranslational modifications (PTMs). These PTMs regulate DNA-related processes in a positive or negative way, influencing the interaction between the histones and DNA or recruiting protein factors involved in gene expression or gene silencing. Finally, cells replace canonical histones with histone variants [4,5] . Unlike canonical histones, whose synthesis is coupled to DNA replication during S-phase, histone variants are typically constitutively expressed through the cell cycle [6] . The physical properties of histone variants can change chromatin structure, or histone variants may localize to specific regions of the genome to confer specialized properties to chromatin.
Toxoplasma gondii chromatin structure has been implicated as an important regulator of gene expression. T. gondii is an opportunistic parasite that infects between 10 and 90% of the population depending on the country [7]. Although T. gondii infection is asymptomatic in healthy people, untreated clinical toxoplasmosis is lethal in immunocompromised individuals such as transplant or HIV patients. In addition, toxoplasmosis is an important cause of congenital defects or fetal demise when women become infected during pregnancy. In humans and other intermediate hosts, T. gondii tachyzoite stages are associated with acute disease, while bradyzoites within tissue cysts persist in a chronic phase. The ability of the parasite to persist as a latent bradyzoite form is critical for disease pathogenesis and the developmental transition from tachyzoite to bradyzoite is likely driven by epigenetic processes involving histone PTM and histone localization [8].
T. gondii expresses four canonical histones, represented by single copy genes, except for H2B that is encoded by two separate genes on different chromosomes [9,10]. H2Ba is expressed mainly in tachyzoite forms while H2Bb expression is restricted to the sexual stages [11,12] . H4 does not have a variant, but interacts with H3 and its variants H3.3 and centromeric CenH3 [13]. H2A has two variants: H2A.X, which shares 91% identity with canonical H2A, and H2A.Z, which is considered the most conserved histone variant across species [14] . To date, no histone H1 has been identified in Apicomplexa.
H3.3 differs from H3 by just four amino acids in T. gondii (RY instead of KF at position 54-55; QAIL in H3. 3 and SAVL in H3 in positions 88-91) [15]. Although there are only minor sequence differences between H3.3 and H3, these few differences can result in different roles in gene activation [16]. CenH3 is a conserved variant necessary for kinetochore formation in several organisms, and TgCenH3 is located at centromeric regions of T. gondii [13]. T. gondii centromeres are arranged in a single apical cluster, represented by a dot that is attached to the centrocone during all cell cycle by immunolocalization [13].
T. gondii H2A.X has the conserved SQ(E/D) φ motif important for DNA repair (φ represents a hydrophobic amino acid and S is a serine that is phosphorylated in response to DNA damage). Serine 132 of H2A.X is phosphorylated in tachyzoite forms of T. gondii and plays a role in DNA repair [14,17]. TgH2A.Z shares 82% identity with human H2A.Z, though TgH2A.Z has an additional 20 amino acids in the N-terminal region. By Chromatin Immunoprecipitation followed by qPCR (ChIP-qPCR), TgH2A.Z was detected at the promoter regions of active genes, a conserved characteristic of H2A.Z in most organisms, whereas TgH2A.X was associated with silent genes [14].
H2B.Z (previously known as H2BV), is a variant unique to Apicomplexa [12,18,19] that has been proposed to be the major histone H2B in T. gondii tachyzoites. H2B.Z interacts with H2A.Z but not H2A.X [14]. Mass spectrometry analysis of T. gondii tachyzoites revealed that both H2B.Z and H2A.Z are hyperacetylated in their N-terminal domains, a classic hallmark of active chromatin [20], and both are also substochiometrically ubiquitinated [21]. By contrast, few or no acetylation marks have been detected on canonical H2A, H2B or H2A.X [20].
In this study, we surveyed the genome-wide location of histone variants across the T. gondii genome, using Chromatin Immunoprecipitation. Our findings implicate histone variant exchange in regulation of stage-specific transcription and expand our understanding of histone variants in eukaryotes.

Localization of histone variants is not uniform within nuclei
We confirmed the location of histones variants in the nuclear compartment by labeling intracellular RH strain tachyzoites with monospecific polyclonal antibodies recognizing T. gondii histone variants or endogenous epitope-tagged proteins, and imaging by deconvolution fluorescence microscopy (Fig. 1). Comparison of localization patterns with the staining intensity of DAPI allowed us to examine the presence of the histone variants in dense and less dense chromatin. Regions of "heterochromatin" stain more intensely with DAPI than the "euchromatic" domains. None of the histone variants localized exclusively to heterochromatic or euchromatic chromatin. H3.3 showed punctuate distribution, generally opposite to DAPI staining, while H2B.Z and H2A.Z exhibited a limited degree of co-localization with dense chromatin (Fig. 1). H2A.X localized most intensely with regions of compact chromatin, but not exclusively to heterochromatin. None of the histones colocalized with DAPI within the apicoplast, a remnant organelle of endosymbiosis present in many Apicomplexan parasites.

Native ChIP-chip and ChIP-seq in T. gondii
ChIP followed by tiling microarray hybridization (ChIPchip) or high throughput sequencing (ChIP-seq) are widely used for studying genome-wide protein-DNA interactions. In previous work, we applied ChIP-chip to T. gondii and showed that the distributions of epigenomic modifications on chromatin are conserved in protozoa [22]. We standardized native chromatin immunoprecipitation (N-ChIP) for T. gondii using native DNA without crosslinking digested with Micrococcal nuclease (MNase) followed by next generation sequencing [23].
The digested fragments of around 150-200 bp represent nucleosomal DNA, as verified by nucleosome mapping of input DNA. As controls, we performed N-ChIP with antibodies specific for H3K4me3 and CenH3. Using N-ChIP-seq, we observed similar locations to those in our previous studies using cross-linked DNA [13,22].
To study the genome-wide localization of histone variants in T. gondii, both ChIP-seq and ChiP-chip were performed using antibodies described above and epitope-tagged lines in the Type I RH strain background. A summary of experiments and platforms used is shown in additional file 1, along with a summary of how the data was processed and analyzed. Histone variants display distinct localization patterns that differ over genes with varying expression profiles such as LDH1 (tachyzoite) and LDH2 (bradyzoite), ENO1 (bradyzoite) and ENO2 (tachyzoite) or SAG1 (tachyzoite) and BAG1 (bradyzoite) ( Fig. 2A). ChIP-seq showed enrichment over the same genes as ChIP-chip ( Fig. 2A and additional file 2A) and positive correlation between H2A.Z, H3K4me3 and H2B.Z ChIP-chip and ChIP-seq replicates (Pearson's r > 0.75) (additional file 2). Correlations between ChIP-seq replicates are shown in additional file 3.

Localization of histones and histone-variants across the genome
The correlation between histone variant locations was examined by correlating the read densities of each histone variant ChIP-seq dataset (Fig. 2B). H2A.Z and H2B.Z signals are strongly correlated, consistent with their reported association as a dimer in nucleosomes [14]. H3K4me3 also correlated with H2A.Z and H2B.Z. Histones H2A.X and H3.3 distributions are positively correlated, and weakly negatively correlated with H2A.Z/H2B.Z (Fig. 2B). This supports the hypothesis that in T. gondii, H2A.X and H2A.Z are found in different nucleosomes [12,14]. The same pattern was observed at promoters, 5'UTRs, and gene bodies (Fig. 2B). Normalized read densities at promoters, 5'UTRs, gene bodies and 3'UTRs are shown in Fig. 2C. At 3'UTRs, little signal was detected for any histone variant (Fig. 2C).
In tachyzoites, a majority of genes were marked by H3K4me3 (~ 75%), H2A.Z (~ 80%) and H2B.Z (~ 90%) (additional file 4A). H3K4me3, H2A.Z and H2B.Z peaks localize predominantly to promoter regions; H3.3 and H2A.X mark ~ 75% and ~ 45% of genes respectively and localize to regions upstream or downstream or overlapping of mapped TSS (additional file 4B, C, D). Together, these findings imply that nucleosomes at promoter regions and gene bodies have distinct compositions. Promoter regions are marked by nucleosomes containing H2A.Z, H2B.Z and presumably canonical H3, while gene bodies are occupied by nucleosomes containing H2A.X and H3.3 and possibly canonical H2Ba [11]. HFF cells were grown on glass slides and infected with tachyzoites for 24 h. Slides were stained with DAPI and antibodies against anti-H2AZ, anti-myc (for H2AX-myc and H2BZ-myc) or anti-HA (H3.3). Images were collected along a Z-series and processed for three-dimensional deconvolution using similar signal thresholds for each of the antibodies. The small extranuclear DAPI staining in some cells labels the genome of the apicoplast. The merged figures correspond to the overlay of DAPI (blue) and histone variants (green). The scale bars correspond to 5 μm and LDH2. B. The overall similarity in ChIP-seq profiles was compared by calculating Pearson correlation coefficients genome-wide or for different genomic features including promoters, gene bodies and 3′ UTRs (as defined by www. toxodb. org annotations). Correlations were visualized in a clustered heatmap (generated using deeptools http:// deept ools. readt hedocs. io); color indicates correlation coefficient as shown in color key. Peaks for each variant represent the intersection of peaks from each biological replicate. C. Average ChIP-seq profiles of each histone variant on different genomic features. Start and end denote boundaries of each feature, with 500 bp upstream and downstream flanking sequence either side. ChIP-seq read pileups were normalized by scaling to library size Since H3K4me3 colocalizes with H2A.Z/H2B.Z, these data imply exchange of H2Ba/H2A or H2Ba/ H2A.X dimers, with double variant dimer composed of H2A.Z/H2B.Z [14,24,25] at promoter regions of active genes.

Centromeres and telomeres have unique histone compositions
Centromeres and telomeres have distinct functions and are often associated with altered histone composition. We examined histone variant localization at these specialized regions. CenH3 localizes almost exclusively to centromeric regions and is accompanied by H3.3, while H2A.X, H2A.Z and H2B.Z are absent from centromeres (Fig. 3A). CenH3 localized to 13 different loci and was absent from Chromosome VIIb, consistent with the revised genome of T. gondii containing 13 chromosomes [26]. Further work is required to determine whether CenH3 and H3.3 are part of the same nucleosome within a double variant heterodimer of CenH3 and H3.3, or if T. gondii centromeres consist of interspersed CenH3 and H3.3-containing nucleosomes that may in turn be enriched at different phases of the cell cycle [27].
At telomeres and telomeric associated sequence (TAS), CenH3 is absent and H2A.X and H3.3 are enriched (Fig. 3B), suggesting that they may play a role in telomeric and TAS maintenance [11]. As in other species, T. gondii telomeres and TAS are notable for repetitive DNA sequences [28] making accurate genome annotation and assembly difficult. Nevertheless, H3.3 and H2A.X showed high levels of signal at the ends of at least 10 of the 14 chromosomes (data not shown), suggesting a role in protecting the extremities of telomeres. T. gondii telomeres are also associated with the H3K9me3 repressive mark in ChIP-Seq [29] and the HP1 orthologue Chromo 1 [30]. H2A.X, along with H3.3, was also found in lower levels spread through the chromosomes, usually opposite to H2A.Z ( Fig. 2A).

H2A.Z and H2B.Z mark promoters of active genes
H2A.Z has a conserved function in transcriptional activation [28]. In T. gondii, 98% of H3K4me3 peaks colocalize with H2A.Z and H2B.Z in or around promoter regions (Figs. 2 and 4 A, B). Even in the asynchronous tachyzoite cultures we sampled, strong nucleosome positioning was evident at the 5′ region of active genes. Several nucleotide motifs were enriched in these regions ( Fig. 4 C), suggesting that H2A.Z and H2B.Z are positioned with H3K4me3 and H4ac in promoter regions [22] that recruit macromolecular complexes consisting of transcription factors and chromatin modifying enzymes.
Distinct profiles of histone variant coverage were observed on mRNA compared to rRNA or tRNA (additional file 5A). rRNA genes, which are highly transcribed in tachyzoites, were relatively depleted of histones, Fig. 3 Composite profiles of histone variant distribution at centromeres and telomeres. Profiles were generated by averaging the ChIP-seq read density over annotated genomic loci for the histone variants as indicated by the colored lines. ChIP-seq profiles were normalized by scaling to 1x library size. Start and end denote boundaries of each feature, with 500 bp upstream and downstream flanking sequence either side. A. Centromeric regions. B. Telomeric regions whereas tRNA genes were flanked by H3K4me3. The presence of H3K4me3, H2A.Z and H2B.Z at promoters and 5'UTR was associated with gene expression (additional file 5B), although there did not seem to be significant differences in ChIP-seq signal when genes with low, medium or high levels of expression (determined from RNA-seq data [31]) were compared. H3K4me3, H2A.Z and H2B.Z were not enriched in the 5′ regions of genes whose expression was in the lowest expression quintile. Normalized read densities were similar for all variants at G1-regulated or SM-regulated genes (additional file 5C).

H2A.Z and H2B.Z mark gene bodies of silent genes
Because H2A.Z also has a role in transcriptional repression [32], we examined its distribution on silent genes. Surprisingly, manual inspection of ChIP profiles suggested that in addition to localizing to active promoters, H2A.Z and H2B.Z are also present at low levels on gene bodies of silent genes, between TSS and transcription termination sites (TTS) (Fig. 5). We examined the genome-wide distribution of H2A.Z and H2B.Z on silent genes and stage-specific genes (Fig. 5A, additional file 6). Enolase 1, which is expressed in bradyzoites, and enolase 2, which is expressed in tachyzoites, Fig. 4 H3K4me3, H2A.Z and H2B.Z overlap at promoters. A. Genome browser view (visualized using IGV) of H3K4me3, H2A.Z and H2B.Z ChIP-seq read density profiles at two genes (IMC1 and IMC4) with opposing orientations (bottom track). RNA-seq from T. gondii RH strain tachyzoites is shown in grey (from [31]). Peaks called using the peak-calling algorithm MACS2 are shown under read density profiles. B. Location of combined H3K4me3, H2A.Z and H2B.Z peaks in relation to promoter regions (defined as 1000 bp upstream of 5'UTR start, defined from RNA-seq data on tachyzoites [31]). C. Sequence motifs identified by HOMER that are significantly enriched in genomic regions where H3K4me3, H2A.Z and H2B.Z overlap (See figure on next page.) Fig. 5 H2A.Z and H2B.Z mark active and silent genes. A. Average histone variant profiles at genes that are active or silent in tachyzoites stages based on RNA-seq data [31]. TSS (transcriptional start site) and TTS (transcriptional termination site) denote boundaries of genes, with 2 kb upstream and downstream flanking sequence. Gene boundaries as defined by genome annotation were used if experimental TSS and TTS were not available. ChIP-seq signals were normalized by scaling to library size. B. Genome browser view of H2A.Z, H2B.Z and H3K4me3 at an active (left panel) and a silent (right panel) gene. Peaks called using MACS2 are shown under read density profiles. C. Statistical approach for investigating link between histone variant gene coverage and stage-specific gene expression. Genes marked by histone variant ChIP peaks were split into groups depending on what percentage of the gene (starting at 5'UTR, to 3'UTR) is covered by a peak. Each group of genes was compared to predefined gene sets containing genes that are upregulated in different parasite stages [12,29] and a measure of enrichment (p-value) was calculated using a hypergeometric distribution test. D. Plots showing -log2(p-value) of enrichment for genes with different amounts of H2A.Z, H2B.Z or H3K4me3 coverage compared to tachyzoite or bradyzoite gene sets. Genes with different levels of H2A.X or H3.3 coverage in tachyzoite or bradyzoites gene sets and significance of enrichment (−log 2 (p-value)) calculated are shown for comparison. p-value of < 0.05 in the hypergeometric test was considered significant and is indicated by the dotted grey line map to a tandem array on chromosome VIII ( Fig. 2A). Both H2A.Z and H2B.Z are observed on promoters of active ENO2 gene but also on the top of the CDS of silent ENO1 gene ( Fig. 2A). Similarly, H2B.Z and H2A.Z are distributed over the CDS in silent LDH2 gene and at the promoter of the active LDH1 gene (Fig. 5 B and additional file 2).
To investigate whether this trend occurs genome-wide, we extracted lists of genes marked by histone variant peaks and performed gene set enrichment analysis using previously defined gene sets [21,33]. This analysis was described previously to analyze T. gondii data and is commonly used for pathway analysis in numerous organisms. The T. gondii gene sets that were used were comprised of genes that are upregulated in different developmental stages, cell cycle time points or subcellular localization. A hypergeometric distribution test (also known as Fisher's Exact Test) was used to calculate the statistical likelihood of enrichment or overrepresentation of genes in a particular gene set compared to a randomized gene set of the same size, yielding a p-value of enrichment which we term 'enrichment score' . Overrepresented gene sets were defined as those with a p-value of < 0.05 in the hypergeometric test [21,34,35]. To test whether genes with different amounts of H2A.Z or H2B.Z coverage were overrepresented in tachyzoite or bradyzoite genesets, we split tachyzoite and bradyzoite genes into groups based on the percentage of the gene that was covered by an H2A.Z or H2B.Z peak (e.g. 0-20%, 20-40% of the gene length beginning at the promoter) coverage ( Fig. 5C and D). Enrichment scores were plotted for histone variant peaks across genes expressed in both developmental stages (Fig. 5 D and additional file 6). While genes with lower H2A.Z and H2B.Z gene coverage (e.g. 0-20%) are statistically enriched in genes upregulated in tachyzoites, bradyzoite genes are overrepresented in high coverage gene sets (e.g. 80-100% coverage).
As a comparison, the same analysis was performed on H3K4me3, H2A.X and H3.3 peaks. As expected, since H3K4me3 is exclusively found at promoters in every organism to date, tachyzoites genes are overrepresented among genes with 0-40% H3K4me3 coverage, and genes with 40-60% H3K4me3 coverage are enriched in tachyzoite and not bradyzoites genes. Genes marked with 80-100% H2A.X and H3.3 are statistically enriched for tachyzoite genes, but not bradyzoite genes (Fig. 5D). Together, these findings suggest that the genome positions of H2A.Z and H2B.Z are linked to control of stageregulated genes.

Cell cycle control and histone variants
Canonical histones play a critical role in cell cycle control. Their synthesis is coupled to DNA replication during S-phase and is essential for chromatin assembly. In contrast, histone variants are synthesized throughout the cell cycle and whether they affect cell cycle progression is unclear [30]. To investigate whether histone variants play a role in regulating cell cycle dependent transcripts in T. gondii, we used the enrichment analysis method described above to compare genes with different amounts of histone variant coverage to cell cycle gene sets defined previously [33] . During the 8 h T. gondii cell cycle, there are two waves of transcription corresponding to G1 and S/M phases [36]. G1-regulated genes mainly have basal eukaryotic functions such as metabolism or biosynthetic processes, while SM-regulated genes are important for parasite maturation, division and biogenesis of parasitespecific organelles.
Genes with low amounts of H2A.Z and H2B.Z coverage (0-20% and 20-40%) are enriched in G1 genes, while genes expressed later in the cell cycle have higher amounts of coverage (Fig. 6). Interestingly, genes with 0-20% coverage are enriched in genes upregulated at a slightly earlier time point (4-6 h) than genes with 20-40% coverage (5-7 h), suggesting that coverage correlates with the timing of expression. Recent work using single cell sequencing [37] demonstrated that T. gondii G1 is divided into G1a followed by G1b, substages of G1 that are characterized by different sets of cell cycle expressed genes. Differentiated bradyzoites are predominantly in G1 [38] and single cell gene expression patterns of bradyzoites are those of G1b cells [37].
Turning to S/M-regulated genes, genes with both low and high levels of H2A.Z and H2B.Z coverage are overrepresented in genes upregulated in mid-S/M (3-5 h) phase. This is the time period during which parasites begin to express bradyzoite markers that differentiate to bradyzoites [38], and many parasite-specific proteins are upregulated at this point. Of note, a similar enrichment pattern of both high and low levels of H2A.Z and H2B.Z coverage were seen in merozoite genes, which are expressed in the feline intestine stages (additional file 6).
A similar analysis showed that genes with high coverage (80-100%) of H2A.X and H3.3 are enriched in G1 genes upregulated between 5 and 7 h of the cell cycle (Fig. 6). Genes with high H3.3 coverage were also enriched for mid-S/M-regulated genes, as well as genes with 20-40% H2AX. G1-regulated genes had higher H2AX and H3.3 coverage, consistent with the observation by manual inspection that these variants are present on gene bodies. Comparison of the genome-wide distributions of these variants suggests that the composition of histones within nucleosomes is dynamic and varies throughout the cell cycle as different sets of genes are transcribed.

Discussion
T. gondii, an opportunistic parasite with worldwide distribution, has several developmental transitions that are regulated by epigenetic processes. Here, we provide a genome-wide overview of localization of T. gondii histone variants, which are linked to chromatin structure and regulation of gene expression. T. gondii has five histone variants: H3.3, H2A.X, H2A.Z, H2B.Z and CenH3. CenH3 has a conserved function at centromeric regions in eukaryotes including T. gondii [13] and we confirmed the localization of CenH3 to centromeres by ChIPseq. To analyze the location of other variants across the genome, we used high-throughput chromatin immunoprecipitation followed by either microarray hybridization [22] or sequencing.
In tachyzoite stages, H3K4me3 is absent on promoters of bradyzoite genes and other silent genes. In the Apicomplexan Plasmodium falciparum, H3K4me3 and H3K9ac were shown to be present on promoters of active genes in schizont stage but on both active and silent genes in ring stages [39]. The difference in results may reflect divergent chromatin biology or may be a result of experimental conditions. Plasmodium varies from T. gondii in having a longer 48-h asexual cell cycle (vs 6-8 h for T. gondii), therefore chromatin state and gene expression can be examined simultaneously in synchronized cells, an approach that is technically challenging given the short life-cycle of T. gondii and the large amounts of biological material required for N-ChIP.
At T. gondii promoters, H3K4me3 colocalizes with H2A.Z and H2B.Z, whereas, when present on gene bodies, H2A.Z and H2B.Z are not accompanied by H3K4me3, suggesting an additional role for these histone variants, perhaps in gene silencing. An open chromatin state, marked by H2A.Z/H2B.Z, H4ac and H3K4me3, is required for expression of genes in asexual stages of both Plasmodium spp. and T. gondii [8,40]. H2A.Z and H2B.Z form a dimer [14] and colocalize at promoter regions of active genes.
H2A.Z is the most evolutionarily conserved histone variant, while H2B.Z is specific to Apicomplexan parasites. Both these variants differ from canonical histones within the N-terminal domain. H2A.Z has an additional ~ 20 aa region that is extensively modified by PTM [20]. In P. falciparum, the H2B.Z/H2A.Z dimer is enriched at TSS and 5′-UTR and is associated with well-known activation marks H3K4me3 and H3K9ac as well as H3K18ac and H3K27ac [41]. The H2A.Z and H2B.Z dimer is excluded from heterochromatin regions such as silent var genes encoding variant antigens, but is present at the promoter of the single active var copy [25].
Bradyzoite genes are silent in tachyzoites and rarely expressed in the RH strain. H2A.Z/H2B.Z coverage of silent bradyzoite genes extends from the promoter along the gene body, but the levels detected are lower than the prominent peaks seen at promoters of active tachyzoite genes. We also observed enrichment of H2A.Z and H2B.Z on gene bodies on a subset of merozoite genes (additional file 6). We hypothesize that genes with H2A.Z/H2B.Z coverage extending along gene bodies are poised for transcription, as occurs in yeast [42,43] and embryonic stem cells [44]. In mammalian embryonic stem cells, poised but silent promoters are marked by H3K4me3 and H3K27me3 as well as H2A.Z [44]. In T. gondii, these broad H2A.Z/H2B.Z peaks may play a role in state-specific gene regulation during the tachyzoitebradyzoite transition.
It is likely that differences in PTM of H2A.Z/H2B.Z or other histones within the same nucleosome are responsible for regulating transcription positively or negatively. In cancer cells, acetylated H2A.Z (H2A.Zac) is enriched at promoters of poised and active genes, but deacetylated H2A.Z [45] is spread throughout the entire promoter when the gene is silent. In T. gondii, H2A.Z and H2B.Z are extensively acetylated in the N-terminal domain [20], a classic mark of transcriptional activation. In contrast, relatively few acetylation sites have been detected on canonical H2A, H2B or H2A.X in tachyzoites [20]. Taken together with the current study, it is likely that acetylated H2A.Z and H2B.Z localize to promoter regions and play a role in transcriptional regulation. Future studies in T. gondii examining the distribution of H2A.Zac/H2B.Zac and deacetylated forms are needed, to confirm whether the acetylation state of H2A.Z/H2B.Z is necessary for transcriptional activation [46]. Several other PTM have been detected on H2A.Z [47,48]. While hyperacetylated H2A.Z is associated with transcriptional activation, ubiquitinated or methylated H2A.Z is associated with gene repression [32]. In T. gondii, H2A.Z and H2B.Z each have seven ubiquitination sites in both N-and C-terminal domains with unidentified functions [21]. Ubiquitination was detected on TgH2AK119Ub and TgH2BK120Ub [21], residues that are important for recruitment of repressive complexes [49] and transcriptional activation [50], respectively. H2A K119 ubiquitination is mediated by E3 ligase activity of the polycomb repressive complex PRC1 [51]. Acetylation and ubiquitination may regulate the differential localization of H2A.Z and H2B.Z on active and silent genes in T. gondii. Acetylated H4K31 is enriched in promoters of active genes, whereas methylated H4K31 is found in silent regions and the body of inactive genes [29]. Collectively, these data suggest that coordinated deposition of other PTM on histones plays a critical role in T. gondii gene expression and developmental transitions.
In mammalian cells, deacetylated H2A.Z localizes to centromeres and co-localizes with the classic centromeric marks H3K4me2, H3K9me2 or H3K9me3 [52,53]. In addition, H2A.Z associates with subtelometric regions and other heterochromatin boundaries, where it suppresses the spread of heterochromatin caused by Sir2 deacetylase activity [54]. H2A.Z/H2B.Z were not enriched at centromeric or telomeric-TAS regions in our study, suggesting that these variants have little or no role controlling heterochromatin architecture in T. gondii.
H2A.X is commonly associated with double strand DNA break repair, including in T. gondii [14], but also functions in transcriptional silencing, chromosome segregation, and senescence [55]. H2A.X localizes to different loci than H2A.Z/H2B.Z, and is enriched at the chromosome ends and centromeres, together with H3.3 [11]. In contrast to the literature, our genome-wide study shows that repressed loci are not significantly enriched with H2A.X, and instead, genes transcribed in tachyzoites are marked by H2A.X within gene bodies.
Functional regions of chromatin, telomeres and centromeres, have distinct structures and composition. This study provides insight into the composition of nucleosomes present in these regions in T. gondii (Fig. 7). In eukaryotes, CEN-PA (homologue of CenH3) replaces canonical H3 in centromeric chromatin [56]. In T. gondii, H3.3 and CenH3 both appear to be present in centromeric regions, as confirmed by mass spectrometry (unpublished data). In contrast to T. gondii, P. falciparum H3.3 is not present at centromeric regions [57], suggestive of different chromatin regulatory mechanisms in these two Apicomplexa species. Our findings mirror those in Drosophila melanogaster and humans, where centromeric chromatin was shown to contain CEN-PA containing nucleosomes interspersed with canonical H3-containing nucleosomes [58].
In mammalian cells, CEN-PA is recruited to centromeric regions by formation of double stranded breaks to which H2A.X is also recruited [59]. Given the localization of H2A.X at centromeric regions, this mechanism may be conserved in T. gondii. In addition, H3.3 and H2A.X were both detected at telomeric regions in T. gondii. In mammalian cells, H2A.X is important for the correct arrangement of telomere ends during meiotic prophase [60], while H3.3 is deposited at telomeric regions where it is important for facilitating heterochromatin assembly [61].
As well as histone variants having a potential role in stage-specific gene expression, we also noted a cell Model of different nucleosomal composition on T. gondii chromatin. Different chromatin regions have distinct configurations that differentially change the chromatin architecture according to DNA-related processes. Nucleosomes in typical heterochromatin regions such as centromeres contain the specific centromeric histone CenH3, together with H2A.X. H3.3 enrichment is also evident in those regions, suggesting an atypical heterodimer with CenH3 and H3.3 or alternate CenH3 and H3.3-containing nucleosomes. H2A.X likely forms a tetramer with an unknown H2B, probably H2Ba. Canonical H2Ba with H2A.X is located at opposite sites to H2A.Z/H2B.Z dimer. H3.3 and H2A.X are also enriched at chromosome ends. Transcription start sites (TSS) of active genes are enriched with the variant dimer H2A.Z/H2B.Z. On the CDS of these genes, H2A.Z/H2B.Z are replaced with H2A.X (with H2Ba), and the canonical histone H3 by H3.3. Transcriptionally silenced genes have a similar composition with H2A.Z, H2B.Z and H3.3 at the TSS. But this configuration extends throughout the CDS of these genes, especially at genes which are specific for bradyzoites. These regions do not contain H3K4me3, a landmark for active transcription, suggesting that these genes could be poised. H4 has no variants, so it should be present throughout the genome cycle-dependent pattern of variant histone distribution. In these experiments with asynchronous cultures, which are primarily in G1 phase, genes with less H2A.Z/ H2B.Z coverage (covering only the 5′ end of the gene) are those upregulated at an earlier time point of cell cycle compared to those with more coverage. Genes where the H2A.Z/H2B.Z dimer colocalizes with H3K4me3 on promoters and also localizes on gene bodies, are genes expressed in mid-S/M phase. The timing of the tachyzoite-bradyzoite developmental transition is also linked to the cell cycle [37,38,62], so distribution and PTM modifications of the H2A.Z/H2B.Z dimer could provide a mechanism by which T. gondii integrates signals to either initiate the differentiation program or begin DNA synthesis and cell replication.

Conclusions
Together, this work suggests that histone variants regulate cell cycle and stage-specific gene expression in T. gondii. We hypothesize that variant histones regulate stage-specific gene expression by demarcating functional regions of chromatin. H2A.Z/H2B.Z dimer is implicated as a critical player in regulation of both gene activation and gene silencing. Future work is needed to determine how the localization of histone variants evolves during T. gondii bradyzoite differentiation and the mechanisms by which H2A.Z/H2B.Z and H2A.X/H2B dimers are exchanged.

Cell culture and parasite infection
HFF cells (Human Foreskin Fibroblasts were originally isolated from pooled human foreskins at Stanford University in 1994 and preserved as low passage primary cell line stocks) were grown in Dulbecco's Modified Eagle Medium (DMEM), supplemented with 10% fetal bovine serum, 2 mM L-glutamine, 100 U/ml penicillin and 100 μg/ml streptomycin and maintained at 37 °C with 5% CO 2 . Confluent HFF cells were infected with tachyzoite forms of T. gondii (strain Type I RHΔhxgprt [63]) at a ratio of 3-5 parasites per cell. At around 40 h post infection, cells were harvested and lysed sequentially with 20G/23G/25G needles.
Endogenous epitope tagged variant H2A.X histone strain in the Type I RHΔhxgprt background [14] was subjected to chromatin immunoprecipitation (IP) without previous purification from host cell debris. The H2B.Zmyc strain was generated by cloning the H2B.Z ORF in the plasmid pTUB8mycGFPtailTy-HX in frame with a c-Myc tag at the N-terminal end using selection for HXPRT into the RHΔhxgprt strain. The H3.3-HA was cloned by placing the H3.3 ORF in-frame with a C-terminal HA tag and was also transfected into the Type I RHΔhxgprt background, and stable clones created using selection for HXPRT. For IPs that required antibodies, parasites were purified from host cell debris using a 3.0 μm Nucleopore filter (Whatman) until no human cells were visualized. After two washes with cold PBS, cell pellets were frozen in liquid nitrogen.
Following incubation, cells were washed twice with PBS and incubated with secondary antibodies (Alexa fluor 488-Molecular Probes Cat No A-11034) and 10 μg/mL of 4′,6-diamidino-2-phenylindole (DAPI) to detect nucleic acids. Slides were mounted with Vectashield (Vector Laboratories) and analyzed with a Delta Vision Elite Core Deconvolution microscope. Imaging was performed on an inverted Olympus IX71 microscope with 60x objective (Albert Einstein College of Medicine Analytical Imaging Facility); Z-stacks of 0.1 μm were collected and deconvolution was performed using SoftWoRx 3.6.0. Images were processed with ImageJ Fiji 1.46 J, Adobe Photoshop CS3 and Adobe Illustrator CS3.

Native chromatin immunoprecipitation (N-ChIP)
Freshly collected or frozen parasites (3- Eluted DNA was recovered using MinElute PCR Purification Kit, Cat. no. 28006, Qiagen). The concentration of the purified DNA was determined with the Quant-iT ™ PicoGreen system (Invitrogen). ChIP samples were submitted to the Einstein Epigenomics Shared Facility (ESF) for sequencing. Adapters were ligated to DNA to create the Illumina library and sequenced using a Hiseq2500 sequencing machine (Illumina). This involves a first step of cluster generation in situ on the flow cell followed by massively-parallel sequencing using fluorescently-labeled reversible terminator nucleotides. Clustered templates were sequenced base-by-base during each read. Images generated by the Hiseq2500 were processed by Hiseq Sequence Control Software (HCS V2 with RTA v1.17 and CASAVA 1.8 software packages). These software components perform the functions of image analysis (Firecrest) and base-calling (Bustard).

ChIP-chip data processing
Pair files (IP and input) were combined and processed using NimbleScan software (NimbleGen) according to the software manufacturer's instructions. A minimum FDR of 0.05 was applied to define peaks. Scaled log transformed ratio files (IP/input) were generated in NimbleScan and visualized with IGV software.

Comparison of ChIP-seq and ChIP-chip data
Correlations between ChIP-seq and ChIP-chip were calculated at the peak level using the Bioconductor R package 'diffbind' (http:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ DiffB ind. html). Peak lengths and numbers were calculated from peak files determined from Nim-bleScan (ChIP-chip), MACS2 (ChIP-seq) and Homer (ChIP-seq).

Collective annotation of ChIP peaks
To obtain a list of reproducibly detected peaks that are common to all ChIP-seq and ChIP-chip experiments for H3K4me3 or histone variants, the findOverlapsOf-Peaks function in the Bioconductor R package 'ChIP-peakAnno' (http:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ ChIPp eakAn no. html) was used to generate a merged BED file containing all overlapping peaks with a minimum of 50 bp overlap. Peaks were annotated using ChIPpeakAnno annotatePeakInBatch function using a GFF genome annotation file downloaded from ToxoDB. org (T. gondii ME49, release 26 and limited to 1000 bp distance to a genomic feature. Plots of peak locations were generated using ggplot2.
Prior to downstream analysis, for each histone mark or variant, called peaks from biological replicates were intersected; only peaks that were present in all analyses were combined for assigning genes to peaks. The correlation plots in Fig. 2B were generated using raw ChIPseq files, and the correlation plots in additional file 2C were generated using called peaks from ChIP-seq or ChIP-chip.

Average ChIP-seq profiles over genomic features
Deeptools plotProfile (http:// deept ools. readt hedocs. io/ en/ latest/ index. html) was used to plot composite ChIPseq profiles by averaging the read density of library sizenormalized BigWig files over various genomic regions, with 500 bp upstream and downstream flanking regions.

Enrichment analysis of genes with different histone variant coverage
Lists of genes where histone variants localize were analyzed by pathway and gene set enrichment analysis in R. Gene sets consisting of genes upregulated at different cell cycle time points, GO and KEGG pathways and localization [33] were compared to the histone variant gene lists by a hypergeometric test as previously described [21]. Briefly, p-values of enrichment were determined and adjusted for multiple hypotheses testing by the Bonferroni method. Adjusted p-values were divided by p-values generated by 1000 random iterations of each gene set, generating a normalized p-value used to generate plots. Scripts for enrichment analysis and annotation files can be found at https:// github. com/ natal iesil mon/ toxot ools. p-value of < 0.05 in the hypergeometric test were considered significant.

Motif analysis
DNA sequences from overlapping H2A.Z, H2B.Z and H3K4me3 peaks were extracted and searched for motifs using the commandline version of Homer (http:// homer. ucsd. edu/ homer/ motif/), compared to a FASTA file containing either the entire T. gondii genomic sequence or promoter regions as the background.