Genomic and transcriptomic evidence for descent from Plasmodium and loss of blood schizogony in Hepatocystis parasites from naturally infected red colobus monkeys

Hepatocystis is a genus of single-celled parasites infecting, amongst other hosts, monkeys, bats and squirrels. Although thought to have descended from malaria parasites (Plasmodium spp.), Hepatocystis spp. are thought not to undergo replication in the blood–the part of the Plasmodium life cycle which causes the symptoms of malaria. Furthermore, Hepatocystis is transmitted by biting midges, not mosquitoes. Comparative genomics of Hepatocystis and Plasmodium species therefore presents an opportunity to better understand some of the most important aspects of malaria parasite biology. We were able to generate a draft genome for Hepatocystis sp. using DNA sequencing reads from the blood of a naturally infected red colobus monkey. We provide robust phylogenetic support for Hepatocystis sp. as a sister group to Plasmodium parasites infecting rodents. We show transcriptomic support for a lack of replication in the blood and genomic support for a complete loss of a family of genes involved in red blood cell invasion. Our analyses highlight the rapid evolution of genes involved in parasite vector stages, revealing genes that may be critical for interactions between malaria parasites and mosquitoes.


Introduction
Species of the genus Hepatocystis are single-celled eukaryotic parasites infecting, amongst other hosts, Old World monkeys, fruit bats and squirrels [1]. Phylogenetically, they are thought to reside within a clade containing Plasmodium species, including the parasites causing malaria in humans [2]. They were originally considered distinct from Plasmodium and have remained in a different genus because they lack the defining feature of asexual development in the blood, known as erythrocytic schizogony [3]. The presence of macroscopic exoerythrocytic schizonts (merocysts) in the liver of the vertebrate host is the most prominent feature of Hepatocystis [3]. Similar to Plasmodium parasites, Hepatocystis merocysts yield many single-celled merozoites. However, unlike Plasmodium, Hepatocystis merocysts appear to be the only replication phase in the vertebrate host [4]. First generation merozoites of Plasmodium spp. are released from liver cells and invade red blood cells, where they multiply asexually, before erupting from red cells as secondary merozoites. These merozoites invade further red blood cells before some develop into stages that can be transmitted to the vector. In contrast, liver merozoites of Hepatocystis spp. are thought to commit to the development of gametocyte transmission stages directly upon invading red blood cells. They are then vectored not by mosquitoes, but by biting midges of the genus Culicoides [5]. After fertilisation, Hepatocystis ookinetes encyst in the head and thorax of the midge between muscle fibres, whereas Plasmodium ookinetes encyst in the midgut wall of mosquitoes. After maturation, oocysts of both Plasmodium and Hepatocystis rupture and release sporozoites that migrate to the salivary glands [1]. These discrete biological differences, in the face of phylogenetic similarity and many shared biological features, make Hepatocystis a potentially powerful comparator for understanding important aspects of malaria parasite biology, such as transmission and host specificity.
A population of red colobus monkeys (Piliocolobus tephrosceles), from Kibale National Park, Uganda were previously shown to host Hepatocystis parasites based on morphological identification of infected red blood cells and DNA sequencing of the cytochrome b gene [6]. In this work we use Hepatocystis sp. genome and transcriptome sequences derived from P. tephrosceles whole blood samples to generate a draft genome sequence and gain insights into Hepatocystis evolution. We go on to use these insights to explore key aspects of malaria parasite biology such as red blood cell invasion, gametocytogenesis and parasite-vector interactions.

Genome assembly and annotation
While examining a published genomic sequence generated from the whole blood of a red colobus monkey (Piliocolobus tephrosceles) in Kibale National Park, Uganda (NCBI assembly ASM277652v1), we noticed the presence of sequences with significant similarity to Plasmodium spp. We hypothesised that this represented genomic material from a bloodborne parasite captured during the sequencing process and for each contig in the assembly we examined the AT-content and sequence similarity to Plasmodium spp. and macaque (Fig 1A). We identified a substantial subset of contigs that appeared to be derived from an apicomplexan parasite. Phylogenetic analysis, using an orthologue of cytochrome b from these contigs, suggested that they represent the first substantial genomic sequence from the genus Hepatocystis (Fig 1B), a parasite previously reported in Kibale National Park that infects at least four species of Old World monkeys [6]. At least four species of Hepatocystis are known to infect African monkeys-H. kochi, H. simiae, H. bouillezi and H. cercopitheci [6] -but with little sequence data currently linked to morphological identification, it was not possible to determine the species. We have thus classified the parasite as Hepatocystis sp. ex Piliocolobus tephrosceles (hereafter Hepatocystis sp.; NCBI Taxonomy ID: 2600580). The extraction of the Hepatocystis sp. sequences from the P. tephrosceles assembly yielded a set of 11,877 scaffolds with a total size of 26.26 Mb and an N50 of 2.4 kb. Automated genome annotation with Companion [7] identified 2,967 genes and 1,432 pseudogenes in these scaffolds. To improve upon this assembly, we isolated putative Hepatocystis sp. reads from the original short read DNA sequencing data. These were assembled into a draft quality nuclear genome assembly of 19.95 Mb, comprising 2,439 contigs with an N50 of 18.3 kb ( Table 1). The GC content (22.05%) was identical to that of Plasmodium spp. infecting rodents, and slightly higher than P. falciparum (19.34%). We identified 5,341 genes, compared to 5,441 in P. falciparum and 5,049 in P. berghei, suggesting a largely complete gene set (Table 1; S1 Table). Despite the fragmented nature of the assembly, we were able to identify synteny with P. falciparum around centromeres (S1 Fig)

Phylogenetic position of Hepatocystis sp. ex. Piliocolobus tephrosceles
There is consensus that Hepatocystis spp. are nested within the Plasmodium genus [2,8,9], however their placement within the genus has not been robustly determined. Indeed, our cytochrome b phylogeny confirms that our assembled genome is that of a Hepatocystis species, but it provides little support for the placement of this genus in relation to Plasmodium spp. A phylogeny generated using all mitochondrially encoded protein sequences also provided little support for key nodes (S3 Fig; see GitHub page for all sequence alignment data). The mitochondrial genome is therefore not reliable for determining the species phylogeny. A phylogeny based on 18 apicoplast proteins was more robust and placed Hepatocystis sp. as an outgroup to the Plasmodium species infecting rodents (Vinckeia; S4 Fig). We wanted to reliably place Hepatocystis sp. relative to other Hepatocystis species. Limited sequence data are available for Hepatocystis outside of this study. However, 11 nuclear genes have been sequenced for H. epomophori, a parasite of bats [2]. Based on the sequence of these genes, we found that Hepatocystis sp. forms a sister group to H. epomophori (S5 Fig). Furthermore, the Hepatocystis genus again forms a sister group to the Vinckeia subgenus of Plasmodium, although the tree contains some ambiguous branch points. To improve the robustness of the placement of Hepatocystis sp. within Plasmodium, we used 2673 orthologous nuclear genes from each of 12 species across the Plasmodium genus, which robustly places Hepatocystis sp. as a sister clade to the Plasmodium species infecting rodents (subgenus Vinckeia; Fig 2). Interestingly, some Vinckeia species (P. cyclopsi) also infect bats, supporting an earlier suggestion that the ancestor of (A) Contigs from the red colobus (Piliocolobus tephrosceles) assembly had a bimodal distribution of AT-content and sequence similarity to Plasmodium spp. (B). A phylogenetic tree of cytochrome b indicated that the closest match for the apicomplexan parasite sequenced from red colobus blood is a Hepatocystis isolate from a monkey host. Parasite cytochrome b sequences derived from RNA-seq assemblies from red colobus blood samples are almost entirely identical to the cytochrome b sequence assembled from Hepatocystis DNA reads from a single monkey. Branches of the tree have been coloured by bootstrap support values from 15 (red) to 100 (green). Some bootstrap support values are also shown next to the nodes as text. Red arrows highlight the Hepatocystis samples from the current study. Blue place names indicate the African continent, green Australia, orange Asia. https://doi.org/10.1371/journal.ppat.1008717.g001

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology Hepatocystis and Vinckeia might have infected bats [10]. However, whole-genome data also suggest that Vinckeia is derived from a group of monkey-infecting parasites [11]. Hepatocystis sp. has a long branch that could indicate rapid evolution after splitting from its ancestor with Vinckeia.

In vivo transcriptome data supports a lack of erythrocytic schizogony
Transcriptome sequencing of blood samples from 29 individuals was performed as part of the red colobus monkey genome sequencing project [12]. We found evidence that each of these individuals was infected with the same species of Hepatocystis as found in the genomic reads, consistent with high prevalence of this parasite in Kibale red colobus monkeys [6]. The extremely low SNP density suggested that parasites from different red colobus individuals were highly related ( Fig 3A). We identified an average of 1.36 SNPs (standard deviation = 0.47) and 0.43 indels (standard deviation 0.15) per 10 kb of genome when calling variants using RNA-seq reads.
Although it is believed that Hepatocystis spp. do not undergo erythrocytic schizogony [3], this has been challenged by limited microscopic evidence for asexual stages in the blood [13]. To determine whether there was transcriptomic evidence for schizonts in the blood we  Table). It is thought that chronic infections (of up to 15 months) may be maintained from continual development in the liver [3]. The presence of early blood stages in these individuals may therefore reflect this continual production of new blood forms, rather than recent infection. Proportions of rings and trophozoites were positively correlated and both these forms correlated negatively with female gametocytes (Fig 3C). Interestingly, the inferred proportions of male and female gametocytes were not strongly correlated suggesting there might be variation in commitment rates of gametocytes to male or female development.

Expanded and novel gene families
The largest gene family in Hepatocystis sp. was a novel family, which we have named Hepatocystis-specific family 1 (hep1; Table 2). These 12 single-exon genes (plus four pseudogenes) each encode proteins of~250 amino acids, beginning with a predicted signal peptide (S7A Fig). We could find no significant sequence similarity to genes from any other sequenced genome (using HHblits [14]). However, they contain a repeat region with striking similarity to that in Plasmodium kahrp, a gene involved in presenting proteins on the red blood cell surface [15]. Three members were highly expressed in in vivo blood stages, with one correlating well with the presence of early stages (S8A Fig; HEP_00211100, Pearson's r = 0.80 with rings).

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology Another novel family, hep2, contained N-terminal PEXEL motifs, suggesting it is exported (S7B Fig). Distinct parts of its sequence showed similarity (albeit with low significance) to exported proteins from P. malariae (PmUG01_00051800; probability: 77.88% HHblits) and P. ovale curtisi (PocGH01_00025800; 78.47%) as well as a gene in P. gallinaceum (PGAL8A_00461100; 69.52%). Three or four members were highly expressed in blood stages in vivo, with one member highly correlated with predicted proportions of early blood stages (S8B Fig; HEP_00165500, Pearson's r = 0.83 with rings; HEP_00480100, Pearson's r = 0.77 with rings). The gene encoding Thrombospondin-Related Anonymous Protein (trap) is involved in infection of salivary glands and liver cells by Plasmodium sporozoites. It is found strictly in a single copy in all Plasmodium species sequenced to date. However, it is present in six copies in Hepatocystis sp., suggesting that trap-mediated aspects of sporozoite-host interactions may be RNA-seq sample, relative to the genome assembly reference, highlight consistently low genetic diversity. Samples SAMN07757853, SAMN07757863, SAMN07757870 and SAMN07757873 have been excluded from the figure due to their low expression of Hepatocystis genes. (B) Deconvolution of RNA-seq samples to identify parasite stage composition shows no evidence for blood schizonts. Ring and trophozoite cells are assumed to relate to early stages of gametocyte development, which are not distinguishable from asexual rings and trophozoites. (C) Proportions of early blood stages (ring and trophozoite) are negatively correlated with mature female gametocytes, however male and female gametocyte ratios are poorly correlated, suggesting that sex ratios vary among samples. https://doi.org/10.1371/journal.ppat.1008717.g003

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology more complex. None of these genes were highly expressed in blood stage transcriptomes, consistent with their known role in sporozoites. Exported protein 1 (e.g. PF3D7_1121600) is a single copy gene in all Plasmodium species. It encodes a Parasitophorous Vacuolar Membrane (PVM) protein and is important for host-parasite interactions in the liver [16]. In Hepatocystis sp. it is expanded to four copies.

Missing orthologues tend to be involved in erythrocytic schizogony
The genomes of Plasmodium spp. each contain large families of genes known or thought to be involved in host parasite interactions [17]. These include, amongst others the var, rif and stevor genes in the Laverania subgenus, SICAvar genes in P. knowlesi and the pir genes across the genus. We find only four intact pir genes and a single pir pseudogene in Hepatocystis sp., compared to~200-1000 in Plasmodium spp. infecting rodents. One was particularly highly expressed in the blood stages from most monkeys that were sampled (S8C Fig;  HEP_00069900). All of these are most similar to the ancestral pir subfamily, present in single copy in Vinckeia species and in 19 copies in P. vivax P01 (  [18]. Given that asexual Hepatocystis parasites are not thought to exist in the blood of monkeys or bats, there would be no need for this function of pir genes. However, in Vinckeia, pir genes are expressed in several other stages, including male gametocytes [19], which do feature in the Hepatocystis lifecycle. The function of the ancestral pir gene subfamily is unknown, although it is expressed at multiple stages of the lifecycle in P. berghei [20]. We wanted to determine, more generally, the types of genes that might have been lost in Hepatocystis sp. relative to Plasmodium spp. This is made difficult due to uncertainty in determining missingness in a draft genome. We overcame this problem using clusters of genes identified as having common expression patterns across the lifecycle of the close Hepatocystis relative P. berghei [20]. Cluster 10 (late schizont expression; Fisher's Exact test odds ratio = 0.09, FDR = 0.0002) tended to contain orthologues shared by P. berghei, P. ovale wallikeri and P. vivax, but not Hepatocystis sp. (Fig 4; S3 Table; S10 Fig). Eleven out of 25 orthologues missing in the late schizont cluster (including two pseudogenes) encoded Reticulocyte Binding Proteins (RBPs). In fact, we could not identify any RBPs in the Hepatocystis sp. genome or Hepatocystis sp. RNA-seq assemblies. Interestingly, counter to previous suggestions [21], we found no evidence for RBP sequences in the Hae. tartakovskyi genome sequence. However, we did find fragmentary sequences with significant sequence similarity to RBPs from malaria parasites with avian hosts in transcriptome assemblies of Leucocytozoon buteonis [22] and Hae. columbae [23] (S11 Fig). Also missing from this cluster was Cdpk5, a kinase that regulates parasite egress from red cells [24]. The other principal gene families involved in erythrocyte binding and invasion by Plasmodium are the erythrocyte binding ligands (eh)/ duffy-binding protein (dbp) and the merozoite surface protein (msp) families [25]. These are largely conserved relative to P. berghei. Thus, orthologues missing relative to Plasmodium spp. tend to be involved in erythrocytic schizogony, the part of the life cycle also absent.

The most rapidly evolving genes are often involved in vector biology and control of gene expression
Our whole-genome phylogeny (Fig 2) showed a long Hepatocystis sp. branch, suggesting that some genes have changed extensively in Hepatocystis compared to Plasmodium spp. This might indicate functional changes important for the particular biology of Hepatocystis. We found previously that the ratio of synonymous mutations to synonymous sites (dS) saturates between Plasmodium clades [45] and therefore considered the ratio of non-synonymous

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology mutations (dN) rather than the more commonly used dN/dS. We first looked for enrichment of conserved protein domain families in genes with the highest 1% of dN values. There was an enrichment for the AP2 domain (Pfam:PF00847.20; Fisher test with BH correction; pvalue = 0.01). Most Plasmodium species possess 27 ApiAP2 transcription factors containing this domain, which are thought to be the key players in control of gene expression and parasite development across the life cycle. AP2-G plays an important role in exiting the cycle of schizogony and commitment to gametocytogenesis in Plasmodium spp. [26,27], whereas, as demonstrated, Hepatocystis sp. lacks erythrocytic schizogony. Hepatocystis spp. also form much larger cysts in the liver (giving the genus its name) and develop in different tissues within a different insect vector compared to Plasmodium spp. Our Hepatocystis sp. assembly contained orthologues of all 27 ApiAP2 genes present in P. falciparum (Table 2). This suggests that life cycle differences between Plasmodium and Hepatocystis spp. are not reflected in gain or loss of these key transcription factors. However, the relatively high rate of non-synonymous mutations suggests there may have been significant adjustment in how these transcription factors act. To determine parts of the life cycle that were enriched for the most rapidly evolving genes, we looked at whether particular gene expression clusters from the Malaria Cell Atlas [19,20] were enriched for genes with high dN (top 5% of values; Table 3; S12A Fig; S4 Table). We found that three clusters (2, 4 and 6) had fewer genes with high dN than expected by chance (Fisher's exact test with Holm multiple hypothesis testing correction, p-value < 0.05) and that these contained genes expressed across much of the life cycle, especially growth phases. These clusters also tended to contain essential genes expected to be highly conserved in Hepatocystis sp. Although there was not a significant trend for gametocyte-associated genes having higher than average dN, the top 5% of genes ranked by dN contained several putative gametocyte genes ( Table 3). Two of these encode putative 6-cysteine proteins P47 and P38, the first required for female gamete fertility [28]. Additionally, Merozoite TRAP-like protein (MTRAP), essential for Plasmodium gamete egress from erythrocytes [29] and two genes (HEP_00254800 and HEP_00195400) with orthologues involved in osmiophilic body formation [30,31] had high dN values. Overall, clusters involved in ookinete (15) and general mosquito stages (16) had significantly higher values than other clusters (Kolmogorov-Smirnov test: Cluster 15 vs all other clusters: D = 0.42, p-value = 1.05e-05. Cluster 16 vs all other clusters: D = 0.52, p-value = 4.50e-12; S12B Fig). This is also reflected in the Hepatocystis sp. genes with the highest dN values, which include oocyst rupture protein 2 (orp2; dN = 1.08), ap2-o, ap2-sp2, secreted ookinete protein (psop7) and osmiophilic body protein (g377). These genes provide clues about changes in the parasite that might relate to its adaptation to transmission by biting midges, rather than mosquitoes.

Discussion
We have taken advantage of parasite reads captured as part of a primate genome sequencing study in order to assemble and annotate a draft quality genome sequence for a species of Hepatocystis. Our nuclear and apicoplast genome phylogenies confirm the recently proposed phylogenetic placement of this genus as an outgroup to the rodent-infecting Vinckeia subgenus of Hepatocystis sp. orthologues of genes which are highly expressed in late schizogony in P. berghei tend to be missing from the genome. P. berghei gene expression clusters from the Malaria Cell Atlas were used to determine whether orthologous genes, conserved with P. vivax and P. ovale, but absent in Hepatocystis, tended to be expressed in particular parts of the life cycle. The only significant cluster was cluster 10, which includes genes most highly expressed in late schizont stages (highlighted by a red box). The top panels show the log 2 observed/ expected ratios for orthologous genes in each cluster which are shared between P. berghei and Hepatocystis and the same ratio for those which are shared between P. berghei, P. vivax and P. ovale, but not Hepatocystis sp. Cluster 1, 18, 19 and 20 were not tested because they contained fewer than 2 expected counts for either ratio. The asterisk indicates a Fisher's exact test false discovery rate of < = 0.05. https://doi.org/10.1371/journal.ppat.1008717.g004

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology Plasmodium [2]. However, a distinct branching pattern and low bootstrap support for many nodes in our mitochondrial genome phylogeny highlights why some previous analyses have come to different conclusions about the placement of Hepatocystis. Thus, the use of mitochondrial genes to infer phylogenetic relationships between species within the Haemosporidia should be approached with caution. We found a long branch leading to Hepatocystis sp., Table 3. Top 15 genes with functional annotations ranked by Hepatocystis sp. dN in comparison of Hepatocystis sp., P. berghei ANKA and P. ovale curtisi. Genes with completely unknown function and genes with very little information on their possible functions have been left out from this table. The rank column indicates the Hepatocystis sp. dN rank of each gene in the complete table (with 4009 genes) that includes genes with unknown function (S4 Table).

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology suggesting a relatively deep split from the rodent-infecting species. In addition, we showed robustly that Hepatocystis sp. clusters, as expected, with Hepatocystis epomophori, which infects bats. This finding supports the polyphyly of the Plasmodium parasites infecting apes, monkeys and rodents but the monophyly of Hepatocystis itself. A close relative of rodent-infecting P. berghei (P. cyclopsi) has been found in bats [10] and thus the Hepatocystis/Vinckeia group represents a relatively labile group with respect to host preference. Indeed, the possibility of crossspecies transmission of Hepatocystis sp. was reported previously [6].
The paraphyly of Plasmodium with respect to Hepatocystis exists because Hepatocystis spp. lack a defining characteristic of the Plasmodium genus, namely erythrocytic schizogony-asexual development in the blood. Thus, the very part of the Plasmodium life cycle which causes the symptoms of malaria is thought to be absent in Hepatocystis. While multiple lines of enquiry have failed to identify these forms (Garnham, 1966) there has remained some doubt, with reports of cells with schizont-like morphology in the blood for some species [13,44]. We were able to take advantage of bulk RNA-seq data collected from blood samples of a number of red colobus monkeys, all apparently infected with very closely related Hepatocystis parasites. By comparing to single-cell RNA-seq data from known cell types, we found no evidence for schizont stages in the blood. The apparent lack of schizonts could be due to sequestration of these stages away from the bloodstream. This is seen in some Plasmodium species, so we cannot rule out schizogony away from peripheral circulation. However, in its current state, this lack of evidence supports the idea that erythrocytic schizogony has been lost three times: in Hepatocystis, Polychromophilus and Nycteria [2]. Ancestrally, gametocytogenesis must have been the default developmental pathway, being required for transmission. However, in Plasmodium, it seems that erythrocytic schizogony became the default developmental pathway with epigenetic control of the ApiAP2-G transcription factor, required for development into sexual stages [26,27]. Perhaps the simplest explanation for a loss of erythrocytic schizogony would be that ApiAP2-G is no longer under control, but is constitutively expressed in parasites leaving the liver. In line with this idea we find ApiAP2-G present and highly expressed in blood stage Hepatocystis sp. Furthermore, all Plasmodium ApiAP2 transcription factors are conserved in Hepatocystis sp., indicating that changes in its life cycle are not associated with their loss or gain.
The lack of erythrocytic schizogony is supported by a tendency for orthologous genes missing in Hepatocystis sp. (relative to Plasmodium spp.) to be those expressed in blood schizonts. The most noticeable example is the complete absence of the Reticulocyte Binding Protein (RBP) family, found across all Plasmodium spp. examined so far, including those which infect birds [45]. RBP proteins are known to function as essential red blood cell invasion ligands in Plasmodium falciparum [46] and multiple copies are thought to provide alternative invasion pathways [25]. However, previous transcriptomic data have suggested that, while rbp genes in P. berghei are highly expressed in schizonts, they are less abundant or have a distinct repertoire in liver stages [20,47,48]. This implies that RBPs are less important, or at least that distinct invasion pathways are used by first generation merozoites. Also missing were cdpk5 (involved in schizont egress) and msp9 (an invasion related gene enriched in blood vs. liver schizonts in Caldelari et al. [47]). Taken together, these results underscore the increasing realisation that first generation merozoites have distinct properties from merozoites that have developed in the blood and suggest the way in which first generation merozoites invade red blood cells may be distinct. Looking more widely across the Haemosporidia, we found scant evidence for RBPs outside of the Plasmodium genus. There were no matches at all in the draft genome sequence of the relatively close Plasmodium outgroup Haemoproteus tartakovskyi. However, there were fragmentary matches to RBPs from Plasmodium species infecting birds in transcriptome assemblies of Hae. columbae and the more distant Leucocytozoon buteonis. The Hae. columbae fragment aligned to a conserved region of this divergent family, suggesting it could encode some core function of RBPs. More complete genome sequences from across the haemosporidians will be needed to understand the evolution of this gene family, which seems to be key for understanding the host specificity of Plasmodium parasites [49].
The genomes of Plasmodium spp. each contain large, rapidly evolving gene families that are known, or thought to be involved in host-parasite interactions, principally in asexual blood stages. The reason for their numbers may be due to a bet-hedging strategy, providing the diversity necessary for evading adaptive immune responses or dealing with unpredictable host variation. Although the Hepatocystis sp. genome contains two novel multigene families, we identified only 10-15 copies of each. The largest gene families in the closely related rodentinfecting Plasmodium species (pir and fam-a) are present here only as their ancestral orthologues, also conserved in the monkey-infecting Plasmodium species. We should be cautious in noting a lack of expansion in such families in Hepatocystis sp., as previous draft Plasmodium genome sequences have been shown to under-represent these genes. However, it is perhaps not surprising that pir genes are poorly represented. They are thought to play a role in the maintenance of chronic infections mediated by asexual stages in the blood [18] and Hepatocystis infection does not involve this stage of development. Given that Hepatocystis can sustain long chronic infections [1], presumably from the liver, this parasite may help us to better understand how Plasmodium survives in the liver.
A striking feature of the Hepatocystis life cycle is its vector-a biting midge rather than a mosquito. We found evidence of rapid evolution amongst orthologues of Plasmodium genes involved in mosquito stages of development, suggesting that adaptation to a new insect vector was a major evolutionary force. These rapidly evolving genes provide insights into parasitevector interactions and may provide avenues for the development of interventions to prevent transmission of the malaria parasite.
Overall, our findings demonstrate the insights that can be gained into malaria parasite biology from relatives of Plasmodium, even with draft-quality genome sequences. In the future, we expect that high-quality genome sequences of Hepatocystis spp., and additional relatives from genera such as Nycteria, Haemoproteus, and Polychromophilus, will be of great value for understanding the evolution and molecular biology of one of humanity's greatest enemies.

Ethics statement
All animal research was approved by the Uganda Wildlife Authority (permit number UWA/ TDO/33/02), the Uganda National Council for Science and Technology (permit number HS 364), and the University of Wisconsin-Madison Animal Care and Use Committee (V005039) prior to initiation of the study. Biological materials were shipped internationally under CITES permit #002290 (Uganda). The animal was anesthetized with a combination of Ketamine (5 mg/kg) and Xylazine (2 mg/kg) administered intramuscularly using a variable-pressure air rifle (Pneudart, Inc, Williamsport, PA, USA). After sampling, the animal was given a reversal agent (Atipamezole, 0.5 mg/kg), and released after recovery back to its social group. All animal use followed the guidelines of the Weatherall Report on the use of non-human primates in research.

Sample collection and data generation
The sequence data used in this study were part of a project originally designed to generate a reference genome for the red colobus monkey (genus Piliocolobus). Biomaterials used were from wild Ugandan (or Ashy) red colobus monkey (Piliocolobus tephrosceles) individuals from

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology Kibale National Park, Uganda. These animals reside in a habituated group that has been a focus of long-term studies in health, ecology, and disease [12,50,51]. Red colobus individuals were immobilised in the field as previously described [52]. Whole blood was collected using a modified PreAnalytiX PAXgene Blood RNA System protocol as described in Simons et al. [12]. Additionally, whole blood was collected into BD Vacutainer Plasma Preparation Tubes, blood plasma and cells were separated via centrifugation, and both were subsequently aliquoted into cryovials and stored in liquid nitrogen. Samples were transported to the United States in an IATA-approved liquid nitrogen dry shipper and then transferred to −80˚C for storage until further processing.
Methods for DNA extraction, library preparation, and whole genome sequencing are described in Simons [53]. Briefly, high molecular weight DNA was extracted from the blood cells of one red colobus monkey individual and size selected for fragments larger than 50,000 base pairs. A 10X Genomics Chromium System library preparation was performed and subsequently sequenced on two lanes of a 150 bp paired-end Illumina HiSeqX run as well as two lanes of a 150 bp paired-end Illumina HiSeq 4000 run.
Methods for RNA extraction and library preparation are described in Simons et al. (2019). Briefly, RNA was extracted from 29 red colobus individuals using a modified protocol for the PreAnalytiX PAXgene Blood RNA Kit protocol. Total RNA extracts were concentrated, depleted of alpha and beta globin mRNA, and assessed for integrity (RIN mean: 8.1, range: 6.6-9.2). Sequencing libraries were prepared using the KAPA Biosystems Stranded mRNA-seq Kit and sequenced on four partial lanes of a 150 bp paired-end Illumina HiSeq 4000 run. These data were uploaded to NCBI as part of BioProject PRJNA413051.

Separation of Hepatocystis and Piliocolobus scaffolds
The Piliocolobus tephrosceles genome assembly (ASM277652v1) was downloaded from the NCBI database. Scaffolds were first sorted by their GC% and Diamond 0.9.22 [54] BLASTX hits against a database of representative apicomplexan and Old World monkey proteomes. The sorting was improved by examining mapping scores of the scaffolds mapped to Plasmodium species and Macaca mulatta genomes (Mmul_8.0.1, GenBank assembly accession GCA_000772875.3) using Minimap2 2.12 [55]. The separation of scaffolds was further verified and refined by running NCBI BLAST of 960 bp fragments of all scaffolds against the NCBI nt database (Jul 18 2017 version) [56]. To predict genes in the apicomplexan scaffolds, Companion automatic annotation software [7] was run with these scaffolds as input and the P. vivax P01 genome as the reference.

Identification of Hepatocystis sequences in Piliocolobus RNA-seq data
Illumina HiSeq 4000 RNA-seq reads from the study PRJNA413051 were downloaded from the European Nucleotide Archive. In order to find out if the RNA-seq data contained apicomplexan sequences, mapping of these reads to apicomplexan scaffolds from Piliocolobus tephrosceles genome assembly (ASM277652v1) was done using HISAT2 2.1.0 [57].

Hepatocystis genome assembly
Filtering of reads for assembly. Minimap2 [55] and Kraken 2.0.8-beta [58] were used to identify the best matching species for each 10x Chromium genomic DNA read (from Illumina HiSeq X and HiSeq4000 platforms). Our Kraken database contained 17 Old World monkey genomes and 19 Plasmodium genomes downloaded from NCBI FTP in June 2018 [56]. The Kraken database also included the contigs of the P. tephrosceles assembly ASM277652v1, separated into P. tephrosceles and Hepatocystis sp. Plasmodium malariae UG01 (from PlasmoDB [59] version 39) and Macaca mulatta (Mmul_8.0.1) assemblies were used as reference genomes for the assignment of reads based on Minimap2 mapping scores. Reads that were unambiguously identified as monkey sequences using Kraken and Minimap2 were excluded from subsequent assemblies. The Supernova assembler manual [60] warns against exceeding 56x coverage in assemblies. Reads selected for Supernova assemblies were therefore divided into 34 batches, with~10 million reads in each batch. Reads were ordered by their barcodes so that those with the same barcode would preferentially occur in the same batch.
Supernova and SPAdes assemblies. We generated 34 assemblies with Supernova v2.1.1 with default settings. In addition, two SPAdes v3.11.0 [61] assemblies with default settings were generated with Hepatocystis reads: one with HiSeq X reads and another with HiSeq 4000 reads. Chromium barcodes were removed from the reads before the SPAdes assemblies.
Canu assembly. Scaffolds from the Supernova assemblies were broken into contigs. All contigs from the Supernova and SPAdes assemblies were pooled and used as the input for Canu assembler 1.6 [67] in place of long reads. Canu assembly was done without read correction and trimming stages. The settings for Canu were as follows: -assemble genomeSize = 23000k minReadLength = 300 minOverlapLength = 250 corMaxEvidenceErate = 0.15 correctedErrorRate = 0.16 stopOnReadQuality = false -nanopore-raw.
Processing of Canu unassembled sequences file. Selected contigs from the Canu unassembled sequences output file ( � .unassembled.fasta) were recovered and pooled with assembled contigs ( � .contigs.fasta). The first step in the filtering of the contigs of the unassembled sequences file was to exclude contigs that had a BLAST match in the assembled sequences output file (with E value cutoff 1e-10). Next, contigs where low complexity sequence content exceeded 50% (detected using Dustmasker 1.0.0 [68]) were removed. Contigs with GC content higher than 50% were also removed. Diamond BLASTX (against a database of Macaca mulatta, P. malariae UG01, P. ovale wallikeri, P. falciparum 3D7 and P. vivax P0 proteomes) and BLAST (using the nt database from Jul 18 2017 and nr database from Jul 19 2017) were then used to exclude all contigs where the top hits were not an apicomplexan species. In total, 0.34% of contigs from the unassembled sequences file were selected to be included in the assembly.
Deduplication of contigs. Initial deduplication of contigs was done using BBTools dedupe [69] (Nov 20, 2017 version) and GAP5 v1.2.14-r3753M [70] autojoin. In addition, BUSCO 3.0.1 [71] was used to detect duplicated core genes with the protists dataset. Two contigs flagged by BUSCO as containing duplicated genes were removed. All vs all BLAST of contigs (with E-value cutoff 1e-20, minimum overlap length 100 bp, minimum identity 85%) was used to find possible cases of remaining duplicated contigs. Contigs yielding BLAST hits were aligned with MAFFT v7.205 [72] and the alignments were manually inspected. Contained contigs were deleted and contigs that had unique overlaps with high identity were merged into consensus sequences using Jalview.
Removal of contaminants after Canu assembly. All Canu assembly contigs were checked with Diamond against a database of Macaca mulatta, P. malariae UG01, P. ovale wallikeri, P. falciparum 3D7 and P. vivax P01 proteomes. The Diamond search did not detect any contaminants. Contigs not identified by Diamond were checked with BLAST against the nt database (Jul 18 2017 version). Contigs where the top BLAST hit was a human or monkey sequence were removed from the assembly.
A subset of contigs in the assembly was observed to consist of short sequences with low complexity, high GC% and low frequency of stop codons. These contigs did not match any sequences by BLAST search against nt and nr databases (with E-value cutoff 1e-10). Due to their difference from the rest of the contigs in the assembly, it was assumed that these contigs were contaminants rather than Hepatocystis sequences. In order to programmatically find these contigs, GC%, tandem repeats percentage, percentage of low complexity content and frequency of stop codons were recorded for all contigs in the assembly. Tandem Repeats Finder 4.04 [73] was used to assess tandem repeats percentage and Dustmasker 1.0.0 [68] was used to find low complexity sequence content. PCA and k-means clustering (using R version 3.5.1) showed that the assembly contigs separated into two groups based on these parameters. The group of contigs with low complexity (189 contigs) was removed from the assembly.

Curation and annotation of the Hepatocystis genome assembly
The assembly was annotated using Companion [7]. The alignment of reference proteins to target sequence was enabled in the Companion run but all other parameters were left as default. A GTF file derived from mapping of Hepatocystis RNA-seq reads of three biological samples (SAMN07757854, SAMN07757861 and SAMN07757872) to the assembly was used as transcript evidence for Companion. To produce the GTF file, the RNA-seq reads were mapped to the assembly using 2-pass mapping with STAR RNA-seq aligner [80] (as described in the "Variant calling of RNA-seq samples" section) and the mapped reads were processed with Cufflinks [81]. All Plasmodium genomes available in the web version of Companion were tested as the reference genome for annotating the Hepatocystis genome, in order to find out which reference genome yields the highest gene density. For the final Companion run the P. falciparum 3D7 reference genome (version from June 2015) was used. The Companion output was manually curated using Artemis [82] and ACT [83] version 18.0.2. Manual curation was carried out to correct the overprediction of coding sequences, add missing genes and correct exon-intron boundaries. Altogether 680 gene models were corrected, 546 genes added and 221 genes

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology deleted. RNA-seq data was used as supporting evidence. Non-coding RNAs were predicted with Rfam [84].
All genes were analysed for the presence of a PEXEL-motif using the updated HMM algorithm ExportPred v2.0 [85]. Distant homology to hep1 and hep2 gene families was sought by using the HHblits webserver with default options [14].

Analysis of other Haemosporidian genomes and transcriptomes
Genomes and transcriptomes of other Haemosporidians. A Haemoproteus tartakovskyi genome assembly was downloaded from the Malavi database (http://130.235.244.92/Malavi/ Downloads/Haemoproteus_tartakovskyi). The Companion annotation tool [7] was used to automatically annotate the assembly, using the P. falciparum 3D7 genome as the reference. The alignment of reference proteins to target sequence was enabled and the rest of the settings were left as default. The genome annotation was further edited manually using Artemis 18.1.0 [82].
BLAST searches for RBPs. A BLAST database was made from Plasmodium RBPs downloaded from PlasmoDB [96] (release 46). NCBI blastx 2.9.0+ with e-value cutoff 1e-5 was run against this database, using the Haemoproteus genome and transcriptome assemblies and the Leucocytozoon transcriptome assembly as queries. The transcripts that yielded BLAST hits were further examined using BLAST against the NCBI nt and nr databases (April 2020 versions) and by sequence alignments with RBPs from PlasmoDB. One of the Haemoproteus transcripts

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology (GGWD01016989.1) matched RBPs in the blastx search (the best match was with PRELSG_0 014300: E-value 1.31e-23, score: 89.0) and did not yield non-RBP BLAST hits in searches against the nt and nr databases. One Leucocytozoon transcript also matched Plasmodium RBPs (top match: PRELSG_0013000, E-value 1.06e-09, score: 56.2). The sequence of the Haemoproteus transcript GGWD01016989.1 was translated to amino acids using ExPASy Translate (April 2020 version) [97] and then aligned with a selection of Plasmodium RBPs from PlasmoDB (release 46) using MAFFT 7 [98] with default settings. The resulting alignment was cropped in Jalview 2.10.4b1 to include only the region that contained the Haemoproteus sequence [99]. A phylogenetic tree was generated from the alignment as described in the next section.

Phylogenetic trees
Haemosporidian sequences were downloaded from NCBI FTP and PlasmoDB (release 43). The phylogenetic tree of cytochrome B and the tree that included 11 Hepatocystis epomophori genes were based on DNA alignments. The cytochrome B tree also included cytochrome B sequences from de novo assemblies of Hepatocystis RNA-seq reads derived from Piliocolobus tephrosceles blood. The trees of mitochondrial, apicoplast and nuclear proteomes were based on protein alignments. For apicoplast proteome and nuclear proteome trees, orthologous proteins were identified using OrthoMCL 1.4 [100]. The OrthoMCL run included the Haemoproteus tartakovskyi proteome that had been derived from the Haemoproteus tartakovskyi genome assembly using Companion, as previously described. All vs all BLAST for OrthoMCL was done using blastall 2.2.25 with E-value cutoff 1e-5. OrthoMCL was run with mode 3. Proteins with single copy orthologs across all the selected species were used for the protein phylogenetic trees. Sequences were aligned with MAFFT 7.205 [98] (with-auto flag) and the alignments were processed using Gblocks 0.91b [101] with default settings. Individual Gblocks-processed alignments were concatenated into one alignment. The phylogenetic trees were generated using IQ-TREE multicore version 1.6.5 [102] with default settings and plotted using FigTree 1.4.4 (https://github.com/rambaut/figtree/releases). Inkscape (https://inkscape. org) version 0.92 was used to edit text labels of the phylogenetic trees generated with FigTree.

Clustering of pir proteins into subfamilies
Sequences of Plasmodium pir family proteins (including bir, cyir, kir, vir and yir proteins) were downloaded from PlasmoDB [59] (release 39). The sequences were clustered using MCL [103], following the procedures described in the section "Clustering similarity graphs encoded in BLAST results" in clmprotocols (https://micans.org/mcl/man/clmprotocols.html). The BLAST E-value cutoff used for clustering was 0.01 and the MCL inflation value was 2. The pir protein counts per subfamily in each species were plotted as a heatmap using the heatmap.2 function in gplots package version 3.0.1.1 in R version 3.5.1.

Mapping and assembly of Hepatocystis sp. RNA-seq data
To separate Hepatocystis reads from Piliocolobus reads, RNA-seq data from the ENA (study PRJNA413051) were mapped to a FASTA file containing genome assemblies of Hepatocystis and M. mulatta (NCBI assembly Mmul_8.0.1), using HISAT2 version 2.1.0 [57], with "-rnastrandness RF". BED files were generated from the mapped reads using BEDTools 2.17.0 [63]. Reads from each technical replicate were merged, resulting in a single set of read counts for each individual monkey. The BED files were filtered to remove multimapping reads and reads with mapping quality score lower than 10. Names of reads that specifically mapped to the Hepatocystis assembly were extracted from the BED file. SeqTK 1.0-r31 (https://github.com/ lh3/seqtk) was used to isolate Hepatocystis FASTQ reads based on the list of reads from the

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology previous step. The Hepatocystis reads were then mapped to the Hepatocystis genome assembly using HISAT2 2.1.0 with "-rna-strandness RF" flag. The SAM files with mapped reads were converted to sorted BAM files with SamTools 0.1.19-44428cd [62]. The EMBL file of Hepatocystis genome annotations was converted to GFF format using Artemis 18.0.1 [82]. Htseqcount 0.7.1 [104] was used to count mapped reads per gene in the GFF file with "-t mRNA -a 0 -s reverse". Htseq-count files of individual RNA-seq runs were merged into a single file.
In order to extract Hepatocystis cytochrome b sequences of each RNA-seq sample, Hepatocystis RNA-seq reads of each sample were isolated from Piliocolobus reads as described above and then assembled with the SPAdes assembler v3.11.0 [105] with the "-rna" flag. Hepatocystis cytochrome b contigs were identified in each of the 29 RNA-seq assemblies using BLAST against Hepatocystis cytochrome b from the DNA assembly (E-value cutoff 1e-10).
In addition to assemblies of individual RNA-seq samples, an assembly of all RNA-seq samples pooled was done. The reads for this assembly were sorted by competitive mapping to P. ovale curtisi GH01 (from PlasmoDB release 45) and Macaca mulatta (Mmul_8.0.1, GenBank assembly accession GCA_000772875.3) genomes with Minimap2 (with the "-ax sr" flag). Reads mapping to the Macaca mulatta genome with minimum mapping score 20 were removed and the rest of the reads were assembled with the SPAdes assembler v3.13.1 [105] with the "-rna" flag. Hepatocystis contigs were identified by comparison of sequences with Plasmodium and Macaca mulatta reference genomes using Diamond, Minimap2 and BLAST, similarly to what is described in the section "Separation of Hepatocystis and Piliocolobus scaffolds". Further decontamination was done using Diamond and BLAST searches against 19747 sequences from Ascomycota and 165860 bacterial sequences downloaded from UniProt (release 2019_10) [106] and 3 Babesia proteomes from PiroplasmaDB (release 46) [107]. Selected contigs were also checked with BLAST against the NCBI nt database. The assembly was deduplicated using BBTools dedupe (Nov 20, 2017 version) and GAP5 v1.2.14-r3753M. Assembly completeness was assessed using CEGMA 2.5. In order to reduce the number of contigs so that they could be used as input for Companion, the assembly was scaffolded with RaGOO Version 1.1 [108], using the Hepatocystis DNA assembly as the reference. The assembly was then processed by the Companion annotation software (Glasgow server, November 2019 version, with P. falciparum 3D7 reference genome, with protein evidence enabled and the rest of the settings left as default). In order to detect proteins missed by Companion, EMBOSS Transeq (version 6.3.1) was used to translate the transcriptome assembly in all 6 reading frames. The output of Transeq was then filtered to keep sequences between stop codons with minimum length of 240 amino acids. Protein BLAST with E-value cutoff 1e-20 was used to detect sequences in Transeq output that were not present in the proteins annotated by Companion. These selected Transeq output sequences were checked for contaminants with BLAST similarly to what was described before. The sequences that passed the contaminant check were combined with the set of Hepatocystis RNA-seq assembly proteins that were detected by Companion. OrthoMCL was run with proteins from Hepatocystis RNA-seq assembly (Companion and selected Transeq sequences combined), Hepatocystis DNA assembly proteins, 20 Plasmodium proteomes from PlasmoDB release 43 and P. ovale wallikeri proteome (GenBank GCA_900090025.2). The settings for OrthoMCL were as described in the "Phylogenetic trees" section. dN analysis P. berghei ANKA and P. ovale curtisi protein and transcript sequences were retrieved from PlasmoDB [59] (release 45). One-to-one orthologs between Hepatocystis, P. berghei ANKA and P. ovale curtisi were identified using OrthoMCL [100]and a Newick tree of the three species was generated with IQ-TREE [102]. The settings for OrthoMCL and IQ-TREE were as

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology described in the "Phylogenetic trees" section. Transcripts of one-to-one orthologs were aligned using command line version of TranslatorX [109] with "-p F -t T" flags, so that each alignment file contained sequences from three species. Gaps were removed from alignments while retaining the correct reading frame. Alignment regions where the nucleotide sequence surrounded by gaps was shorter than 42 bp were also removed. In addition, the script truncated alignments at the last whole codon if a sequence ended with a partial codon due to a contig break. The alignments and the Newick tree of the 3 species were then used as input for codeml [110] in order to determine the dN and dN/dS of each alignment. The codeml settings that differed from default settings were: seqtype = 1, model = 1. P. berghei RNA-seq cluster numbers from Malaria Cell Atlas [20] were assigned to each alignment based on the P. berghei gene in the alignment. Transcriptomics-based gametocyte specificity scores of Plasmodium genes were taken from an existing study on this topic [111] (transcripts S2 Table of "Transcriptomics_all_studies" tab). The P. falciparum genes in the gametocyte specificity scores table were matched with equivalent Hepatocystis genes using OrthoMCL (run with the same settings as when used for phylogenetic trees). Statistical tests with the dN results (Kolmogorov Smirnov test, Fisher test and Spearman correlation) were performed using the stats library in R.

Variant calling of RNA-seq samples
SNPs and indels were called in Hepatocystis RNA-seq reads that had been separated from Piliocolobus reads as described above. Four technical replicates of each RNA-seq sample were pooled. Variant calling followed the "Calling variants in RNAseq" workflow in GATK [112] user guide (https://gatkforums.broadinstitute.org/gatk/discussion/3891/calling-variants-inrnaseq). First, the reads were mapped to the reference genome using 2-pass mapping with the STAR RNA-seq aligner [80] version 2.5.3a. 2-pass mapping consisted of indexing the genome with genomeGenerate command, aligning the reads with the genome, generating a new index based on splice junction information contained in the output of the first pass and then producing a final alignment using the new index. GATK [112] version 4.0.3.0 was used for the next steps. The mapped reads were processed with GATK MarkDuplicates and SplitNCigarReads commands. GATK HaplotypeCaller was then run with the following settings:-dont-use-softclipped-bases-emit-ref-confidence GVCF-sample-ploidy 1-standard-min-confidencethreshold-for-calling 20.0. Joint genotyping of the samples was then done using GATK Combi-neGVCFs and GenotypeGVCFs commands. This was followed by running VariantFiltration with these settings: -window 35 -cluster 3-filter-name FS -filter 'FS > 30.0'-filter-name QD -filter 'QD < 2.0'. SNPs were separated from indels using GATK SelectVariants. Samples SAMN07757853, SAMN07757863, SAMN07757870 and SAMN07757873 were excluded from further analysis due to their low expression of Hepatocystis genes (htseq-count reported below 50,000 reads mapped to the Hepatocystis assembly in each of these samples). The average filtered SNP counts per 10 kb of reference genome for each sample were calculated as the number of filtered SNPs divided by (genome size in kb � 10).

RNA-seq deconvolution
Deconvolution of a bulk RNA-seq transcriptome sequence aims to determine the relative proportions of different cell types in the original sample. This requires a reference dataset of transcriptomes from "pure" cell types. To create this, we used single-cell P. berghei transcriptome sequences from the Malaria Cell Atlas [20]. For each cell type, single-cell transcriptome sequences were combined by summing read counts per gene to generate a set of pseudobulk transcriptome sequences (see our GitHub repository). The aim of summing across cells is to reduce the number of dropouts which are common in individual single-cell transcriptome sequences. Bulk Hepatocystis RNA-seq transcriptome sequences, mapped and counted as above, were summed across replicates and filtered to exclude those with fewer than 100,000 reads. Hepatocystis and P. berghei pseudobulk read counts were converted to Counts Per Million (CPM) and Hepatocystis gene ids were converted to those of P. berghei one-to-one orthologues. Genes without one-to-one orthologues (defined by orthoMCL analysis) were excluded. CIBERSORT v1.06 [113] was used to deconvolute the Hepatocystis transcriptomes with the MCA pseudobulk as the signature matrix file. To test the accuracy of this deconvolution process we generated mixtures of the pseudobulk resulting in e.g. equal representation of read counts from male gametocyte, female gametocyte, ring, trophozoite and schizont pseudobulk transcriptomes (see our GitHub repository). We also deconvoluted bulk RNA-seq transcriptomes from Otto et al. [89] processed as in Reid et al. [19].

Enrichment of missing genes in Malaria Cell Atlas gene clusters
We wanted to determine whether there were functional patterns common to orthologues missing from the Hepatocystis genome relative to Plasmodium species. To do this we looked for orthologous groups (orthoMCL as above) containing genes from P. berghei, P. ovale wallikeri and P. vivax P01, but not Hepatocystis. Genes from P. berghei have previously been assigned to 20 clusters based on their gene expression patterns across the whole life cycle [20]. We looked to see whether missing orthologues tended to fall into particular clusters more often than expected by chance (see our GitHub repository). We used Fisher's exact test with Benjamini-Hochberg correction to control the false discovery rate. We reported clusters with FDR > = 0.05. Exons are shown in coloured boxes with introns as linking lines. '//' represents a gap. The shaded/grey areas in P. knowlesi and P. falciparum mark the start of the conserved, syntenic regions to other Plasmodium species. The presence of genes that are subtelomeric in Plasmodium species, i.e. PHIST proteins, suggests that the Hepatocystis scaffolds are also subtelomeric. A complete subtelomere that includes telomeric repeats is missing in our Hepatocystis assembly. Thus, whether Hepatocystis chromosomes retain the organisation common to most Plasmodium species remains unclear. Genes of Hepatocystis sp. ex Piliocolobus tephrosceles are highly similar to Hepatocystis epomophori genes sequenced in a different study [2]. The tree is based on the following genes: splicing factor 3B subunit 1, tubulin gamma chain, DNA polymerase delta catalytic subunit, eukaryotic translation initiation factor 2 gamma subunit, T-complex protein 1 subunit alpha, pantothenate transporter, ribonucleoside-diphosphate reductase large subunit, aminophospholipid-transporting P-ATPase, GCN20, transport protein Sec24A and RuvB-like helicase 3. (MCA) gene cluster 10 represents genes highly expressed in late schizonts. 25 genes from this cluster were conserved in P. ovale wallikeri and P. vivax, but were missing from our Hepatocystis genome assembly. Genes were clustered here by expression pattern and single-cells were ordered by pseudotime as in [20]. (B) MCA cluster 4 represents genes highly expressed across much of the life cycle-liver stages, trophozoites, female gametocytes and ookinetes/oocysts. 27 genes from this cluster were conserved in P. ovale wallikeri and P. vivax, but were missing from our Hepatocystis genome assembly.  Table. Summary of gene properties. For each gene in the assembly, the following is listed: annotation, number of exons, gene length (bp), the presence or absence of start and stop codons (reflecting the completeness of the assembly of the gene) and RNA-seq expression level (mean FPKM with standard deviation) in sample SAMN07757854 (RC106R). For the proteins encoded by the genes, the table shows the number of transmembrane segments predicted by TMHMM, ExportPred 2 score, 1 to 1 orthologs in P. berghei ANKA and P. ovale curtisi GH01 (based on OrthoMCL), PFAM domains, the number of matches to PEXEL motif (RxLxE/Q/ D) and SignalP-5 signal peptide prediction. (XLSX) S2

PLOS PATHOGENS
Hepatocystis genome illuminates malaria parasite biology