Identifying transcript 5′ capped ends in Plasmodium falciparum

Ethics statement

Human erythrocytes and serum were obtained from donors after providing informed written consent, following a protocol approved by the Ethics Committee, National Science and Technology Development Agency, Pathum Thani, Thailand, document no. 0021/2560.

Parasite culture and mRNA enrichment

Plasmodium falciparum strain K1 (NCBI Taxonomy ID: 5839) was cultured in vitro in human O+ erythrocytes and medium containing pooled human serum and RPMI 1640 as described previously (Shaw et al., 2007). Cultured parasites were synchronized to ring stage by Percoll gradient enrichment of mature stages followed by sorbitol treatment (Lambros & Vanderberg, 1979). Parasites were harvested immediately (ring stage), 12 h (trophozoite stage) and 24 h (schizont stage) after sorbitol treatment. Parasites were liberated from the host cell by treatment with 0.1% (w/v) saponin. Parasite total RNA was obtained using Trizol reagent according to the manufacturer’s instructions (Invitrogen). Purified total RNA was stored in 75% ethanol at −80 °C before use. On the day of library preparation, total RNA was dried and then resuspended in nuclease-free water. Total RNA concentration was estimated by Nanodrop ND1000 measurement assuming A260 = 1.0 is equivalent to 40 µg/mL RNA. Up to 100 µg of total RNA was used for mRNA enrichment using a Dynabeads oligo dT25 mRNA kit (Thermo Scientific) or mRNA magnetic beads kit (New England Biolabs) followed by treatment with 1 unit of XRN-1 nuclease (New England Biolabs) for 30 min at 37 °C. XRN-1 enriched mRNA was purified by acidic phenol:chloroform extraction and ethanol precipitation. The mRNA was redissolved in 20 µL of nuclease-free water and genomic DNA was removed using a Turbo DNA-free kit following the manufacturer’s instructions (Ambion).

Synthesis of cDNA for HiSeq libraries

In vitro synthesized spike-in RNA SIRV-Set 3 (Iso Mix E0/ERCC, Lexogen) was 5′ capped with Vaccinia capping system (New England Biolabs) following the manufacturer’s instructions. A 100–200 ng sample of parasite mRNA from ring or trophozoite stage (estimated by Nanodrop measurement) was mixed with 1.5 ng of 5′ capped SIRV-Set 3 RNA. Approximately the same amounts of parasite mRNA from schizont stages were used for reverse-transcription without the addition of spike-in RNA. The mRNA samples were incubated at 94 °C for 3 min in 1x first-strand cDNA synthesis buffer supplied with SuperscriptIV reverse transcriptase enzyme (Invitrogen) to fragment the RNA. Fragmented mRNA was chilled on ice for 5 min before reverse-transcription with 50 pmol of 5′ tagged random primer (RTNGS1 (Table S1) for adapter ligation method1, RTCIRC (Table S1) for adapter ligation method2; see below) with SuperscriptIV as described by the manufacturer in a final volume of 20 µL. After the addition of SuperscriptIV enzyme, reactions were incubated at 25 °C for 5 min, followed by heating to 37 °C for 30 min, 50 °C for 30 min, and then 70 °C for 15 min. First-strand cDNA was diluted to 100 µL with 10 mM Tris 1 mM EDTA and digested with 1 µL of RNaseA/T1 mix (Thermo Scientific) for 10 min at 37 °C to remove incompletely reverse-transcribed RNA. The mRNA/cDNA hybrid was purified using an equal volume of Agencourt^® AMPure^® XP beads (Beckman Coulter) and recovered from beads in 100 µL nuclease-free water. A 20 µL sample of purified mRNA/cDNA was kept for processing as the unenriched control sample.

Synthesis of cDNA for MiSeq libraries

Samples of approximately 1 µg of mRNA (estimated by Nanodrop measurement) were incubated at 94 °C for 5 min in 1x first-strand cDNA synthesis buffer supplied with SuperscriptIII reverse transcriptase enzyme (Invitrogen) to fragment RNA. Fragmented mRNA was chilled on ice before reverse-transcription with 50 pmol of 5′ tagged random primer (RTNGS1 for method1, RTCIRC for method2) with SuperscriptIII as described by the manufacturer in a final volume of 20 µL. Method1 reverse-transcription reactions also contained 40 µM biotin-11-dUTP (Thermo Scientific) to label first-strand cDNA with biotin. After the addition of SuperscriptIII enzyme, reactions were incubated at 25 °C for 5 min, followed by heating to 37 °C for 30 min, and then 50 °C for 30 min. First-strand cDNA was diluted to 100 µL with 10 mM Tris 1 mM EDTA and digested with 1 µL of RNaseA/T1 mix (Thermo Scientific) for 10 min at 25 °C to remove incompletely reverse-transcribed RNA. The mRNA/cDNA hybrid was purified by proteinase K treatment, phenol:chloroform extraction and desalting on a Amicon-30 column (Ambion). The volume of desalted mRNA/cDNA hybrid was adjusted to 100 µL with nuclease-free water. A 20 µL sample of desalted mRNA/cDNA was kept for processing as the unenriched control sample

Enrichment of full -length cDNA

P. falciparum eIF4E protein fused to glutathione-S transferase (GST-PfeIF4E) (Shaw et al., 2007) was used for enriching full-length cDNA. GST-PfeIF4E was expressed as recombinant protein in Escherichia coli and purified by SP-sepharose cation exchange chromatography as described previously (Shaw et al., 2007). Purified protein was buffer-exchanged using desalting columns and stored in a solution of 15% glycerol/RNase-free phosphate-buffered saline (PBS, Ambion). 5′ cap-binding activity of GST-PfeIF4E protein was determined by m⁷GTP pulldown assay as described previously (Shaw et al., 2007), except that γ-aminophenyl-m⁷GTP (C10-spacer)-agarose (Jena Bioscience) was used as the affinity support. For each full-length cDNA enrichment, approximately 100 µg of SP-sepharose purified GST-PfeIF4E protein was immobilized on 50 µL of glutathione magnetic beads (Thermo Scientific) in PBS. The glutathione magnetic beads were washed four times with 0.5 mL of PBS to remove unbound protein and then resuspended in 100 µL of PBS. The remainder of the purified mRNA/cDNA sample was incubated with bead-bound GST-PfeIF4E protein for 20 min at 25 °C with agitation (750 rpm Thermomixer). The beads were then washed thrice with 0.5 mL PBS and the mRNA/cDNA eluted in 100 µL of 1% sodium dodecyl sulfate/0.2 M sodium chloride. The eluted mRNA/cDNA was purified with an equal volume of Agencourt^® AMPure^® XP beads (Beckman Coulter) and recovered from beads in 20 µL of nuclease-free water. RNA was removed from enriched and unenriched cDNA samples by alkaline hydrolysis (15 min treatment at 65 °C with 0.2 N NaOH). Alkaline-treated cDNA was neutralized with an equal volume of 1 M HEPES pH 7.4 solution. First-strand cDNA was purified with an equal volume of Agencourt^® AMPure^® XP beads (Beckman Coulter) and recovered from beads in 20 µL of nuclease-free water.

DNA adapter ligation to cDNA

In some experiments, sequencing adapter was added by cDNA tailing followed by double-stranded adapter ligation (hereafter referred to as method1). A ribo-G tail was added to the purified single-stranded cDNA primed with RTNGS1 (Table S1) using terminal transferase (TdT, New England Biolabs). The tail length is limited to four Gs (Schmidt & Mueller, 1996). The TdT tailing reactions contained 2 mM GTP, 0.25 mM CoCl₂, 1 unit of TdT enzyme in 1x TdT enzyme buffer and first-strand cDNA. TdT reactions were performed in 20 µL reaction volumes for 20 min at 37 °C. Ribo-tailed cDNA was purified with an equal volume of Agencourt^® AMPure^® XP beads (Beckman Coulter) and recovered from beads in 20 µL of nuclease-free water.

Double-stranded DNA adapter was made by combining 4 nmol each of NGS1 and NGS1COMP oligonucleotides (Table S1) in a volume of 100 µL with 10 mM NaCl and 10 mM Tris pH 8.0. The NGS1 oligonucleotide is 5′-phosphorylated to act as a donor in ligation and blocked with a three-carbon spacer at the 3′ end to prevent concatemerization of adapters. Annealing was accomplished by heating the mixture for 3 min at 80 °C followed by slow cooling to 25 °C (−0.1 °C/min). Double-stranded DNA adapter was ligated to first-strand cDNA using T4 RNA ligase 2 (RNL2, New England Biolabs). The ligation reactions contained 40 pmol of DNA adapter, 7.5% (w/v) PEG₆₀₀₀, 1x RNL2 enzyme buffer, 2.5 units of RNL2 enzyme and first-strand cDNA in a volume of 20 µL. The ligation reactions were performed for 99 cycles of 37 °C for 30 s, 22 °C for 30 s. After the completion of the RNL2 adapter ligation reactions, cDNA for HiSeq sequencing was diluted to 100 µL with nuclease-free water and purified with an equal volume of Agencourt^® AMPure^® XP beads (Beckman Coulter). Purified cDNA was recovered in 50 µL of nuclease-free water. For adapter-ligated cDNA samples sequenced on the MiSeq platform, ligation products were purified using 25 µL Streptavidin M280 magnetic beads following the manufacturer’s recommendations (Invitrogen). Synthesis of second-strand cDNA was performed on the beads using T7 DNA polymerase following the manufacturer’s recommendations (New England Biolabs). The beads were resuspended in 50 µL of T7 reaction mix (0.3 mM dNTPs, 2.5 units T7 DNA polymerase, 1x T7 DNA polymerase buffer). The reaction was performed for 15 min at 37 °C. The reaction was terminated by washing the beads in 0.4 mL of washing solution (40 mM Tris pH 8.0, 10 mM MgCl₂). The second-strand cDNA was eluted from the beads by heating for 3 min at 99 °C in a suspension of 50 µL 1x Sodium Saline Citrate (150 mM sodium chloride, 15 mM trisodium citrate pH 7.0).

An alternative adapter ligation strategy (hereafter referred to as method2) was also employed. Purified first-strand cDNA primed with RTCIRC in which RNA had been removed by alkaline treatment was circularized with 100 U of CircLigaseII enzyme, 2.5 mM MnCl₂, 1 M betaine, and 1x CircLigaseII buffer in a reaction volume of 20 µL as recommended by the manufacturer (Epicentre). The reaction was incubated at 60 °C for 1 h and the reaction terminated by heating to 80 °C for 10 min. The reaction was diluted to 100 µL with nuclease-free water and cDNA purified using half the volume of Agencourt^® AMPure^® XP beads (Beckman Coulter). Purified cDNA was recovered in 50 µL of nuclease-free water.

Library purification and Illumina sequencing

Adapter-ligated cDNA was used as a template for PCR amplification to make sequencing library. For samples sequenced using the HiSeq platform, PCRs contained adapter-ligated cDNA template (1 µL), 2.5 pmol each of PE1 and ScriptSeq primers (Table S1), 200 µM dNTPs, 0.625 U PrimeSTAR^® GXL DNA polymerase (Takara) and 1x buffer supplied with the enzyme in a reaction volume of 25 µL. The PCR program used was 98 °C for 30 s, followed by 18–25 cycles of 98 °C for 10 s, 68 °C for 60 s and a final extension of 68 °C for 5 min. For samples sequenced using the MiSeq platform, PCRs contained adapter-ligated cDNA template (1 µL), 12.5 pmol each of PE1 and ScriptSeq primers, 200 µM dNTPs, 2 mM MgCl₂, 0.5 units of Platinum Taq enzyme (Invitrogen) and 1x Platinum Taq buffer in a reaction volume of 50 µL. The PCR program used was 94 °C for 90 s, followed by 18–25 cycles of 94 °C for 10 s, 60 ° C for 60 s. The optimal number of PCR cycles for each sample was determined empirically by visualization of products separated by 1.5% agarose gel electrophoresis from pilot reactions conducted with varying numbers of cycles. Amplified DNA fragments 400–600 bp in size were excised from the gel and purified using a MinElute gel extraction kit (Qiagen), with the modification that agarose was solubilized in QG buffer at room temperature to prevent DNA denaturation. DNA was quantified using a Qubit™ dsDNA HS Assay kit (Invitrogen). Libraries were pooled in equimolar ratios according to the ScriptSeq index primer recommendations (Illumina).

Pooled libraries were submitted to Novogene AIT (Singapore) for Illumina standard protocol sequencing on one lane of a HiSeq2000 flow-cell (Illumina), or to the Faculty of Medicine, Chulalongkorn University, Bangkok, Thailand for MiSeq sequencing. For MiSeq sequencing, pooled libraries were processed with 1% Illumina phiX174 control and loaded onto MiSeq v3 flow cells at 10 pM. 150 bp paired-end sequencing was performed using the standard sequencing primers as recommended by the manufacturer (Illumina). For MiSeq sequencing of libraries prepared using adapter ligation method 1, the sequencing “recipe” for the MiSeq instrument was modified to perform the first four sequencing cycles as “dark”, meaning no data are recorded for the purpose of identifying sequencing clusters. This modification was necessary since the homopolymer tail added by TdT during library construction provides insufficient diversity for generating a high density of clusters. Without this modification, there will be data loss in the cluster identification step (clusters passing filter quality control) of Illumina sequencing.

Analysis of transcriptomic data from 5′ capped nucleotide enrichment experiments

The 5′ end of read1 from cDNA adapter ligation method1 data obtained using the HiSeq platform was trimmed of homopolymer sequence added by TdT using Cutadapt 1.18 (Martin, 2011). Read trimming was not performed for method 1 data obtained using the MiSeq platform, as signals for bases 1–4 were not recorded for this platform. Adapter sequence in the 3′ end of reads was trimmed and low-quality reads were removed using fastp with default settings (Chen et al., 2018). Preprocessed paired read data were aligned to the combined P. falciparum 3D7 v3.2 (Böhme et al., 2019) / SIRVomeERCCome (Lexogen) reference genome using HISAT2 (Kim et al., 2019), in the spliced mode guided by the P. falciparum / SIRVomeERCCome genome annotation with maximum intron length limited to 5,000 bp and other settings as default. Alignment summary statistics are shown in Table S2. To assess biological reproducibility, gene expression was determined for annotated P. falciparum genes using unenriched .bam files from each experiment as input to StringTie (Pertea et al., 2016). Gene expression values reported as transcripts per million were log₁₀ transformed and used for sample pairwise Pearson correlations using the ggcorrplot package in R (Kassambra, 2019).

For analysis of 5′ cap nucleotide enrichment, .bam alignment files were used as input for the ToNER program (Promworn et al., 2017) with default settings, i.e., depth value calculated from the read start position, pseudo-read added to depth values at all transcribed nucleotides and no filtering of positions for statistical modeling. This program calculates enrichment scores from paired feature-enriched and unenriched transcriptomic data. The enrichment scores are then transformed and fitted to a normal distribution by the Box–Cox procedure, and statistics of enrichment for each nucleotide reported. To increase statistical power, combined P-values from independent enrichment experiments were calculated by Fisher’s method. The enrichment signals at individual nucleotides were clustered into genomic intervals representing the 5′ capped ends of major transcripts using the CAGEr package (Haberle et al., 2015) as follows. Clustering of nucleotide signals in CAGEr requires integer counts, which are referred to in the CAGEr manual as tags. To generate tags for clustering, the reciprocal of enrichment P-value was rounded to the nearest integer. To militate against program abort owing to the constraint of maximum allowed integers in R, the maximum tag count was set as 1e7. Tags were clustered with the distclus function in CAGEr with no normalization and maximum distance between merged clusters of 20 nt. The dominant_ctss position in each cluster reported by CAGEr was annotated as a 5′ capped nucleotide. CAGEr clusters overlapping known 5′ ends in SIRV spike-in RNAs (assigned as true positives) were identified with BEDTools (Quinlan & Hall, 2010). The sum of tags (tpm) in each CAGEr cluster mapping to SIRVomeERCCome reference was used for Receiver Operator Characteristic (ROC) analysis with the plotROC package (Sachs, 2017). The optimal tpm value for filtering CAGEr clusters was identified using the OptimalCutpoints R package (López-Ratón et al., 2014) with default settings. For quantitative analysis of External RNA Controls Consortium (ERCC) RNA spike-ins, the CAGEr clusters reported from tag input of combined P-values of enrichment from four experiments were used. An RNA 5′ end was called as detected if a CAGEr cluster overlapped its known location. The empirical cumulative distribution functions of detected and undetected 5′ ends with respect to natural log of ERCC RNA concentration reported by the manufacturer (Lexogen) were determined using the base R ecdf function. The distributions of detected and undetected 5′ end groups were compared by two-sample Kolmogorov–Smirnov test in R. The natural log tpm values for CAGEr clusters overlapping detected 5′ ends were plotted against the natural log concentration of each RNA using the ggplot2 package (Wickham, 2009) in R.

Analysis of genomic context in the vicinity of 5′ capped nucleotides

5′ capped nucleotide signals detected as CAGEr clusters above threshold (tpm >32) in separate stages of P. falciparum development (including cluster signals within 100 bp of each other) were merged into 17,961 non-overlapping genomic intervals using BEDTools (Quinlan & Hall, 2010). Within each interval, the dominant_ctss nucleotide with highest integer count (tpm.dominant_ctss reported by CAGEr) across the sampled developmental stages was annotated as the representative 5′ capped nucleotide for genomic analysis. Genomic sequences flanking 100 bp on either side of the 17,961 5′ capped nucleotide (test set) and 17,961 random positions selected using BEDTools (Quinlan & Hall, 2010) (control set) were obtained from the P. falciparum 3D7 v3.2 reference genome (Böhme et al., 2019) using the seqPattern package in R (Haberle, 2020). Sequence analysis was performed using the kpLogo program (Wu & Bartel, 2017) using the control set sequences as background and other settings as default. For analysis of epigenetic features, publicly available datasets from published studies were processed to construct genome coverage files (bigwig format) as described below in sub-section analysis of external epigenetic data. Plots of average coverage for each feature were made using the genomation package (Akalin et al., 2015) in R. To mitigate the effect of extreme values, the top and bottom 5% of scores were clipped using the winsorize function in the genomation package.

Mismatched base analysis of -1 positions

The first base of all uniquely aligned reads (with secondary alignments removed from alignment .bam files using SAMtools) was extracted from all datasets. For reads mapping to the reference (–) strand, the complement of the first base was taken. Bases mismatched to the reference sequence (corresponding to soft-clipped bases in the alignment files) were counted for two groups. The first group comprised reference positions one base upstream of dominant 5′ capped nucleotides in non-overlapping genomic intervals (17,961 positions), and the second group all other reference positions. Contingency tables were constructed for each dataset for the count of the tested mismatched bases and all other mismatched bases for both groups of nucleotides (Table S3). One-tailed Fisher’s exact tests were performed of the alternative hypothesis that the true odds ratio is greater than 1. One-tailed tests were performed because mismatched bases under-represented in the first position of reads are not of interest. Bonferroni-corrected P-values less than 0.001 were considered significant. For analysis of reference sequence patterns at -1 positions with high counts of mismatched reads, positions with 10 or more aligned reads and >50% of reads with mismatched first base in any dataset were selected.

Unsupervised cluster analysis of 5′ capped nucleotides

17,961 5′ capped nucleotides in non-overlapping genomic intervals annotated from our data (see above) were clustered using epigenetic data, including occupancies of H2A.Z, H3K9 acetylation and H3K4 trimethylation (Bártfai et al., 2010) and H2B.Z (Hoeijmakers et al., 2013). For details of how chromatin mark occupancies were determined, see below in sub-section analysis of external epigenetic data. Principal Component Analysis (PCA) scores were calculated using the MacroPCA package in R (Hubert, Rousseeuw & VandenBossche, 2019). The data were pre-processed by transformation with natural logarithms, scaling by unit variance, means-centering and removal of rows with too many missing values. The first three principal components explained 83.3% of the variance; hence, the PCA scores from the first three principal components were used for clustering. To determine if more than one cluster existed in the data, statistical tests of unimodality were performed on PC1 scores using the multimode package in R (Ameijeiras-Alonso, Crujeiras & Rodríguez-Casal, 2019). Test P-values less than 0.05 were considered significant. Clustering was performed using the cross-entropy clustering (CEC) package in R (Tabor & Spurek, 2014; Spurek et al., 2017). The algorithm employed in this program is a hybrid of k-means and Gaussian mixed model, and can thus separate clusters with a variety of shapes. The settings used for CEC clustering were: Gaussian distribution models unconstrained, nstart =1000, initial clusters =5 and card.min =20%. To determine independently how many relevant clusters are present, the PCA scores from PC1, PC2 and PC3 were analysed using the NbClust package in R (Charrad et al., 2014). Cluster indices were determined using the Euclidean distance matrix and k-means method for n = 2 to n = 6 clusters. All indices were calculated except Gplus and Tau, which could not be determined owing to computational constraint, and Gamma which is only applicable for hierarchical clustering.

Analysis of external transcriptomic data

To assess the reproducibility of 5′ capped nucleotides in P. falciparum, data were analyzed from independent published studies using different methods for 5′ capped nucleotide enrichment. P. falciparum 5′ cap sequencing data reported in Adjalley et al. (2016) were downloaded from the NCBI GEO database series under accession number GSE68982 using sratoolkit.2.10.0-ubuntu64 (SRA Toolkit Development Team, 2014). The read1 fastq files were pre-processed with UMI-tools (Smith, Heger & Sudbery, 2017) to move the 8 bp unique molecular index to the read header. Pre-processed paired data files were filtered and trimmed of 3′ adapter sequence using fastp (Chen et al., 2018) with default settings. Data from P. falciparum SMART-enrichment experiments reported in Kensche et al. (2016), European Nucleotide Archive accession numbers ERR861692 and ERR861693, were uploaded to the public server at usegalaxy.org. Read1 were filtered and trimmed using the Barcode Splitter tool to retain reads that started with smartseq2 adapter (AAGCAGTGGTATCAACGCAGAGTACATGGG), allowing up to three mismatches or indels. The filtered read1 and read2 from 5′ cap sequencing and SMART-enrichment experiments were processed using fastp (Chen et al., 2018) with default settings. Processed paired 5′ cap sequencing and SMART-enrichment data from fastp were aligned to the P. falciparum 3D7 v3.2 reference genome (Böhme et al., 2019) using HISAT2 (Kim et al., 2019) in the spliced mode guided by the genome annotation (PlasmoDB release 44) with maximum intron length limited to 5,000 bp and other settings as default. PCR duplicates in 5′ cap sequencing data were removed using UMI-tools (Smith, Heger & Sudbery, 2017) with the–unmapped-reads =use option. Unmapped and read2 reads were removed from .bam files and a merged .bam file was created from all experiments in each study using SAMtools (Li et al., 2009).

The filtered .bam files from individual experiments were used to obtain read depth at each position in the genome with BEDTools (Quinlan & Hall, 2010), with -strand and -5 options. The read depth values were used as ctss input for CAGEr (Haberle et al., 2015) for detection of 5′ capped nucleotide signals. CAGEr was run with the distclus function, no normalization and maximum distance between merged clusters of 20 nt. The CAGEr analysis results of 5UTR-seq experiments in P. falciparum were obtained from Table S9 reported in a previous study (Chappell et al., 2020). A merged .bam file of aligned data from all 5UTR-seq experiments was provided by Dr. Lia Chappell.

The dominant_ctss nucleotide genomic locations and the corresponding tpm.dominant _ctss values reported for CAGEr clusters in each dataset were concatenated and used to construct .bed and .bedgraph files of 5′ capped nucleotides for each study. In addition to CAGEr clusters, the read depth at each position in the genome was obtained from the merged .bam file from all experiments in each study using BEDTools (Quinlan & Hall, 2010) with -strand and -5 options. To assess the agreement of 5′ cap nucleotide assignment in our data with other studies, agreement was scored if the depth value from the combined .bam file of all experiments in the same study was two or greater.

Direct RNA sequencing data from P. falciparum obtained using the Oxford Nanopore MinION platform reported in Lee et al. (2021) under SRA accession number SRR11094274 were uploaded to the public server at usegalaxy.org and aligned to the P. falciparum 3D7 v3.2 reference genome (Böhme et al., 2019) using Minimap2 (Li, 2018) with settings long-read spliced alignment, maximum intron 5 kb, and search GT-AG on the transcript strand only. The .bam alignment file was used as input to the Full-Length Alternative Isoform analysis of RNA (FLAIR) program (Tang et al., 2020a) with the–nvrna option using reference genome sequence and annotation (PlasmoDB release 44) for correction of isoform sequences. The corrected and inconsistent isoform sequences outputted by FLAIR were concatenated to create a single transcript annotation .gtf file. To assess whether transcripts annotated by FLAIR were 5′ capped, the transcript 5′ end locations were compared with 5′capped nucleotides annotated from 5′ cap enrichment data using BEDTools (Quinlan & Hall, 2010). Genomic locations in .bed file format were constructed for each 5′ cap enrichment study from the dominant_ctss nucleotide reported in all CAGEr clusters. FLAIR transcripts were considered full-length (5′ capped) if the transcript 5′ end was 20 nt or closer to a dominant_ctss nucleotide.

Analysis of external epigenetic data

Epigenetic data for analysis of 5′ capped nucleotide genomic context were obtained from published studies, including chromatin immunoprecipitation sequencing (ChIP-seq) data reported in Bártfai et al. (2010), Hoeijmakers et al. (2013), Karmodiya et al. (2015), Lu et al. (2017), Tang et al. (2020b), Bhowmick et al. (2020) and micrococcal nuclease-digested chromatin sequencing (MNase-seq) reported in Kensche et al. (2016) were obtained from the NCBI database under series accession numbers GSE23867, GSE39702, GSE63369, GSE85478, PRJNA612099, GSE142803 and GSE66185, respectively. The .csfasta files in accession GSE63369 were converted to .fastq using Cutadapt 1.18 (Martin, 2011) and aligned to the P. falciparum 3D7 v3.2 reference genome (Böhme et al., 2019) using bowtie−1.2.3-linux-x86_64 (Langmead et al., 2009) with the following options: -C -v 2 –best –strata m 3. The data from the other studies were uploaded to the public server at usegalaxy.org and processed using fastp (Chen et al., 2018) with default settings. For MNase-seq data, read length was trimmed with fastp to a maximum of 72 bp. Processed read data were aligned to the P. falciparum 3D7 v3.2 reference genome (Böhme et al., 2019) using bowtie2 (Langmead & Salzberg, 2012) with default settings. The .bam files for experimental replicates or multiple runs of the same library were merged using SAMtools (Li et al., 2009).

In order to calculate nucleosome occupancy from MNase-seq data at each sampled timepoint of development, sonicated genomic DNA control datasets with matching insert size distributions were created. The count of fragments with the same insert size in the alignment file generated from each MNase-seq experiment was obtained using SAMtools (Li et al., 2009). Next, insert size distributions were determined from the count of inserts in 31 bins of varying insert size (32–52, 53–72, 73–82, 83–92, 93–102, 103–12, 113–122, 123–132, 133–142, 143–152, 153–162, 163–172, 173–182, 183–192, 193–202, 203–212, 213–222, 223–232, 233–252, 253–272, 273–292, 293–312, 313–332, 333–352, 353–372, 373–392, 393–412, 413–432, 433–452, 453–472, and 473–500 bp). Alignment .bam files of sonicated genomic DNA control with defined insert sizes were made using BamTools (Barnett et al., 2011). The alignment files of sonicated genomic DNA control with defined insert sizes were randomly sampled to the equivalent depth in the MNase-seq experimental dataset and merged to create a genomic DNA control with matching insert size distribution using SAMtools (Li et al., 2009).

Normalization factors for each sample dataset from all epigenetic experiments were calculated by dividing the reference genome size (23,332,839 bp) by the total depth of coverage in .bam files obtained using SAMtools (Li et al., 2009). The read depth at each position in the genome was obtained using BEDTools (Quinlan & Hall, 2010). Nucleotide occupancy for each epigenetic feature of interest was calculated as the ratio of normalized read depth in each test dataset to developmental stage-matched control. To prevent division by zero, a pseudo depth value equal to 0.001 was added to every position for test and control datasets. For ChIP-seq data reported in Lu et al. (2017), RNA polymerase II (RNApolII) occupancy was calculated by subtraction of timepoint-matched normalized IgG negative control read depth from normalized ChIP-seq read depth. Nucleotide occupancy values were used to create genome coverage files in .bigwig format using the bedGraphToBigWig tool (Kent et al., 2010).

Genomic intervals annotated as P. falciparum accessible chromatin from Assay for Transposase-Accessible Chromatin (ATAC-seq) experiments were obtained from the supplementary information files reported in Toenhake et al. (2018), Ruiz et al. (2018) and Yin et al. (2020). The accessible chromatin intervals reported in each study were intersected using BEDTools (Quinlan & Hall, 2010).

Data availability

Custom scripts written for analysis of transcriptomic data are available from the GitHub repository: https://github.com/BSI3/5CAPture-seq. Sequencing data generated in this study are available from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database under series accession number GSE103036.

Outline of the method for identifying 5′ capped nucleotides

We present a new method for identifying 5′ capped nucleotides from mRNA transcripts. mRNA is purified from the biological sample and used to synthesize cDNA. A portion of the cDNA is reserved as unenriched for generation of a control library and the remainder is enriched for full-length cDNA (Fig. 1A). Full-length cDNA fragments are enriched using recombinant eIF4E protein, which is capable of binding to the 5′ cap nucleotide in first-strand cDNA products (Edery et al., 1995). We enriched full-length cDNA using recombinant P. falciparum eIF4E protein, which has an affinity for the 5′ cap nucleotide similar to mammalian eIF4E (Shaw et al., 2007). To mitigate 5′ end biases when adding adapter, we employ two different adapter ligation protocols, namely cDNA tailing followed by double-stranded adapter ligation (adapter ligation method1) and intramolecular ligation (adapter ligation method2) (Figs. 1B, 1C). 5′ end bias patterns were observed among the unaligned reads, in which C was markedly lower for the first base of method1 HiSeq libraries compared with downstream bases (Fig. S1), and G markedly lower for the first base of all method2 libraries (Figs. S1 and S2).

Accuracy and sensitivity of the method

To assess the accuracy and sensitivity of the proposed method, we performed four experiments with P. falciparum mRNA spiked with 5′ capped synthetic RNAs, including two replicates for each adapter ligation protocol. The reads mapped to the in silico SIRVomeERCCome reference sequence were used to determine enrichment signals for synthetic RNAs with known 5′ ends (Fig. 2A). The performance for detecting 5′ capped ends was assessed by Receiver Operator Characteristic (ROC) plot (Fig. 2B). As the same batch of synthetic RNA was used for the four spike-in experiments, we also tested whether combining P-values from replicate experiments could improve performance. Individual experiments with method2 gave slightly higher performance than method1 as shown by the greater area under ROC curves. The greatest performance was achieved by combining P-values from all four experiments. From the ROC curve of all four experiments combined, the optimal cutpoint for classifying 5′ capped nucleotide signals was determined (Fig. 2B). We assessed the sensitivity of the method for detecting 5′ capped ends in absolute terms by comparing the distributions of the ERCC spike-in RNA concentrations for detected versus undetected 5′ ends, which were significantly different (Fig. 2C). Moreover, the 5′ capped nucleotide signals for ERCC spike-in RNAs showed a significant positive correlation with the known RNA concentrations indicating that the enrichment signals are quantitative and reflect RNA abundance (Fig. 2D). From the analysis of spike-in RNA enrichment data, we determined the performance of the method using different adapter-ligation protocols and showed that combining P-values from all available replicate experiments gives the best performance.

Annotation of 5′ capped nucleotides in P. falciparum

We performed 12 experiments in total for identifying 5′ capped ends in P. falciparum mRNA isolated from different stages of the IDC. Six experiments were conducted with high sequence coverage on the HiSeq platform (including the four with synthetic RNA spike-in) and six at lower coverage on the MiSeq platform. To assess variability among samples, gene expression profiles were compared. All pairwise combinations of samples from the same stage of the IDC (including samples prepared using different library and sequencing protocols) showed Pearson r values greater than 0.5, whereas some pairwise combinations of samples from different stages showed Pearson r less than 0.5 (Fig. 3A). Enrichment P-values were thus combined from four replicate experiments of each stage of development to increase the power of detection. 5′ capped nucleotide signals were annotated for ring, trophozoite and schizont stages of P. falciparum IDC based on the threshold of 5′ capped nucleotide signal determined from analysis of spike-in RNA (Table S4). The 5′ capped nucleotide signals detected for P. falciparum encompass genomic regions of variable width (Table S4). In other transcriptomic studies of eukaryotic 5′ capped nucleotides, the majority of signals correspond with core promoters, or genomic regions of varying width in which multiple transcription start sites (TSS) are used at different frequencies (Zhao et al., 2011). The interquantile signal width reported by CAGEr can be used for the analysis of 5′ capped nucleotide distribution (Haberle et al., 2015). Bimodal distributions comprising sharp (<10 bp) and broad (>10 bp) peaks were observed for each sampled stage of the P. falciparum IDC; however, the majority of signals showed broad peak widths (Fig. 3B). Many of the 5′ capped nucleotide signals from the three stages of development overlapped the same genomic regions, which made it difficult to separate signals of alternative 5′ ends used dynamically throughout the IDC. Therefore, to simplify the genomic analysis we collapsed the signals across different stages into 17,961 non-overlapping intervals. The dominant 5′ capped nucleotide (5CN) with the strongest signal across the IDC in each genomic interval was investigated in more detail.

Sequence patterns in the vicinity of dominant 5′ capped nucleotides

Eukaryotic promoter regions are characterized by A/T richness and the presence of position-specific short motifs such as upstream TATA boxes and pyrimidine/purine at the −1/ + 1 position (where +1 denotes TSS) (Müller & Tora, 2014). To test whether TSS sequence patterns existed at 5CN, base compositions were investigated. The average base composition at 17,961 5CN and 100 flanking genomic bases showed elevated A composition, in particular at +1 and positions immediately downstream compared with random (Fig. 4A).

Aligned reads may possess non-reference (mismatched) 5′ terminal nucleotides that can affect the accuracy of annotating transcript 5′ ends. For example, analysis of CAGE data requires correction for the non-reference G present on the first base of the majority of aligned reads, which originates from the cap signature (Haberle et al., 2015). The odds ratio of mismatched base counts from read1 starts aligned one base upstream of 5CN (−1 position) was compared with all other transcribed positions. Odds ratios significantly greater than 1 were observed for mismatched A(8/12 datasets), C(2/12 datasets), G(7/12 datasets) and T(1/12 datasets) (Table S3). Although mismatched bases are more common at −1 than at other transcribed positions, it should be noted that the majority of aligned reads do not start with a mismatch (Table S3). To test whether base mismatching of aligned reads could be influenced by reference sequence, sequence patterns were investigated among genomic regions flanking 5CN with high frequencies (>50%) of aligned reads with mismatched first bases (Fig. 4B). Among genomic regions in the vicinity of all 5CN, we observed significant over-representation of A-rich motifs at the 5CN (+1) and downstream positions. CA was the most significant over-represented motif at the −1 position. Different over-represented motifs were identified among genomic regions with high frequencies of aligned reads with mismatches at the -1 position. T, AA, TA and CA were identified as the most significant over-represented motifs among sequences with high frequencies of non-reference A, C, G and T at -1 positions, respectively.

Genomic features in the vicinity of dominant 5′ capped nucleotides

9,334 (52%) of the 5CN are located within P. falciparum protein-coding exons (Table S4). Given that eukaryotic TSS are typically located outside of protein-coding regions, many of the exonic 5CN may represent transcriptional noise unrelated to TSS. We hypothesized that epigenetic information could be used to classify TSS in a manner naïve to gene annotation. Eukaryotic core promoter regions possess a distinctive chromatin architecture, characterized by high occupancy of nucleosomes with variant histone H2AZ, and histone modifications including H3K9 acetylation and H3K4 trimethylation (Jiang & Pugh, 2009; Müller & Tora, 2014). We used published P. falciparum epigenetic data of H2A.Z, H3K9 acetylation and H3K4 trimethylation (Bártfai et al., 2010) and H2B.Z (Hoeijmakers et al., 2013) for unsupervised clustering of 5CN. H2B.Z is an apicomplexan-specific histone variant that colocalizes with H2A.Z (Hoeijmakers et al., 2013; Petter et al., 2013). The P-values from unimodal tests were below the threshold of significance (Silverman critical bandwidth test P = 0.01; Hall and York critical bandwidth test P < 2.2e⁻¹⁶; Ameijeiras-Alonso et al. excess mass test P < 2.2e⁻¹⁶; Cheng and Hall excess mass test P < 2.2e⁻¹⁶; Fisher and Marron Cramer-von Mises test P < 2.2e⁻¹⁶; Hartigan and Hartigan dip test P < 2.2e⁻¹⁶). Hence, we rejected the null hypothesis that the data are unimodal and conclude that more than one cluster exists. Two clusters of 5CN (C1, 6,450 positions; C2, 11,075 positions) were resolved by unsupervised clustering (Fig. S3; Table S4). Furthermore, the number of relevant clusters by the majority vote of clustering indices was two (Fig. S3). C1 5CN overlap genomic regions with locally high occupancies of H2AZ, H2BZ, H3K9 acetylation and H3K4 trimethylation chromatin marks across the IDC, whereas C2 5CN show low occpancies of these chromatin marks. (Fig. 5A). C1 5CN overlap high occupancies of other histone acetylation and methylation marks, with the notable exception of H3K4me1, which is elevated among C2 5CN (Fig. S4). C1 and C2 5CN are distinguished by other epigenetic features, including elevated nucleosome, PfGCN5 histone acetyltransferase and RNApolII occupancies for C1 (Fig. S5). A higher proportion of 5CN are located outside of protein-coding exons for C1 than C2 (Fig. 5B), and the average distance to the nearest accessible chromatin feature is less for C1 than C2 (Fig. 5C). Analysis of genomic sequence in the vicinity of 5CN showed different patterns of motifs among C1 and C2 (Fig. 5D). T-rich motifs are over-represented upstream of C1 5CN, whereas they are under-represented in the vicinity of C2. The most significant over-represented motifs at the −1 position are TA and CA for C1 and C2, respectively. Previous studies of transcript 5′ ends in P. falciparum highlighted patterns of genomic features similar to the C1 group, including elevated H2AZ, H3K9 acetylation and H3K4 trimethylation chromatin marks and high local A/T contents (Adjalley et al., 2016; Kensche et al., 2016; Chappell et al., 2020) Correspondingly, a greater proportion of 5CN in the C1 group than the C2 group are supported by data in the other studies (Fig. 5E; Table S4).

Patterns of 5′ capped nucleotides among genes and associated transcripts

The transcriptomic surveys from data generated in this study and others (Adjalley et al., 2016; Kensche et al., 2016; Chappell et al., 2020) provide comprehensive information of transcript 5′ capped ends. In order to annotate transcripts, data of complete (end-to-end) transcript structures are required, which can be difficult to obtain from short-read cDNA data owing to breaks in coverage at highly AT-rich or unmappable regions of the P. falciparum genome (Chappell et al., 2020). To our knowledge, the largest available dataset of complete transcript structures was obtained by Nanopore long-read direct RNA sequencing (Lee et al., 2021). In direct RNA sequencing, the transcript is sequenced from the 3′ end. However, the 5′ cap modification cannot be identified from the data as a variant base call, and it is therefore not possible to determine if the transcript sequence is complete (extends to the original transcript origin) or not. Transcripts inferred from direct RNA sequencing are often truncated because of RNA degradation, electronic noise generated during the sequencing process and low signals near to the 5′ end (Soneson et al., 2019; Workman et al., 2019). We created P. falciparum transcript annotations by FLAIR analysis of the direct RNA sequencing data (Data S1). Full-length transcripts were identified from the correspondence of transcript 5′ end locations with 5CN from cDNA data (Table S5; Fig. 6). A total of 4,512 transcript 5′ ends (42%) corresponded with 5CN from one or more studies. The greatest number of transcripts annotated by FLAIR corresponded with 5CN from our data (2482), followed by 5CN from (Kensche et al., 2016; Adjalley et al., 2016; Chappell et al., 2020) with 1,720, 1,580, and 1,230 transcripts, respectively. Some transcript 5′ ends corresponded with 5CN from more than one study, but only 148 transcript 5′ ends corresponded with 5CN from all four studies (Fig. 6A). The majority of full-length transcripts overlap annotated genes/coding regions, with minorities of antisense and novel transcripts (Fig. 6B). Multiple full-length transcript isoforms with alternative 5′ ends are apparent for P. falciparum genes, including the PF3D7_1434200 (calmodulin) gene shown in Fig. 6C. The 5′ ends of two full-length coding transcripts map to a transcription initiation region upstream of the calmodulin gene mapped previously by RNA-ligase mediated 5′ rapid amplification of cDNA ends (Polson & Blackman, 2005). The untranslated regions of P. falciparum gene transcripts can extend for several hundred bases beyond the terminal coding exons (Chappell et al., 2020), such that the transcripts can overlap coding exons of adjacent genes. In the example shown in Fig. 6D, the 3′ untranslated region of a full-length coding transcript for the PF3D7_1214600 gene overlaps, and is antisense to the adjacent PF3D7_1214500 gene. Although the P. falciparum genome is gene-dense with median intergenic distance less than 2 kb (Russell et al., 2013), novel full-length transcripts not overlapping coding exons are present. In the example shown in Fig. 6E, a novel transcript initiates near to a full-length coding transcript of the PF3D7_0419600 gene on the opposite strand.