A comprehensive evaluation of rodent malaria parasite genomes and gene expression

Rodent malaria parasites (RMP) are used extensively as models of human malaria [1],[2]. Four different species that infect African rodents have been adapted for laboratory use: Plasmodium berghei, P. yoelii, P. chabaudi and P. vinckei. Small differences exist in the biology of the different RMP in laboratory mice and this makes them particularly attractive models to investigate different aspects of human malaria. Specifically, P. chabaudi is a model to investigate mechanisms of drug resistances and immune evasion, in particular antigenic variation [3],[4]. It invades normocytes and reticulocytes and mostly produces chronic, non-lethal, infections. In contrast, P. berghei preferentially invades reticulocytes and usually produces infections in mice that induce severe pathology [2]. In combination with different mouse strains it has been used as a model to study immunopathology, experimental cerebral malaria, pregnancy-associated malaria and lung pathology [2]. P. yoelii is widely used in studies on the biology of liver stages and on innate and acquired immunity against liver stages [5],[6]. Blood stage P. yoelii parasites of some lines are restricted to reticulocytes whereas others can invade all red blood cells and have been used to study receptors for erythrocyte binding [7],[8]. The availability of efficient reverse genetics technologies for P. berghei and P. yoelii [9]-[11] and the ability to analyse these parasites throughout the complete life cycle have made these two species the preferred models for analysis of Plasmodium gene function [12]-[14]. For these two species more than 600 different genetically modified mutants have been reported [15].

The first draft RMP genome was published in 2002 for P. yoelii yoelii 17XNL [16]. This was followed by publication of draft genomes of P. berghei ANKA (PbA) and P. chabaudi chabaudi AS (PcAS) in 2005 [17]. Comparisons with the genome of the human parasite P. falciparum and other primate malaria species defined a large set of core genes that are shared between RMPs and primate malarias [18]-[20]. Although availability of draft RMP genomes made a significant impact in applying post-genomic technologies for understanding malaria biology [18] and were used in many follow-up functional genomics studies to analyse gene regulation and function [9],[10], these RMP genomes were highly fragmented and were annotated with little or no manual curation. The fragmented nature of the genomes has hampered genome wide analysis of gene regulation and function, especially of the (subtelomeric) multigene families. To utilise RMP models to their full potential, we therefore undertook production of high quality reference genomes: for PbA and PcAS large-scale improvement of their existing genomes, with re-sequencing, re-analysis and manual re-annotation, and for P. y. yoelii a genome sequence was produced de novo from the virulent YM line using the latest sequencing technologies and computational algorithms. In addition, we have utilised comprehensive RNA-seq data derived from a number of life-cycle stages to both improve gene model prediction and to provide genome-wide, quantitative data on gene expression. By sequencing additional isolates/lines of P. berghei, P. yoelii and P. chabaudi (including the subspecies P. c. adami) we have documented genotypic diversity that exists within different RMP species. The availability of RMP reference genomes in combination with the RNA-seq and genotypic diversity data will serve as excellent resources for gene-function and post-genomic analyses and, therefore, better interrogation of Plasmodium biology and development of anti-malaria interventions.

The genomes of RMP contain a number of multigene families located in the subtelomeric chromosomal regions. These include a large family of so-called `Plasmodium interspersed repeat genes' (pir) [16], that are present also in other human/primate Plasmodium species [20]-[23]. Most of these gene families are expressed in blood stages and these proteins show features that have been reported to contribute to immune evasion through antigenic variation [24]-[26] and may play a role in the sequestration of infected red blood cells and virulence [26],[27]. As a result of the improved annotation, we have been able to define all multigene families in the RMP genomes. Comparative phylogenetic analyses of the pir genes and analyses of pir expression patterns in blood stages of P. berghei provide evidence of functional diversification within this gene family. The improved classification of multigene families will enhance studies on the role of (variant) exported proteins in virulence and evasion and modulation of the immune system.

By extensive re-sequencing and annotation we have generated three high quality RMP reference genomes with nearly all core genes as complete gene models and a much improved and almost complete representation of chromosomal subtelomeric regions. These reference genomes will greatly enhance the use of RMP as model organisms in malaria research. We provide full-length gene models for more than 98% of predicted protein-coding genes. The approximately 60% of genes with functional annotation is comparable to the percentage of functionally annotated genes in the P. falciparum 3D7 reference genome and a high percentage (approximately 90%) of the predicted RMP proteins have orthologs in primate malaria species. It is this high level of orthology between RMP and primate malaria genomes that strongly supports RMPs as models in experimental approaches to characterize the Plasmodium gene function. Similarly, the genome-wide RNA-seq data from different RMP developmental stages is a valuable resource to further analyse Plasmodium gene function and the regulatory networks underlying the multiple differentiation pathways of Plasmodium. The RNA-seq studies presented here provide information on gene expression at an unprecedented depth and breadth of coverage of multiple blood stages and ookinetes. Previously, only a few large-scale transcriptome (microarray) analyses of P. berghei blood stages and ookinetes had been performed [17],[46]. These studies were based on a highly fragmented draft P. berghei genome and, therefore, expression data were only generated for about half of all P. berghei genes. In addition, important/valuable large scale transcriptome studies have been performed on RMP life-cycle stages, such as sporozoites and liver stages [47]-[49]. These life-cycle stages would also benefit from re-examination using the latest RMP genome assemblies we provide in this study.

Our studies reveal that large scale changes in gene expression occur in ookinetes between 16 and 24 hours after fertilization, possibly required for the differentiation of (retort-form) zygotes into the mature ookinetes. The strong up-regulation in mature ookinetes of transcripts involved in ribosome biogenesis and protein translation suggest that the mature ookinete generates transcripts for proteins required after the ookinete has traversed the mosquito midgut wall and starts its rapid transition into the oocysts, possibly using mechanisms of translational repression similar to those in gametocytes [50],[51] and sporozoites [52],[53]. What these three stages have in common is that they are fully differentiated cells that will undergo rapid cellular differentiation and/or growth expansion upon entering a new environment. Whether mature ookinetes store repressed transcripts requires further investigation.

The additional sequence data from multiple RMP isolates will help to further unravel gene function and establish relationships between phenotypic traits and genotypic diversity. The near absence of polymorphisms within the genomes of P. berghei isolates was unexpected. Low sequence diversity of a limited number of genes of P. berghei isolates had been reported previously and it was proposed that this may result from cross-contamination of P. berghei isolates in the laboratory after isolation [54],[55]. However, this seems unlikely as one line would have needed to be mislabelled with the names of all other isolates, then all these mislabelled lines would have had to be sent to all the different laboratories worldwide replacing the `correct' isolates that may have existed in their recipient laboratories. However, sequencing of additional stocks from these isolates, which were frozen in different laboratories soon after isolation from the natural host, may reveal whether low sequence diversity is due to cross-contamination. The P. berghei isolates we have sequenced were obtained from other laboratories (SP11, NK65) and they also show a similar lack of sequence polymorphism. In contrast, the isolates of P. chabaudi exhibit considerable genotypic diversity. These P. chabaudi isolates exhibit differences in virulence, RBC invasion, growth rates and immunogenic profiles [7],[56]-[60] and further studies, for example using linkage or quantitative trait loci (QTL) analyses [33],[61]-[63], will facilitate identification of genes associated with defined phenotypes. For RMP species there is evidence that differences in virulence are associated with differences in RBC invasion [2],[7],[56],[60],[64]. For example, P. yoelii virulence has been associated with mutations in proteins involved in RBC invasion [7],[65],[66]. Interestingly, we found extensive sequence polymorphism in P. chabaudi genes encoding such proteins. While much attention is given to the role of exported proteins of multigene families and virulence in both human and RMP, further analysis of RMP proteins that regulate invasion phenotypes may reveal novel mechanisms that underlie virulence.

The new sequence data allowed for a much improved annotation of chromosomal subtelomeric regions and to better define the different subtelomeric multigene families. In addition to the large pir gene family, all three RMP contain an expanded gene family encoding exported proteins, fam-a, with orthology to a single-copy gene in primate malarias, which contains a START-domain (steroidogenic acute regulatory-related lipid transfer domain; [67]). START-containing proteins of eukaryotes are involved in the transfer of phospholipids, ceramide or fatty acids between membranes [68]. A START domain has also recently been identified in an exported, PEXEL-containing, P. falciparum protein that was shown to transfer phospholipids [69]. The single-copy RMP orthologs of this gene (PF3D7_0104200) also contain a PEXEL-motif, indicating that phospholipid-transporting proteins are exported into the RBC in both primate malarias and RMP. P. chabaudi contains an additional, highly expanded, gene family that contains domains involved in phospholipid/fatty acid metabolism. These genes, encoding putative lysophospholipases, lack a PEXEL motif; however, for several P. berghei orthologs as well as lysophospholipases of P. falciparum there is evidence for their export into the host RBC [29],[70]. Combined, these observations indicate the importance of phospholipid/fatty acid metabolism/transport mediated by Plasmodium proteins exported into the RBC cytosol. Why such genes have been differentially expanded into multigene families in different species remains to be investigated.

The pir family is the largest RMP multigene family and is shared with human and non-human primate species P. vivax, P. knowlesi and P. cynomolgi [20]-[23]. PIR proteins are exported into the RBC in the absence of a PEXEL-motif, and there is evidence that they are located on, or close to, the RBC surface or dispersed in the RBC cytoplasm [24]-[26],[29],[71],[72]. The function of pirs is unknown and no functional domains have been identified so far. Recently, it has been shown that in P. chabaudi a change in virulence was associated with differential expression of members of the pir multi-gene family [27]. It has been suggested that PIRs are transported to the surface of infected RBC and play a role in RBC sequestration comparable to the role of the Pfemp1 gene family of virulence factors in P. falciparum. However, for several P. berghei PIRs a direct role in RBC sequestration is unlikely since no evidence was found for their location on the RBC surface although they were exported into the RBC cytoplasm of both sequestering asexual blood stages and non-sequestering gametocytes [29]. For P. vivax PIRs it has been shown that different members have distinct subcellular locations in the infected RBC [26]. These observations indicate that functional differences may exist between members of the PIR family. Phylogenetic analyses support the possibility of functional differences between the PIRs. A recent phylogenetic analysis of the newly annotated PcAS pirs identified two distinct pir sub-families (A and B), which contain distinct amino acid sequence motifs [44]. Our phylogenetic analyses included pirs from all three RMP species and resulted in the identification of a number of different clades. The presence of clearly distinguishable clades indicates that structural differentiation exists among pirs and that this evolved prior to the separation of the RMP species. Our observations of the stage-specific up- or down-regulation of expression of clusters of structurally different pirs in different blood stages supports the hypothesis that there is functional diversification within the pir family and that purifying selection plays a role in shaping this family. By including multiple species in the pir phylogeny it is clear that this gene family is subject to rapid turnover, that is, gene gain and loss, indicating the absence of strong selective forces that would result in distinct orthologous groups/clades that are shared and maintained in different species for functional reasons. Gain of pir genes in different species is evident in the multiple species-specific expansions of clades. Assuming that the common ancestor had a pir family equal in abundance and diversity, the relatively limited instances of orthology (12 clades) indicates significant losses of ancestral sequence types. A plausible explanation for both the abundance of species-specific sequences and the paucity of ancestral sequences is a continual process of gene turnover driven by gene conversion, a mechanism that has been proposed for pirs of P. chabaudi [44] and which was evident in each of the clades revealed in this study (data not shown). The effect of frequent gene conversion is the replacement of ancestral sequence types with species-specific sequences, which results in distinct species-specific clades without orthology. Loss of orthology is only resisted when selective forces maintain structurally distinct pirs, which we propose, explains the presence of the (limited) orthology between pir clades of the different RMP species. The improved annotation and phylogeny demonstrating clusters of structurally different pirs in all RMP combined with expression profiles are powerful data that can help to further delineate function, the relationship of expression with virulence and how the (species-specific) expansion of the pirs is related to distinct selective pressures.

To maximise the utility of RMP we have greatly improved the genome assemblies of P. berghei and P. chabaudi, comprehensively sequenced the P. yoelii YM genome, sequenced multiple RMP isolates and generated in-depth expression data from multiple RMP life-cycle stages. Comparison of the RMP and P. falciparum genomes and RNA-seq mapping permitted gene annotation at base-pair resolution and has defined the level of orthology between RMP and human parasite genomes. The very high level orthology between RMP and human malarias (both in genome structure and gene content) supports the use of highly tractable RMPs as experimental models to characterize the function of the very many Plasmodium genes that remain uncharacterised.

Only a few large-scale transcriptome (microarray) analyses of different P. berghei life-cycle stages had previously been performed. Moreover, these studies were based on highly fragmented draft RMP genomes and consequently, for example, for one of the most well studied RMP, P. berghei, gene expression data was only mapped to about half of all P. berghei genes that have now been characterised. The RNA-seq studies we present provide information on gene expression, at an unprecedented depth and breadth of coverage, of multiple life cycle stages and provide the foundational data needed for the performance of large-scale analyses of gene regulatory networks that underlie cellular differentiation.

We show that extensive genotypic diversity exists between P. chabaudi isolates making this species an excellent organism to study genotype-phenotype relationships. Differences in virulence red blood cell (RBC) invasion, growth rates and immunogenic profiles exist between parasites of these isolates. Therefore, studies, such as linkage or quantitative trait loci analysis, are now possible to help identify genes associated with these defined phenotypes. For RMP species there is evidence that differences in virulence are associated with differences in RBC invasion, and, indeed, we find extensive sequence polymorphism in P. chabaudi genes encoding proteins involved in RBC invasion. Much attention is given to the role of exported proteins of multigene families and virulence in both human and RMP (for example, var, pirs), and analysis of differences between RMP proteins, that regulate invasion phenotypes, may reveal novel mechanisms that underlie virulence.

Full-length chromosomal annotation has permitted a comprehensive classification of all RMP subtelomeric multigene families. Our analyses indicate that most, if not all, RMP subtelomeric genes (apart from the RBP family) encode proteins exported out of the parasite; however, most lack a canonical PEXEL-motif. The presence of large numbers of PEXEL-negative exported proteins indicates alternative export mechanisms possibly common to all Plasmodium species. Investigations with highly tractable RMP species can therefore be used to understand these mechanisms better.

Our analyses of the phylogeny and expression of the largest RMP multi-gene family, the pirs, indicates functional diversification between members of the pir multigene family (this gene family is conserved between human/primate and RMP malaria species). Our new pir annotation and phylogeny demonstrates that clusters of structurally different pirs are differentially expressed. This is powerful data that can help to better understand their function, the relationship of pir expression with virulence and how the (species specific) pirs expansion is related to different selective pressures.

TDO, UB, MH, LR, SAO, CIN, AP, MB contributed to the genome assembly, annotation and analysis. TDO, BFF, WAMH, AAR, OB, SMK, HGS, APW, CJJ generated and analysed RNA-seq data of P. berghei. DC, JL, AAH provided materials for genome sequences and/or RNA-seq data of different isolates/lines of P. yoelii and P. chabaudi. AE provided materials for genomes of isolates of P. berghei. APJ performed the phylogenetic analysis. TDO, UB, APJ, BFF, WAMH, SMK, AAH, CIN, AP, MB, CJJ conceived the study, participated in its design and coordination and drafted the manuscript. All authors read and approved the final manuscript.