Characterization of the accessible genome in the human malaria parasite Plasmodium falciparum

Abstract Human malaria is a devastating disease and a major cause of poverty in resource-limited countries. To develop and adapt within hosts Plasmodium falciparum undergoes drastic switches in gene expression. To identify regulatory regions in the parasite genome, we performed genome-wide profiling of chromatin accessibility in two culture-adapted isogenic subclones at four developmental stages during the intraerythrocytic cycle by using the Assay for Transposase-Accessible Chromatin by sequencing (ATAC-seq). Tn5 transposase hypersensitivity sites (THSSs) localize preferentially at transcriptional start sites (TSSs). Chromatin accessibility by ATAC-seq is predictive of active transcription and of the levels of histone marks H3K9ac and H3K4me3. Our assay allows the identification of novel regulatory regions including TSS and enhancer-like elements. We show that the dynamics in the accessible chromatin profile matches temporal transcription during development. Motif analysis of stage-specific ATAC-seq sites predicts the in vivo binding sites and function of multiple ApiAP2 transcription factors. At last, the alternative expression states of some clonally variant genes (CVGs), including eba, phist, var and clag genes, associate with a differential ATAC-seq signal at their promoters. Altogether, this study identifies genome-wide regulatory regions likely to play an essential function in the developmental transitions and in CVG expression in P. falciparum.


INTRODUCTION
The chromatin structure defines the scenario where the interactions between transcription factors (TFs) and their cognate regulatory regions take place, and this is dynamically shaped in a cell and stage-specific manner. To successfully interact with promoters, enhancers, insulators and non-coding RNAs (ncRNAs), TFs must induce chromatin remodeling of nucleosomal structures, which results in different levels of chromatin accessibility (1,2). Based on this principle, open chromatin profiling can be used to identify regulatory elements associated with specific transcriptional and epigenetic states (3). The identification of these regulatory sequences in their native chromatin environment is important for understanding how gene expression is coordinated throughout development and in response to the environment.
Malaria parasites show a complex life cycle in the human and the mosquito hosts and face dynamically changing environments during development. Plasmodium falciparum is the most virulent and prevalent human malaria parasite species in Africa (4). In order to successfully develop within the host, this parasite switches transcriptional programs between stages (5)(6)(7) and has evolved mechanisms to adjust its phenotype through transcriptional heterogeneity within populations. The ability of the parasite to survive in the face of changes in host conditions is tightly linked to the variant gene expression of a number of genes involved in processes such as antigenic variation, red blood cell invasion, solute transport and sexual differentiation (8,9). These genes show clonally variant gene (CVG) expression, such that individual parasites having identical genomes and under the same environment can maintain a variant gene in a different transcriptional state, and this state can be transmitted to the next generation by epigenetic mechanisms. Multiple gene families showing clonally variant expression have been characterized (6,10), many of which encode anti-gens expressed at the surface of infected erythrocytes, such as the multicopy var, rifin, stevor, surfin and Pfmc-2TM CVG families (10,11). Among them, the var genes encode the Erythrocyte Membrane Protein 1 (PfEMP1), which is a critical virulence factor for malaria. Each individual parasite has ∼60 different var genes, only one of which is expressed at a time, i.e. mutually exclusive expression (12)(13)(14). Another clonally variant family is clag, in which two of the genes (clag3. 1 and clag3.2) are also expressed in a mutually exclusive manner (15). They encode for membrane proteins involved in solute transport and are linked to drug resistance (16)(17)(18)(19)(20).
Despite the implications for malaria pathogenesis, the role of chromatin-mediated processes in the regulation of gene expression remains poorly understood. This lack of knowledge is in part due to unique particularities of P. falciparum in relation to genomic architecture. For example, one key feature is its extreme AT richness, which is the highest compared to all genomes sequenced to date: >80% on average and often reaching 90-95% in intergenic regions (21). Recent models suggest that highly AT-rich sequences may inhibit the formation of higher-order structures maintaining chromatin accessible and the genome in a transcriptionally permissive state (22,23). However, during development, the parasite undergoes drastic switches in the pattern of gene expression (5,24,25). Active chromatin regions are marked primarily by H3K4me3 and H3K9ac (26,27), whereas the non-active fraction of the genome corresponds to blocks of heterochromatin that are marked by histone 3 tri-methylation at lysine 9 (H3K9me3) and the heterochromatin protein HP1 (28). Importantly, CVGs locate primarily in facultative heterochromatin regions occupying central and subtelomeric clusters. The question then arises: how these genes are switched on/off in such a heterochromatic environment?
A growing hypothesis in the field is that promoter-specific temporal regulation of transcription in P. falciparum occurs through site-specific recognition of DNA sequences by TFs that are uniquely produced in each stage of development. Several cis-regulatory elements involved in the temporal regulation of transcription have been predicted bioinformatically, but in general they are still poorly understood (29)(30)(31)(32). Furthermore, despite some work having been done on trans-acting regulatory DNA sequences (33), the characterization and function of enhancer elements remain a poorly explored issue in the field. In P. falciparum there is only one large family of specific TFs characterized to date, the ApiAP2 family. Current evidence points to an essential role of these DNA regulatory proteins in driving P. falciparum developmental progression (34)(35)(36)(37)(38)(39)(40)(41)(42)(43)(44). Accordingly, the genes encoding these TFs are expressed in a stage-specific manner (45,46). The cognate motifs of several ApiAP2 TFs have been determined experimentally (47), but only a few of these have been tested in vivo by chromatin immunoprecipitation-sequencing (ChIP-seq) (41,42), mainly because of the difficulties in obtaining antibodies against Plasmodium proteins (47). Profiling of accessible chromatin can be used to reveal gene regulatory networks and to identify novel trans-acting factors and cis-regulatory elements. MNase-seq has been employed to characterize nucleosome-depleted regions at different time points dur-ing the intraerythrocytic developmental cycle (IDC) (48)(49)(50). This technique is useful for identifying genomic regions occupied by nucleosomes or accessible to nuclease cleavage. However, this data alone provides limited information regarding the identity and function of regulatory elements (51). FAIRE-seq is another method for mapping chromatin accessibility that has been used in Plasmodium but has low resolution and limited accuracy in identifying DNA-protein binding events. In consequence, many of the TF binding sites in P. falciparum have only been predicted bioinformatically (29)(30)(31)(32). In this context, the Assay for Transposase-Accessible Chromatin by sequencing (ATACseq) has the advantage of allowing simultaneous mapping of nucleosomes and DNA-binding proteins protected sites, being much faster and more sensitive and requiring orders of magnitude less starting material than other assays (3).
In this work, we use ATAC-seq to characterize the chromatin accessibility landscape during the intraerythrocytic cycle in two transcriptionally variant subclones. Our results show that differences in chromatin accessibility recapitulate transcriptional changes during development and are predictive of clonal differences in epigenetic and expression profiles. The analysis of differentially accessible regions allowed us to identify the landscape of binding events, regulatory DNA sequences and putative TFs that are likely responsible for these changes. Altogether, our work shows the dynamics of the P. falciparum regulatory genome and its association with developmentally regulated and CVG expression.

Parasite cultures
10G and 1.2B correspond to P. falciparum subclones derived from the same 3D7-A stock, which have been characterized in previous studies (6,15,52,53). The experiments were performed at different times but on the same laboratory and using the same parasite stocks. Parasites were thawed from a cryopreserved stock and cultured for two or three generations in B+ erythrocytes (3% hematocrit) under standard conditions with media containing Albumax II and no human serum. Cultures were tightly synchronized to a defined 5 h age window by purification of parasites at the schizont stage using Percoll gradients (63% Percoll), followed by sorbitol lysis 5 h later to eliminate erythrocytes infected with late asexual stages (trophozoites and schizonts). Parasites were collected at the following times post invasion: 10-15 hpi (T10, early ring stage), 20-25 hpi (T20, late ring stage), 30-35 hpi (T30, trophozoite stage) and 40-45 hpi (T40, schizont stage). Two independent biological replicates were performed. gDNA was isolated from parasites at the schizont stage using the DNeasy Blood & Tissue Kit (Qiagen), according to the manufacturer's instructions.

ATAC-seq
The ATAC-seq protocol was performed using 10 million cells for early ring, late ring and throphozoite cell stages, and 1 million cells for schizonts, which were obtained after saponin red blood cell lysis. Parasite cells were resuspended in lysis buffer in order to permeabilize membranes. The nuclei pellet was immediately resuspended in the transposition reaction mix (25 l of 2× TD Buffer, 1.25 l of Tn5 Transposase and 23.75 l of nuclease free water), and incubated for 30 min at 37 • C. As a control, 50 ng of gDNA were used for the standard ATAC-seq procedure as described in the Nextera library preparation protocol by Illumina (5 min at 55 • C). All transposed samples were purified using the Qiagen MiniElute Kit. Following purification, ATAC-seq library amplification was carried out with 2× KAPA HiFi mix and 1.25 M of Nextera primers (3). The use of this enzyme eliminates mapping biases due to the high AT content of the P. falciparum genome (54,55). The optimal cycle number was determined using quantitative polymerase chain reaction (qPCR) with PCR conditions as originally described in (3). ATAC-seq libraries were sequenced at BGI (China) using an Illumina HiSeq2000 sequencer to obtain 20-40M of 2 × 50 bp paired-end sequencing reads per sample (Supplementary Table S1).

Data analysis
All sequenced ATAC-seq paired reads were trimmed 10 bases from each read 3' end (-3 10) and aligned to the P. falciparum reference genome Pf3D7-28 (PlasmoDB) with Bowtie2 (v.2.3.1) using default parameters except for -nounal -no-mixed -X 2000. Alignment files were filtered applying a quality threshold of 10 in the MAPQ score, sorted and deduplicated using Samtools (v.1.4). To adjust for fragment size and to make sure Tn5 cutting sites were mapped, aligned reads were shifted +4 bp for + strands and −5 bp for − strands. The extremely high AT bias of the P. falciparum genome might reduce mapping quality and introduce coverage biases. Furthermore, the presence multiple copy gene families, such as var and clag3, can lead to multiple sites mapping and the removal of non-unique reads.
To check the impact of this mapping biases, we ran Bowtie2 with no quality threshold and found a high similarity in the percentage of ATAC-seq coverage at var and clag3 families, with <10% difference between the two quality filtered alignments.
ATAC-seq data analysis was conducted following the prototype ATAC-seq ENCODE Uniform Processing Pipeline prototype (https://www.encodeproject.org/atacseq/). First, nucleosome-free fragments were extracted by applying a size threshold of 130 bp. Peak calling on nucleosome-free reads was performed by MACS2 (v.2.1.1) (63) using the module 'callpeak' with the following parameters: -c -g 2.41e7 -keep-dup all -q 0.01 -nomodel -shift -35 -extsize 75 and -B. Transposed genomic DNA (gDNA) was used as a control for the peak calling. For this analysis, bam files were downsampled (Samtools v.1.4) to equal the number of reads between treatment and input (gDNA control) for each ATAC experimental sample (Supplementary Table  S1). We applied a replicability filter to the output of the MACS2 peak calling analysis using the 'intersect' tool from the BEDTools software suite (v.2.25.0) (64). In addition to the independent peak calling for each biological replicate, we performed peak calling on the pooled sample and the pseudoreplicates, and kept only consensus ATAC-seq peaks if present in the three datasets (replicate, pooled and pseudoreplicate). Peaks located in mitochondrion and apicoplast chromosomes were removed. The resulting set of ATAC-seq peaks, that we considered high confidence Tn5 hypersensitive site (THSS) regions, are included in Supplementary Table S2. For visualization purposes, tracks of input (gDNA) corrected ATAC-seq signal per stage were built with MACS2 'bdgcmp' module (-m ppois) on each pair of fragment pileup and control lambda bedGraph files from the peak calling analysis.
Annotation of THSS regions to genomic features was performed using the HOMER software (v.3.12) (65). We used the gene ensemble available from PlasmoDB (Pf3D7-28) and considered the translation start site ATG as the reference point and −1 Kb upstream as the putative promoter region. This assumption is based on a recent study that mapped transcription start sites (TSSs) in P. falciparum at high resolution and revealed that 81% of the TSSs are positioned less than 1 Kb upstream of the start codon, with more than half of these located at distances less than 0.5 Kb (66). For the ATAC-seq enrichment analysis we used BEDTools to obtain the number of reads overlapping with: known TSSs (48,66), promoters (−1 Kb of the translation start codon, ATG), introns, exons, etc. The resulting read counts were normalized, ATAC to noise signal was corrected (ATAC/gDNA ratio) and log2-transformed using R (v.3.4.1) (http://www.r-228 project.org/). We added a pseudocount (0.1) to avoid infinite values when control/input correcting (dividing by 0 in ratios) or when computing log2 values.
In order to assess the global relationship between the ATAC-seq signal at promoters and the transcription levels of the corresponding gene, we categorized genes into high, medium or low groups, based on their mRNA levels using previously published RNA-seq data (48). For this purpose, we applied a threshold value determined by dividing the mRNA values (RPKM normalized values) in three quantile groups according to their means (cut2 function in the 'Hmisc' R package) (67). Average profile plots representing ATAC-seq and histone modifications enrichments (RPKM normalized and input-corrected) centered on TSSs or ATAC peaks coordinates, were built using ngs.plot (v.2.61) (68). Violin plots, boxplots and line diagrams were produced using 'ggplot2' R package (69). For comparative and visualization purposes, ATAC-seq, RNA-seq and histone ChIP-seq enrichment data at various genomic features are shown on a log2 scale. Correlation tests were performed between mRNA levels, ATAC-seq and H3K9ac/H3K4me3 ChIP-seq enrichment data using a Spearman rank correlation test in R.
Unless if otherwise specified, interval operations like intersect, shuffle, merge, flank or slop were performed using BEDTools and statistical tests and plots were performed in R using Bioconductor (http://www.bioconductor.org).

Characterization of novel ATAC-seq regulatory regions
We classified novel THSS elements, ATAC-seq peaks located in intergenic regions, into TSS or enhancer-like regions using published H3K4me1, H3K4me3 and H3K9ac ChIP-seq data (59). Since H3K4me1 data exists only for trophozoite and schizonts, we restricted our analysis to these two stages. We used ChromHMM (70) to compute chromatin state predictions on the entire P. falciparum genome based on relative enrichment levels of these three histone modification marks. For the binarization of the genome, default values were used for all parameters except for -b 75. We chose a four states model as a prior assuming the following chromatin states: unspecific (high levels of enrichment for all hPTMs), depleted (low levels of enrichment for all hPTMs), TSS (H3K9ac/H3K4me3 enrichment and H3K4me1 depleted) and enhancer (H3K4me1 enrichment and H3K9ac/H3K4me3 depleted). According to the ChromHMM segmentation, most of the genome displays low enrichment signal (depl). In order to obtain a high-confidence set of regulatory regions, we first filtered out THSSs located in depleted and unspecific enrichment signal genome segments. We then ran again ChromHMM model on the filtered dataset and discarded those regions not classified as TSS-like or enhancer-like in the second prediction.

Soft clustering analysis
In order to analyze chromatin accessibility dynamics among different stages, we conducted a soft clustering analysis using the 'Mfuzz' R package (71) on THSS-annotated promoters. Using a standard m fuzzy c-means parameter of 1.7, a total of 30 pre-clusters were created that summarize the variability of temporal profiles in our data. Promoters showing a pertinence membership >51% and a unique peak of chromatin accessibility at a given stage were grouped into four Mfuzz clusters with promoters showing maximum ATAC-seq enrichment levels at early rings, late rings, trophozoite and schizont parasite stages, respectively (Supplementary Figure S4). For visualization purposes, the promoters in each cluster were ordered by ATAC-seq enrichment levels at corresponding stage, and we then used this clustering order to show the association between the dynamics in chromatin accessibility, gene expression and histone modifications enrichment. For this analysis we used the microarray normalized counts obtained in a previous study of the same subclones and stages of development (6). ChIP-seq enrichment data for H3K4me3 and H3K9ac was obtained from (26).

Differential ATAC-seq enrichment analysis
Of the 5085 total peaks, 1053 were unique to 10G and 1797 to 1.2B parasites (Supplementary Figure S7A). However, many of the ATAC-seq peaks unique to one of the two subclones are likely attributable to some variability in the performance of the MACS2 peak calling algorithm (e.g., some low enriched peaks may be similar between the two subclones but pass the peak calling criteria in only one) (Supplementary Figure S7B). To account for this issue, differential ATAC-seq accessibility analysis between 1.2B and 10G subclones was conducted using the 'DESeq2' R package (58). This analysis takes all the ATAC-peaks called in 10G and 1.2B, and tests for normalized read count differences at the peak region. In analyzing clonal differences for each developmental stage separately, we detected that a fraction of the clonal variation was due to slight differences in age between subclones. For this reason, we repeated the analysis treating the stage of development as a co-variable and we used this set of clonally variant peaks for all subsequent analyses. Heatmaps showing ATAC-seq enrichment at promoters and microarray gene expression data for the clonally variant genes were built using the 'iheatmapr' R package (72) and the heatmap.2 function in the 'gplots' R package (73). In order to visualize quantitative differences in gene expression patterns between subclones we used microarray data from (6) and considered the log2 ratio relative to a reference pool. Due to very low transcription levels for most gene copies in multicopy gene families like var, the expression of a given gene in the reference pool could be higher than in the sample. To control for this, we obtained semi-quantitative estimates of var gene expression patterns considering microarray values for the sample channel (F635) instead of the log2 ratio relative to a reference pool when building the heatmaps. Furthermore, the telomeric PF3D7 1100100 var gene was excluded from the analyses because there is a deletion in the sequence in the 10G subclone. For visualization purposes, microarray expression was log2-transformed. Promoter windows for var genes were defined as 2 Kb upstream of the ATG sites if possible, or as the intergenic region not overlapping with exons or promoters of nearest adjacent genes. In cases where over-lapping adjacent genes were of opposing directions, half the intergenic distance was assigned to each gene.

Motif discovery
We carried out de novo motif analysis using HOMER software (v.3.12) (65) on the stage-specific THSS regions obtained by the soft clustering analysis. For this analysis, ATAC peaks specific of each time point (T10, T20, T30 and T40) were merged and the highest scoring summit was selected as summit for the merged peak. We then used the 'slop' tool from BEDTools software suite to obtain equal size 200 bp peaks from the summit coordinate. This analysis allowed us to identify DNA motif sequences enriched in the set of active regulatory sites during the intraerythrocytic development. Next, in order to investigate whether predicted binding sites for ApiAP2 TFs are among the motifs identified by HOMER, we performed similarity domain analysis with the TomTom tool (http://meme.nbcr.net/meme/doc/ tomtom.html) using the previously published database of validated ApiAP2 motifs (47). Only motifs enriched in more than 5% of the targets sequences and below a threshold Pvalue of 10e-10 were considered, and results corresponding to low complexity motifs and offsets or degenerate versions of highly enriched motifs were discarded. For the Tom-Tom analyses, those hits that aligned to the core and not to the edge of the known ApiAP2 motif, and met a minimum P-value of 10e-3, were selected. We used the annotate-Peaks.pl module in HOMER to find motif occurrences in each stage-specific ATAC peaks sets. A Cytoscape network analysis (http://www.cytoscape.org/) was conducted to investigate links between significantly enriched DNA motif sequences at stage-specific regulatory active sites and their target genes.

Real-time quantitative PCR
Microarray expression data of 10G and 1.2B parasite lines from a previous study (6) used for comparison with the ATAC-Seq data, was validated for a subset of CVGs (Supplementary Figure S9). Total RNA from in vitro cultures of the same parasite stocks was obtained by TRIZOL extraction followed by reverse transcription using the Tetro cDNA Synthesis Kit (Bioline), according to the manufacturers' instructions. Amplification was performed using Roche Ligh-Cycler 480. Reactions (12.5 l) were performed in duplicate and consisted of 0.5 l of 1/10 cDNA dilution, 0.2 M of forward and reverse primers and 5 l of KAPA SYBR FAST qPCR Master Mix (KapaBiosystems). Cycling conditions were 10 min initial denaturation at 95 • C, followed by 40 cycles of 15 s denaturation at 95 • C, 30 s annealing at 57 • C and 30 s extension at 60 • C. Melting curve analysis was performed to guarantee specificity of the template. The amount of cDNA was normalized using the housekeeping gene S-adenosylmethionine synthetase (PF3D7 0922200). The 2 − CT method (74) was used to estimate relative mRNA abundance.

Chromatin accessibility by ATAC-seq is predictive of active transcription and epigenetic states in P. falciparum
We generated ATAC-seq libraries using tightly synchronized parasites at four different developmental stages during the IDC: early rings (10-15 h, T10), late rings (20-25 h, T20), trophozoites (30-35 h, T30) and schizonts (40-45 h, T40) (Supplementary Table S1). For each stage, we obtained two independent ATAC-seq replicates and performed experiments in parallel with two different parasite lines, 1.2B and 10G ( Figure 1A). These lines correspond to two subclones of the same 3D7 genetic background that differ in the expression of several CVGs (6,15). The similarity between replicates is high compared to the similarity between subclones and between stages (Supplementary Figures S1A  and B). The similar distribution of fragment sizes suggests that chromatin is accessible to Tn5 transposase to the same degree in all samples independently of the stage of development ( Figure 1A). The identified THSS regions in all samples are listed in Supplementary Table S2. The number of peaks and the intensity of the ATAC-seq signal is higher in the trophozoite stage (10G: 2229, 1.2B: 2158), followed by schizonts (10G: 1062, 1.2B: 1080), late rings (10G: 950, 1.2B: 2223) and early rings (10G: 613, 1.2B: 1042) (Supplementary Figure S1C and Table S2). The number of peaks in 1.2B ring stage parasites is higher compared to 10G, but the ATAC-seq enrichment levels are similar. Such differences are possibly due to slightly lower peak calling performance in the 10G subclone ( Supplementary Figures S1C and 7B).
Distance analysis of ATAC-seq data shows that THSS regions are mostly located in a window of 1 Kb upstream from the translation start codon (ATG) and are scarce at distances greater than 2 Kb of the nearest gene ( Figure 1B and Supplementary Figure S1D). This is in agreement with previous studies in P. falciparum showing that TSSs preferentially localize 0.5-1 Kb upstream from the ATG (48). The ATAC-seq signal is significantly enriched at previously annotated TSSs (48,66) and matches the typical profiles described in higher eukaryotes (3,75,76) (Figure 1C and Supplementary Figure S1E). Based on these findings, we annotate THSSs as promoters if the peak is located less than 1 Kb upstream from the corresponding ATG ( Figure 1D and Supplementary Figure S1F); and use this annotation for all subsequent analysis (Supplementary Table S2).
It is known that chromatin features such as histone modifications exert a significant impact on the compaction of chromatin (77,78). In Plasmodium, as in other eukaryotes, active histone modifications such as H3K9ac and H3K4me3 have been shown to be enriched in intergenic regions and colocalize with H2A.Z (26). The association of these histone modifications with transcription varies through the IDC: H3K9ac has been shown to strongly correlate with transcription all over the asexual cycle, whereas H3K4me3 is enriched in schizonts regardless of the temporal expression dynamics. Thus, we expected a genomewide stage-dependent relationship between these active chromatin features and chromatin accessibility by ATACseq. Using publicly available data (26), we observed that H3K9ac, H3K4me3 and H2A.Z, are enriched around the THSSs. Furthermore, the association between these epigenomic features is stage-specific, being more marked in schizonts ( Figure 1E and Supplementary Figure S2A).
Following the observation that the majority of THSS regions overlap with active histone marks at TSSs, we examined the quantitative relationship between chromatin accessibility and transcription genome-wide. To this end, using a public RNA-seq dataset (48), we computed the ATACseq signal level at known TSSs (48) in groups of genes with low, medium, and high mRNA abundance. We observe that the level of enrichment in ATAC-seq signal at the promoter is positively associated with mRNA levels of the annotated gene ( Figure 1F and Supplementary Figure S2B). The same pattern is observed when categorizing promoter regions by the strength of the ATAC-seq signal into low, medium and high groups (same grouping procedure as above for the RNA-seq data above) and representing the mRNA levels of THSS-annotated genes in each ATAC-seq enrichment group ( Figure 1F and Supplementary Figure S2B).

Chromatin accessibility footprints identify novel regulatory elements in P. falciparum
Chromatin accessibility represents a novel and useful tool to identify novel regulatory sequences including TSSs and enhancers (79). Then, different histone post-translational modifications (hPTMs) can be used to distinguish between these elements. For example, in higher eukaryotes, histone marks H3K4me3/H3K9ac are characteristic of TSS/promoter regions, whereas H3K4me1 is a hallmark of enhancer regions, which are typically depleted in H3K4me3/H3K9ac (80)(81)(82). Based on this assumption, we investigated whether the THSS regions correspond to TSS-like or to enhancer-like regions. To this end, we first used a chromatin-state segmentation (ChromHMM) approach (70) to partition the P. falciparum genome based on the relative enrichment of several reported histone marks (59). Since available data for H3K4me1 exists only for trophozoite and schizonts (59), we restricted our analysis to 2273 merged THSSs from these two stages in both subclones (Supplementary Table S2 Table S3).
The majority of the ATAC-seq peaks with TSS-like signatures (523/558) are located <2 Kb from a downstream transcript ( Figure 2E), whereas there are 35 TSS-like regions that locate more than 2 Kb upstream (Supplementary Table  S3). As expected many of the TSS-like elements (269 out of the 558) coincide with a TSS previously described by others (48,66). As for the rest, we find 250 THSSs located in the promoter of genes with a known TSS, but that does not coincide with the ATAC-seq peak. We also find 12 regions that are located in the promoter of genes with no previous anno-tated TSS, and 27 putative novel TSS-like elements located in non-promoter regions (>2 Kb upstream of the ATG). For the 289 putative novel TSSs, the positive and significant correlation between the ATAC-seq signal enrichment and the mRNA levels of the associated gene in trophozoites suggest that these THSSs are likely to be previously unidentified TSSs ( Figure 2F and Supplementary Figure S3A).
Contrary to the TSS-like elements, 24 out of the 50 enhancer-like regulatory elements are located at distances >2 Kb upstream of the ATG ( Figure 2E and Supplementary Table S3). There are 12 out of the 50 enhancer-like regions that coincide with known TSSs (48,66). Amongst the remaining enhancer-like elements, we find 35 that are located upstream of genes with one or several known TSSs. Most of the known TSSs (161/172) that are associated with one of the novel enhancer-like element are not accessible. The more distal enhancer-like regions are commonly found at the ends of chromosomes, occupying subtelomeric repeats such as the TAREs and sometimes very far from the nearest gene. This is likely because TAREs are highly transcribed (and thus presumably contain promoters and enhancers) (42,83). Interestingly, we find that 22 subtelomeric var genes have enhancer-like THSS regions. At last, the study by Ubhe et al. (33) proposed a set of enhancer-like sequences based on the localized enrichment of H3K4me1 versus H3K9ac and H3K4me3 in intergenic regions. However, we do not observe a noticeable ATACseq enrichment around those sites, and none of the intergenic regions classified as enhancer-like coincide with our enhancer-like dataset (Supplementary Figure S3B). Contrary, 12% (57) of the enhancers described by Ubhe et al. intersect with 72 THSSs classified as depleted (44), unspecific (6) and TSS-like (22) by the ChromHMM analysis. We consider that our analysis combining chromatin accessibility profiling and histone modification enrichment is a much more precise and powerful approach in identifying new trans and cis-regulatory sequences in P. falciparum.

Temporal variation in chromatin accessibility is predictive of stage-specific transcription
Plasmodium falciparum undergoes drastic switches in transcription linked to life cycle progression. Our results indicate that chromatin accessibility is dynamic and varies between stages of development ( Figure 3A and Supplementary Figure S1C). To further investigate the relationship between temporal changes in chromatin accessibility and transcription through the cell cycle, we applied a soft clustering approach using Mfuzz that group promoters based on their temporal accessibility through the cell cycle. In this analysis we observe clusters of dynamic promoters accessible at different stages during development (Supplementary Figure S4). We then determined the temporal dynamics in mRNA levels of the corresponding promoter-annotated gene (6) Table S4). The Gene Ontology terms enriched at each   (6) for the set of THSS-annotated genes ordered by stage-specific cluster and ATAC-seq enrichment as in (C). Data in the heatmap is log2-scaled and mean centered. The Spearman rank correlation coefficient (rho) and the corresponding P-value are shown above columns for the association between the ATACseq enrichment at the promoter and mRNA levels (6) of the nearest downstream genes at each stage. (E and F) Heatmaps based on levels of H3K4me3 (E) and H3K9ac (F), using ChIP-seq data obtained from (26), at promoters for the set of THSS-annotated genes ordered by stage-specific cluster and ATAC-seq enrichment as in (C). Enrichment data in the heatmaps is log2-scaled and mean centered. The Spearman rank correlation coefficient (rho) and the corresponding P-value are shown above columns for the association between ATAC-seq levels and enrichment of H3K4me3 (E) and H3K9ac (F) at promoters at each stage. (G) Changes in chromatin accessibility and gene expression in the aarp gene (PF3D7 0423400). Tracks show normalized and input-corrected ATAC-seq signal. The location of TSSs (48,66) and the binding sites for AP2-I and BDP1 (41,62) are shown. All tracks are shown at equal scale. The plot at the right shows temporal changes in microarray mRNA levels (log2 scale) (6) during the IDC for this gene. stage-specific cluster agree with the gene function prediction expected based on previous transcriptomic and proteomic analyses (25,84) (Figure 3C; Supplementary Figure  S5B and Table S4). On the other hand, we also detect an association between accessibility and the level of enrichment of H3K9ac throughout the cell cycle. The coefficient is lower but significant for the correlation between ATAC-seq and H3K4me3 enrichment in all stages except at T40 (Figures 3E and F; Supplementary Figures S5D and E). There are multiple examples of this nature, such as the gene encoding the apical asparagine-rich protein (PF3D7 0423400) with an ATAC region only in schizonts that correlates with the expression of the gene that peaks at this developmental stage ( Figure 3G). We also report promoters that contain multiple accessible regulatory sequences opening at different times during development (Supplementary Table S4). For example, the gene encoding the early transcribed membrane protein 5 (ETRAMP5) displays an invariant peak less than 1Kb upstream, and two peaks: one distal >4 Kb upstream and one proximal at the 5' start of the gene, that correlate with the transcriptional status of the gene (Supplementary Figure S5F). All this considered, we show that by combining ATAC-seq and expression data during the IDC cycle, we can identify functionally related genes that show dynamic and coordinated chromatin accessibility and gene expression at different developmental stages.

Stage-specific motif inference at THSS regions identifies key developmental regulators
The dynamics in chromatin accessibility at TSSs and enhancers has been shown to capture DNA-protein binding events that lead to promoter activation in vivo (2,85). To characterize the regulatory networks controlling developmental transitions during the IDC, we first performed de novo motif discovery and similarity analysis on the set of developmentally dynamic promoters identified above (Mfuzz clusters, Supplementary Table S4). Amongst the top enriched motifs, we find potential binding sites for several ApiAP2 TFs including AP2-SP/AP2-EXP (10 and 40 h) and AP2-I (30 and 40 h), among others (48,66). However, we also report de novo predicted motifs ( Figure 3C; Supplementary Figure S5B and Tables S4 and 5). In general, there is a correlation between the developmental abundance of the predicted TF and temporal enrichment of the corresponding binding site ( Figure 3C; Supplementary Figure  S6A and Table S5) (47). Next, we computed the occurrences of each stage-specific motif in the set of ATAC-seq peaks that are associated to each Mfuzz cluster (Supplementary  Table S5), and used Cytoscape to construct gene regulatory networks that connect TF predicted binding sites with the target genes ( Figure 4A and Supplementary Figure S6B). This analysis reveals co-regulated genes and interactions between predicted TFs. For example, looking at the regulatory network of late blood stage parasites (T40), we identify a number of invasion genes like rap, ama, rama and members of the msp, clag, rhop and Rh gene families that contain one or multiple T40-specific DNA motifs in the promoter (Figures 4A-C). We also observe that the AP2-I cognate motif GTGCACGC (named 40.1 in Supplementary Table S5) plays a central role in this gene network displaying the high-est number of connections, i.e. 178 genes contain the motif (Figures 4A and B). These observations suggest that different TFs might act coordinately to control their gene expression.
As a validation for our strategy to investigate TF binding events in vivo, we report the overlap between the ATAC-seq peaks and the binding sites of the TFs: AP2Tel, SIP2, AP2-I and other proteins such as BDP1, which have been determined experimentally (35,41,42,62) (Supplementary Table  S6). In the case of AP2-I, we find that 86% of the binding sites intersect with ATAC-seq peaks (Supplementary Table  S6). Importantly, the enrichment of ATAC-seq peaks at the binding site of AP2-I occurs only in late stages, when the activation of invasion genes takes place ( Figure 4B and Supplementary Figure S6C). Together with AP2-I, the Bromodomain BDP1 has been described as a transcriptional activator that binds to chromatin at invasion gene promoters (62). We find 402 ATAC-seq peaks that intersect with 273 out of 799 mapped binding sites (Supplementary Table S6). An example of a target gene co-regulated by BDP1 and AP2-I is the merozoite surface protein 2 (MSP2) gene, that in our study displays a THSS region in the promoter with the GT-GCA motif ( Figure 4C).
All this considered, the temporal analysis of ATAC-seq data allows us to infer the regulatory DNA sequences, the TFs and the target genes that constitute the regulatory networks controlling gene expression during blood stage developmental transitions. Our results support a model in which the promoter region can be bound by one or multiple stagespecific TFs that would operate coordinately to activate different functional sets of genes.

Differential chromatin accessibility associated with clonally variant gene expression
The 10G and 1.2B subclones of 3D7 have previously been shown to differ in the expression of several CVGs (6). We used DESeq to test for clonal variation in ATAC-seq enrichment using a 2-fold threshold for the difference between subclones. Using this criterion, we identified 140 statistical significant differentially-accessible regions between the 1.2B and 10G lines, of which 38 ATAC-seq peaks (15 specific for 1.2B and 23 for 10G) are located in gene promoters (−2 Kb from the ATG, Supplementary Table  S7). Genes with differential ATAC-seq peaks include wellcharacterized CVGs such as clag2, clag3, eba140 and members of the var and phista families (Figures 5A-C; Supplementary Figures S7C-E and Tables S7 and 8) (6,8).
Taking the list of the top-30 most differentially expressed genes between 1.2B and 10G (6), we observe a clear and significant association between differential gene expression and differential accessibility, such that an increase in mRNA levels in one subclone is generally associated with an increased accessibility in the promoter ( Figure 5B). For example, the promoters of the well-characterized clonally variant invasion genes eba140 and clag2 genes, or a clonally variant phista (15,86), display various open regulatory regions only in the subclone in which the gene is highly expressed ( Figure 5C; Supplementary Figure S7E and Table  S7). However, we also find examples of top CVGs in which chromatin accessibility remains invariant despite the gene being variantly expressed. This is a set of 14 genes from multicopy gene families: nine rif, four var and one pfmc-2TM genes and five genes from other families. Except for the var genes (see below) the majority of these cases are genes that are not accessible in either clone (Supplementary Table S8), suggesting that for these genes differential expression is attributable to activation of the gene in only a small fraction of the parasites in one of the subclones.

Dynamics of regulatory regions in mutually exclusively expressed clag3 and var genes
The clag3 genes are an example of CVGs that are expressed in a mutually exclusive manner, so only one copy is expressed at a time in individual parasites, either clag3.1 or clag3.2 (15,16,(86)(87)(88). These genes are also developmentally regulated, thus when active, they are maximally expressed at the schizont stage. Our results show that the variant expression of clag3 genes is associated with differential  (6). Tracks show normalized/inputcorrected ATAC-seq signal in the 1 Kb upstream region. The location of TSSs (48,66) and the binding sites for AP2-I and BDP1 (41,62) are included. All tracks are shown at equal scale. The plot at the right shows mRNA levels (log2 scale) (6) for each subclone during the IDC. chromatin accessibility in the promoter region. The clag3.1 paralog (PF3D7 0302500), which is expressed in the 1.2B line and silenced in 10G, shows a clear opening in the promoter that coincides with a known TSSs block (48,66,86) only in the 1.2B subclone ( Figure 6A). The considerably smaller ATAC-seq peak observed in the 10G subclone likely corresponds to a subpopulation of parasites that switched from clag3.2 to clag3.1 expression during normal growth (88). The clag3.1 active regulatory region has motifs similar to binding sites for the ApiAP2 PF3D7 1222400, and also overlaps with ChIP-seq peaks for AP2-I and BDP1 (Figure 6A and Supplementary Tables S4 and 5). In an analogous situation, the clag3.2 paralog (PF3D7 0302200) is expressed only in 10G and accessibility at the promoter region is only observed in this subclone. Of the two THSS regions in the clag3.2 promoter, the proximal peak coincides with a  (48,66,86) and the binding sites for AP2-I and BDP1 (41,62) are shown. A DNA motif present in a stagespecific regulatory active site is included. A track with THSS regions with characteristics of novel TSS-like elements as in Figure 2 is shown. All tracks are shown at equal scale. The plots at the right represent microarray mRNA levels (log2 scale) (6) at each subclone during the IDC. (C) Chromatin accessibility profiles in the region of chromosome 7 containing the active var gene (PF3D7 0711700) that is expressed predominantly in 1.2B parasites. Two clonally variant accessible RUF6 encoding genes (PF3D7 0711600, PF3D7 0711800) flanking the 1.2B specific var are shown. Tracks show normalized/inputcorrected ATAC-seq signal. The location of known TSSs (48,66) and the binding sites for SIP2 and AP2Tel (35,42) are included. All tracks are shown at equal scale. The plot at the right corresponds to mRNA levels data (log2 scale of microarray data sample channel, F635) (6) of the PF3D7 0711700 var at each subclone during the IDC.
Another gene family that shows mutually exclusive expression is var. Previous studies using microarrays showed that the 1.2B subclone expresses only one predominant var gene, in contrast to the 10G subclone in which a dominantly expressed var was not identified (6) (Supplementary Figure  S8A). The active var gene in 1.2B PF3D7 0711700 displays several ATAC-seq peaks in the putative promoter region, the intron, exon 1 and the downstream region ( Figure 6C).
Nucleic Acids Research, 2018, Vol. 46, No. 18 9427 Only the two more distal peaks occur specifically in 1.2B, and both coincide with annotated members of the RNA of Unknown Function 6 (RUF6) gene family. This family includes 15 genes that encode for long ncRNAs known as GC-rich elements (89,90). Of these, only ruf6-1 and ruf6-2, which are located flanking the 1.2B active var, become uniquely accessible in 1.2B trophozoites. All other members of the RUF6 family are non-accessible in 1.2B except for ruf6-15, located downstream of var PF3D7 1240900 and in which the upstream ruf6-14 gene is not accessible and ruf6-10 that flanks the 1.2B active rif PF3D7 0808900, one of the top-30 differentially expressed genes between 10G and 1.2B (6) (Supplementary Figure S8B).
Other var genes, which are silenced in the majority of parasites in both the 10G and the 1.2B subclones, also display numerous ATAC-seq peaks, some of which intersect with the binding sites of SIP2 and AP2Tel characterized in previous studies (Supplementary Figure S8C and Table S6) (91,92). The peaks that do not correlate with the active state of var genes may be related with the constitutive promoter activity of the intron (93,94).
Collectively, the analysis of ATAC-seq data supports a model of transcriptional regulation for clag3 and var genes based on the differential use of multiple regulatory sequences. In the case of var genes, differential accessibility is observed at stages different from peak full-length var expression. Furthermore, we identify two ncRNAs (RUF6 1-2) flanking an active var gene that show clonally variant accessibility linked to the active state of the gene, and therefore could be involved in the mechanism of var mutually exclusive expression.

DISCUSSION
Despite the high morbidity associated with P. falciparum infection, we still lack precise knowledge of how chromatin structure and dynamics define transcriptional regulatory networks in this human malaria parasite and which are the regulatory elements involved in differentiation and responses to the external environment. In this study, we performed genome-wide mapping of chromatin accessibility using ATAC-seq during the intraerythrocytic cycle and compared these profiles in two different laboratory subclones that differ in their epigenomic and transcriptional profiles.
Our first objective was to obtain a 'proof of concept' of the utility of ATAC-seq to map open chromatin regulatory regions in malaria parasites in vivo. Our results conform to the pattern reported in previous studies in other organisms by showing that most accessible regions localize up to 2 Kb upstream of the ATG and that open chromatin at these sites is predictive of active transcription (3,85,95). The integration of the ATAC-seq data with previously published ChIP-seq data for various histone marks has also allowed us to classify these regions in sites of transcription initiation (TSSs) and enhancers. Approximately half of the THSSs identified here coincide with known TSSs (48,66). In contrast, we find limited overlap with the enhancer regions identified by others (33). This is not unexpected because the method used previously to discover enhancers relies solely on a localized enrichment of H3K4me3/H3K4me1 in in-tergenic regions, whereas we define enhancers on the basis of both characteristic histone modifications and chromatinaccessibility patterns, as recommended by ENCODE (82). At last, our analysis also led to the high-throughput discovery of novel regulatory regions not only associated to IDC stage transitions but also linked to the clonally variant expression of virulence genes involved in host-parasite interactions.
Accessibility at regulatory regions is modulated by chromatin remodeling processes, such as histone modifications and TF binding, and it is involved in the regulation of transcriptional activity (77,78,82). The analysis of ATAC-seq data in P. falciparum reveals temporal alterations in the chromatin accessibility of the genome during development and a strong association between these changes and stagespecific gene expression. We also show that in general, the switch from the closed to the open state leads to gene activation. This is suggestive of a predominance of activators over repressors in P. falciparum. We also report a major overlap between the active regulatory sites and TF binding sites previously determined by others. For example, we show that ChIP-seq peaks of AP2-I (41), an ApiAP2 TF that has been shown to activate invasion genes in late stages, are significantly enriched in THSSs in a similar stage-specific manner (41). But how the temporality in the transcriptional regulation is achieved remains puzzling. Based on what we know from other organisms, in some cases TFs are permanently bound to their sites and poised before the target gene is actively transcribed. Another possible scenario is that binding occurs during DNA replication, and the promoters stay poised until the stage at which transcription of the gene is activated (96,97). Our results indicate that none of the two models generally apply in P. falciparum. Instead, our data shows very clear stage-specific chromatin accessibility profiles and changes in accessibility at stage transitions correlate with transcription changes. At last, among the predicted cis-regulatory sites that are responsible for inducing stage-specific gene expression, we captured the cognate motifs of several ApiAP2 TFs already known to play important roles in P. falciparum (47) and predict novel motif sequences. During the preparation of our manuscript, another study has been published that also conducted a similar ATACseq experiment in blood stage P. falciparum parasites, showing similar temporal accessibility dynamics and reporting the discovery of known and novel TF motifs (98). Collectively, this data supports the notion that the combined action of multiple TFs is most likely an important mechanism responsible for transcriptional regulation during P. falciparum blood-stage development.
In Plasmodium, as in other organisms, cell-to-cell heterogeneity in gene expression has been proposed to serve as a bet-hedging mechanism, allowing subpopulations of parasites to vary their phenotype and survive in a changing environment (6,99,100). Changes in chromatin structure involving nucleosome occupancy and histone marking, mainly H3K9ac/me3, have been associated with clonal variation in gene expression (12,13,15,28). However, the precise mechanism for the generation and maintenance of such variability remains poorly understood, as are the regulatory elements involved. To this end, in this study, we examined and compared the chromatin accessibility landscape between two subclones of the same genetic background that display clonal variation in gene expression at several virulence genes with important functions in host-parasite interactions (6). This allowed us to identify genes that show differential accessibility in the promoter between subclones and in which the switch between the open/closed state is associated to a pattern of differential gene expression (6).
Similar to what we observe for stage-variant genes, a gain in accessibility in CVGs in one subclone is generally linked to the active state, but the accessibility status changes dynamically during the IDC development. These findings suggest that the alternate use of one or more cis-regulatory sequences in a subclone-specific manner underlies transcriptional heterogeneity (101). In agreement, we report overlap between clonally variant ATAC-seq peaks and known binding sites for putative repressors like SIP2 (35), and activators like AP2-I (41). We also predict binding sites for AP2-EXP which has been proposed to function as a repressor in Plasmodium berghei (102,103). It is generally accepted that heterochromatin/euchromatin transitions are an important mechanism in determining promoter accessibility in CVGs (8,28,86). But these transitions can also be driven by a transcription-based mechanism (104). Multiple binding of sequence-specific transcription factors in repetitive regions has been shown to modulate the production of regulatory ncRNAs that control heterochromatin formation in plants, yeast and mammals (105,106).
The var and clag3 families comprise CVGs that show mutually exclusive expression and are differentially expressed between the two subclones studied here. In these two families, the silent transcriptional state has been linked to a loss of the activating marks H3K9ac and H3K4me3 in the promoter, and a gain of H3K9me3 and HP1 (28,86,87,107,108). Our data add new insights into the TF binding events linked to clonal transcriptional variation in these families. For clag3 genes, we report subclone-specific opening of the promoter region in late-stage parasites that coincides with the active expression of the corresponding gene at this stage of development (41,62). Regarding var genes, we observed several accessible peaks in the upstream region, the exon and the intron of the single var gene expressed in the 1.2B subclone (6), but this pattern seems uncorrelated with the active state since it is similar in the 10G subclone (in which the gene is silenced) and resembles the profile observed at several other silent var genes. In contrast, we observed clonally variant accessibility of two RUF6 genes flanking the active var in 1.2B, such that they are accessible only when this var gene is active. Furthermore, other members of this gene family that are positioned adjacent to silent var genes in various chromosomes are not accessible by ATAC-seq. RUF6 genes encode GC-rich long ncRNAs have been proposed to be responsible for the var counting mechanism that drives mutually exclusive expression. In particular, it has been shown that these GC-rich elements are located at the perinuclear expression site of central and subtelomeric var genes, and that over-expression of RUF6 ncRNAs can activate var genes in trans and disrupts the monoallelic var expression pattern (89). Although with the data available this is still speculative, our data suggest the intriguing possibility that the binding of still unidentified factors to RUF6 genes may be a key event for var activa-tion and play a role in controlling mutually exclusive expression. Future research should identify the factors involved and reconcile the reported activity in trans (89) with our observation that only the RUF6 genes flanking the active var gene are accessible.
In summary, open chromatin profiling by ATAC-seq represents a novel strategy to study accessibility dynamics and transcriptional regulation in malaria parasites. The application of this technique in P. falciparum allows us to identify the regulatory elements that are likely responsible for the temporal regulation of transcription in vivo. Furthermore, our ATAC-seq data contribute to close the loop for the association between TF binding, chromatin and transcriptional states, providing new insights into mechanisms of CVG regulation and mutually exclusive expression. The challenge ahead is to apply this approach to other parasite life-stages: mosquito, liver and sexual stages, and identify ways in which these regulatory landscapes can be manipulated for malaria eradication.

DATA AVAILABILITY
ATAC-seq data are deposited in the GEO database under accession number GSE109599.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.