Plasmodium parasites, the causal agents of malaria, result in more than 1 million deaths annually. Plasmodium are unicellular eukaryotes with small ∼23 Mb genomes encoding ∼5200 protein-coding genes. The protein-coding genes comprise about half of these genomes. Although evolutionary processes have a significant impact on malaria control, the selective pressures within Plasmodium genomes are poorly understood, particularly in the non-protein-coding portion of the genome. We use evolutionary methods to describe selective processes in both the coding and non-coding regions of these genomes. Based on genome alignments of seven Plasmodium species, we show that protein-coding, intergenic and intronic regions are all subject to purifying selection and we identify 670 conserved non-genic elements. We then use genome-wide polymorphism data from P. falciparum to describe short-term selective processes in this species and identify some candidate genes for balancing (diversifying) selection. Our analyses suggest that there are many functional elements in the non-genic regions of these genomes and that adaptive evolution has occurred more frequently in the protein-coding regions of the genome.
Half of the world' s population is at risk of contracting malaria from Plasmodium species [1], so an understanding of their biology has considerable potential to influence human health. An understanding of evolution and natural selection are particularly important, because malaria control is limited by the evolution of resistance to anti-malarial drugs [2]–[5] and high levels of genetic variation in parasite surface proteins, which hinder natural and vaccine-induced immunity [6], [7]. An understanding of selection in these genomes can also contribute to our understanding of their function. For example, intronic and intergenic regions have been shown to be more conserved than neutrally evolving sites [8], [9], suggesting that purifying selection has been acting to conserve functional elements in non-genic regions of the genomes. In this study we describe selection in Plasmodium genomes over the long term using alignments of the genomes of seven Plasmodium species, showing that there is considerable constraint outside protein-coding regions and identifying 670 non-genic conserved elements. We describe selection in the short term using genome-wide polymorphism data from 13 strains of P. falciparum. This analysis suggested that protein-coding exons were more likely to be subject to non-neutral (adaptive) than non-genic regions.
The analysis of the divergence between distantly related species facilitates analysis of purifying selection (constraint). However, there are several types of selection acting simultaneously on genomes that are less well described from this data, particularly in non-genic regions that are difficult to align between highly divergent species. Firstly, adaptive (or directional) evolution, which acts to change the nucleotide sequence from its ancestral state [23]. There is also evidence that a few genes are subject to balancing selection, which acts to maintain multiple alleles (different versions of a gene) in a population [24]. Frequency-dependent balancing selection, where rare alleles have a selective advantage, is thought to be particularly important in Plasmodium genes encoding surface-exposed proteins that are targets of acquired immunity [25]–[29]. To examine these types of selection, we generated genome-wide genetic diversity (polymorphism) data from thirteen P. falciparum isolates (strains) and the divergence between P. falciparum and the chimpanzee parasite P. reichenowi, using publicly available ABI capillary reads (Table S5, NCBI dbSNP accessions: ss# 159747249–159815961). We called 69,805 SNPs within P. falciparum and 190,631 fixed differences between P. falciparum and P. reichenowi. We were able to obtain a coarse minor allele frequency (MAF) estimate for 54,641 SNPs and a derived allele frequency estimate for 24,573 SNPs (see Materials and Methods). See Table S6 for further summary statistics.
In summary, we show that there is considerable constraint in intergenic regions and introns. We identify 670 conserved non-genic elements and our analysis suggest that only a minority of these are un-annotated protein-coding exons, or structured ncRNAs. We suspect that many more functional non-genic elements remain undiscovered. Our analysis is consistent with the majority of non-neutral (directional or balancing) selection events having occurred in P. falciparum exons. Genetic diversity data collected from within populations and divergence data from more closely related Plasmodium species, both of which will soon be available, will be required to confirm this prediction.
The following genome versions were used for the alignment: P. falciparum version 2. 1. 4, July 2007, from ftp: //ftp. sanger. ac. uk/pub/pathogens/Plasmodium/falciparum/3D7/3D7. version2. 1. 4 Plasmodium knowlesi version PK4, October 2007, from ftp: //ftp. sanger. ac. uk/pub/pathogens/P_knowlesi/Archive/PK4. annotation/ Plasmodium vivax is as published in Carlton et al. [60], with ordered and orientated contigs in pseudo-chromosomes. Because many of the subtelomeric contigs could not be assigned to the pseudo-chromosomes, they are not present. However, these contigs are extremely AT-rich and contain mainly Vir genes, and so they do not align with chromosome regions of other Plasmodium species. Pseudo-chromosomes and annotation are available upon request. Plasmodium berghei, obtained April 2007, from ftp: //ftp. sanger. ac. uk/pub/pathogens/Plasmodium/berghei/ Plasmodium chabaudi, obtained April 2007, from ftp: //ftp. sanger. ac. uk/pub/pathogens/Plasmodium/chabaudi/ Plasmodium yoelii, PlasmoDB version 5. 4, from http: //www. plasmodb. org/common/downloads/release-5. 4/Pyoelii/ Plasmodium reichenowi A predicted Plasmodium reichenowi genome was created by aligning 78,442 P. reichenowi ABI capillary reads to the P. falciparum genome with SSAHA2, as described previously [9]. Fixed differences were discovered with a minimum phred score of 25. Deletions in the P. reichenowi sequence were identified, requiring at least two read alignments identifying an identical deletion before we accepted it. Insertions in P. reichenowi could not be included without manipulating the overall alignment. We then assumed that P. reichenowi was identical to P. falciparum except in locations of fixed differences or deletions in the P. reichenowi sequence. We excluded regions of the genome that were non-unique (as described below for SNPs) or lacked read coverage. The six assembled genomes were repeat masked using Repeat Masker Open-3. 0, from http: //www. repeatmasker. org, then aligned in two phases. The first phase was synteny mapping. Sub-sequences of the genomes were grouped into syntenic blocks. Each of the six genomes contributed at most one subsequence to a given block, and each block contained sequence from at least two species. Then, in the second phase, a nucleotide-level multiple alignment was constructed within each block. The synteny map was generated using Mercator [61]. Mercator requires anchor sequences along each genome. In addition, for each anchor, one must specify which anchors on other genomes are strongly similar. We chose as anchors the known and predicted exons in each genome, with the annotations obtained as described. Two anchors, each from a different genome, were deemed similar if the BLAT [62] score of the pair was below 1×10−50. The BLAT scores were computed in protein space. Mercator used the anchors and BLAT-similar pairs in a modified k-way reciprocal best hit algorithm [63]. The non-draft genome sequences (P. falciparum, P. knowlesi, and P. vivax) served as scaffolds for the draft species. The result was 170 syntenic blocks. The largest blocks covered parts of all six species and contained entire chromosomes from some of them, while the smallest contained small fragments of just two species. Mercator also produced alignment constraints for each block. A nucleotide-level multiple alignment within each block was generated with MAVID [64], using the alignment constraints as well as a phylogenetic tree relating the six species. Branch lengths for the tree were estimated with PAML [65], fixing the known topology for these Plasmodium species. Working upwards from the leaves, MAVID associates to each branch node a maximum-likelihood alignment of the sequences in the subtree rooted at that node. The alignment at the root of the full tree is a multiple alignment of all the input sequences. The accuracy and coverage of both the syntenic map and the block alignments were validated manually according to various descriptive statistics. The predicted P. reichenowi genome was then added to the alignment assuming complete synteny with P. falciparum. Alignment sites corresponding to exon, intron, and intergenic categories were extracted based on the annotation for P. falciparum, P. knowlesi and P. yoelii (one species from each of the three main clades), requiring identical annotation in all three species. For calculations of constraint we used two different models: The HKY85 model [13], with transition and transversion rates estimated globally for the whole tree, and a non-time reversible 12 parameter model (Nonrev), with all substitution parameters estimated individually for each of the three main branches in the tree (for a total of 36 substitution rates). To further ensure robustness we estimated the model parameters from three different data sets: A) The full alignment (representing overall selective pressures), B) introns, intergenic, and FFD sites (high divergence and AT content), and C) exons (lower divergence and AT content). Model parameters and branch lengths for the topology shown in Figure 1A were estimated with the maximum likelihood based package Hyphy [66]. The resulting six parameterized models were then used to estimate branch lengths based on the full alignment. Branch lengths for the different genomic regions (exon, intergenic, intron, FFD) were then estimated with each model as scalings of the full alignment tree (the relative divergence over all branches of the tree), by keeping the relative branch lengths within the tree fixed. Relative constraint for each region was calculated as, as described in the text. Variance of the constraint estimates was evaluated by 200 bootstrap replicates of alignment columns, with replacement. This variance is shown as error bars in Figure 2A, for each of the six models, and for each of the four categories (full alignment, exon, intergenic, intron). This measures the statistical variation of the constraint estimates, but does not address the biological variation. The latter is difficult to assess in a meaningful way, due to variation in neutral rate across the genome (Figure S4). The precision of estimates of neutral rate in specific regions of the alignment is limited by the relative scarcity of FFD sites, meaning that for smaller, biologically relevant window sizes, e. g. 20 kb, the neutral rate estimates will be highly variant (or non-existent) due to small FFD sample size, leading to uninformative constraint estimates. To estimate the constraint near gene and exon boundaries, we extracted alignment columns upstream of start codons (termed promoters), downstream of stop codons (terminators) and between exons (introns). In each case, we extracted a region corresponding to one third of the observed median length of the given feature. Introns were then further divided into donors (5′) and acceptors (3′). To reduce the effect of mis-annotations, we required all species to be present in the alignment, and identical, overlapping annotation features in all the six species for which annotation exists. When the exact boundaries for a feature on the alignment varied in different species' annotation, we chose the maximum start and minimum stop positions (shortest consensus). Results were very similar when we used the longest consensus (minimum start, maximum stop), a comparison can be seen in Figure S5. To examine whether divergence estimates from Plasmodium species alignments was negatively correlated with DAF, we located each SNP position in the alignment of Plasmodium species. We then calculated the median number of rejected substitutions estimated by GERP in both the 5nt and 11nt window around the SNP. Estimates of constraint within each main clade were calculated in three ways: First by the same procedure as described above for the full tree, but with each of the three main clades (Rodent, Primate, Ape) scaling independently. Second by estimating the models independently for each of these clades, excluding the long branches to the root. These clade-specific models were then used to estimate the sum of branch lengths (again excluding the long branches) for each of the categories (exon, intron, intergenic, FFD). Both these measures were calculated using regions that were consistently annotated (as exon, intron etc.) in P. falciparum, P. yoelii and P. knowlesi. Constraint was then calculated as described above. Third, we calculated constraint within clades as above using regions of the genome that were consistently annotated in all species within the clade (P. knowlesi and P. vivax for the primate clade, P. yoelii, P. berghei and P. chabaudi for the rodent clade, and only P. falciparum for the ape clade). CEs were identified with gerpcol and gerpelem programs of GERP [17], version 2. 1b (from http: //mendel. stanford. edu/sidowlab/downloads/GERP/index. html). Parameters for running gerpcol were estimated using Hyphy using a HKY85 model, fitted to the full alignment from FFD sites. These parameters were; a neutral rate of 3. 067, and transition/transversion rate of 2. 63, and the following phylogenetic tree ( ( (PlaBer: 0. 0839783, PlaYoe: 0. 127063): 0. 126016, PlaCha: 0. 147872): 0. 59425, (PlaViv: 0. 295231, PlaKno: 0. 213155): 0. 74449, (PlaFal: 0. 0141867, PlaRei: 0. 0420219): 0. 678922); . The gerpelem program of GERP was run using 10% false discovery rate and default parameters. The average proportion of substitutions that were separated by a multiple of 3 bases (PMOT) was calculated by comparing each sequence in the CE alignment slice to each other sequence. In each pairwise comparison we record the distance between each substitution, and then the proportion of these distances that were multiples of 3. The PMOT score for a CE is the average of these proportions from all pairwise comparisons. We determined the average tBLASTx bit score for each CE as follows. These BLAST searches are designed to detect protein-coding exons, and the median exon length for P. falciparum is 192 nt. So for each Plasmodium species represent in the CE alignment slice, we extract 192 nt of sequence extending out from the midpoint of the CE. We then used this 192 nt sequence to search the Alveolate genomes (Table S4) with NCBI tBLASTx (searched translated nucleotide databases using a translated nucleotide query), accepting only the best hit (irrespective of E-value). For each CE, we record the average bit score for each Plasmodium species represented in the alignment slice. Because all query sequences are the same length bit scores should be comparable between CEs. CEs were scored for RNA structure using RNAz version 1. 0 (http: //www. tbi. univie. ac. at/~wash/RNAz) [21]. Since some CEs were longer than the maximum length that RNAz can process, alignments were pre-processed with rnazWindow. pl, (default parameters except no reference sequence), before running RNAz. RNAz was run with default parameters, except predictions were done for both strands). Since rnazWindow. pl has a minimum length cut-off of 50 nt, predictions were not produced for the smaller elements. To give an estimate of the false discovery rate in our data, one simulated alignment was produced for each long CE with SISSIz (Version 0. 1) [22]. There were 8726 long exonic CEs where we could produce an RNAz prediction from both the native and the simulated alignment. From these, 98 of the SISSIz alignments had a support vector probability from RNAz >0. 95, and 527 of the native alignments had a SV probability >0. 95. So false discovery rate (FDR) is 98/527 = 0. 19. Similarly non-exonic FDR is 13/33 = 0. 4. Control elements for RNASeq expression analysis were chosen to match the length and GC content of each intergenic CE, but otherwise at random. This produced a set of 406 intergenic ‘non-conserved’ element controls whose length and GC content distributions did not differ significantly from a set of 489 P. falciparum intergenic CE elements (Mann-Whitney tests P>0. 1). We located a set of 120 intronic ‘non-conserved’ element controls in the same way to compare to the 80 intronic CEs in P. falciparum. There was no significant difference between the RNAseq expression levels of intergenic controls and intergenic CE elements, or between intronic controls and intronic CE elements (Mann-Whitney tests, both P>0. 05). Reads from 13 P. falciparum isolates and P. reichenowi were used to identify single nucleotide polymorphisms (SNPs) and fixed differences, as follows. See Table S5 for detailed information on SNP calls. All non-WTSI derived reads were downloaded from the NCBI Trace Archive ftp server (ftp. ncbi. nih. gov/pub/TraceDB/plasmodium_falciparum). Fastq files were mapped to the 3D7 version 2. 1. 4 reference sequence (ftp: //ftp. sanger. ac. uk/pub/pathogens/Plasmodium/falciparum/3D7/3D7. version2. 1. 4/) using SSAHA2, with the parameters (-seeds 15, -score 250, -tags 1, -diff 15, -output cigar, -memory 300, -cut 5000) and ssaha_pileup code. Ssaha2, ssaha_pileup and associated documentation are available from http: //www. sanger. ac. uk/Software/analysis/SSAHA2/. We excluded non-unique and tandem repeat regions of the genome. Non-unique regions were identified using the SSAHA2 read mapping score (mapping scores range from 0–50). We calculated the mean mapping score for each 2kb window of the genome (with a 100 nt step), using all reads aligned from all isolates (Smean). The 10th percentile of all Smean values was 49. 5 (Smean values were frequently 50). We excluded any region within 2kb of a window with Smean <49. 5. We identified tandem repeats using the Emboss application etandem (flags -minrepeat 2 -maxrepeat 10), and excluded all tandem repeats with ≥70% identity (to the other repeat). Excluding non-unique regions and tandem repeats left a remainder of 19,664,344 nt of unique non-repeat sequence, 84% of the genome. As an aid to determining suitable thresholds at which to accept SNP calls, we identified a set of 788 ‘reliable SNPs’ that were called in ≥6 isolates (of the 13 that we examined). Since it is unlikely that we call a false positive 6 times in an identical position, with an identical base, this set of SNPs is enriched for true calls. We then examined these reliable SNPs in each isolate (i) in turn. Each isolate will call a subset of these SNPs. Each SNP call may be from one or more reads, we determine the average phred score (for the SNP base) from all reads, PSNP (= total phred score/number of reads calling this base) for each reliable SNP in isolate i. We then examine the distribution of reliable SNP PSNP scores for isolate i, comparing it to the distribution of PSNP scores from all other SNPs called in isolate i. The distribution of PSNP scores from reliable SNPs were always significantly higher than the distribution of PSNP scores from all non-reliable SNPs (all other SNPs) due to more true calls in the reliable SNP set. Assuming that 95% of reliable calls are correct, we set the minimum phred scores required to call a SNP in isolate i as the 5th percentile of the distribution of reliable SNP PSNP scores. We refer to this value as min (PSNP, i). At each site in the genome where ssaha_pileup calls a SNP (in any isolate) we accept all SNP or reference base calls from each isolates that is supported by ≥2 reads with a PSNP ≥min (PSNP, i). When two alleles satisfied these criteria in an isolate, we accept both alleles, since some samples may contain >1 clone. In practice all min (PSNP, i) scores were ∼25, so SNPs are supported by at least 2 reads with cumulative phred ≥∼50. This method identifies reliable SNPs without a bias towards common SNPs. Fixed differences in P. reichenowi were accepted if supported by ≥2 reads with an average phred score for an allele ≥25. Derived allele frequencies were calculated from polymorphic sites with a P. reichenowi base call and ≥three isolate base calls, minor allele frequencies from polymorphic sites with ≥four isolate base calls (24,573 DAF calls, 54,641 MAF calls). For gene-specific analysis SNPs were assigned to a gene if they lay within the exons, introns or the half of adjacent intergenic region closer to the gene. We estimate that the false discovery rate for SNP calling is 1–2%, as follows. We aligned 30,840 reference 3D7 isolate reads from chromosome 12 onto the assembled 3D7 genome and called SNPs as above. Using the thresholds and filters we described above (including only unique regions and excluding tandem repeats), we accept 47 SNPs from 3D7 in chromosome 12 (Table S5). With a similar number of aligned reads, we accept 2,160 SNPs from the PFCLIN isolate and 2,583 SNPs from the IT isolate. If we assume that 47 of the 2160 SNPs from the PFCLIN isolate are false discoveries, then the false discovery rate (FDR) is 47/2160 = 2. 1%. With the same reasoning the IT isolate has 1. 8% FDR (Table S5). This is probably an overestimate of the FDR because a) some of the 3D7 calls may be correct (i. e. : errors in the reference sequence), and b) 3D7 reads are calling SNPs from a larger proportion of chromosome 12 (95% coverage at ≥2 read depth vs. ∼75% for IT and PFCLIN isolates, see Table S5). We also estimated the error rate of SNP allele calling using some Illumina data that was available for the IT and PFCLIN isolates (D. Jeffares, unpublished data). Briefly, 90 genes were chosen from primarily polymorphic but unique regions of the genome, and PCR-amplified from various isolates including PFCLIN and IT. Amplicons were sequenced to high depth with Illumina technology, and mapped to the same reference genome with MAQ. We examined how many of the PFCLIN and IT calls from the SSAHA2-mapped ABI capillary reads matched the alleles called (this study) matched those from the MAQ-mapped Illumina reads, using only sites covered by either 10 or 20 Illumina reads. Differences, which may be false SNP calls in either data set, were of the order of 1–2%, as predicted above. In general, error rates in different regions of the genome (exon, intron, intergenic, FFD, non-synonymous sites) were not significantly different. The exception was that for both isolates intergenic sites had significantly higher error rates than exonic regions (see Table S5). We expect this to result in an artifactual shift of intergenic sites to a lower allele frequencies, because artifactual alleles will be rare. We take this into account by comparing only intergenic vs. other intergenic sites. For the comparison of DAF/MAF this would be expected to diminish any affect of lower MAF distribution in exons vs. intergenic sites. Intronic, exonic, FFD and non-synonymous sites did not differ in error rates. For each isolate (i), a predicted genome was created, for each site in the genome we accept all SNP or reference base calls that were supported by ≥2 reads with a PSNP ≥min (PSNP, i). Sites without sufficient quantity coverage were denoted ‘N’ and not used in the analysis. Tajima' s D was calculated using Variscan (Version 2. 0, [67]), using a fixed number of alleles (4) for each SNP (Variscan chooses a random selection if >4 are available at a site), using only polymorphic sites (Variscan parameters FixNum = 1, NumNuc = 4, UseMuts = 0). The McDonald-Kreitman test neutrality index was calculated as NI = (Pn/Ps) / (Dn/Ds), where Pn and Ps are non-synonymous and synonymous SNPs and Dn and Ds are non-synonymous and synonymous fixed differences (between P. reichenowi and P. falciparum). It has been shown that the McDonald-Kreitman test can be used to estimate the average proportion of non-synonymous substitutions (α) that have been fixed by adaptive evolution [32], according to the formula, where Ds and Dn are the average number of synonymous and non-synonymous substitutions per gene and Pns is the average of per gene where Pn and Ps are the numbers of synonymous and non-synonymous polymorphisms respectively. This test can be generalised to use other classes of sites as the selected test, in place of non-synonymous sites in the original MK test [33]. We calculated a using four-fold degenerate sites as the neutral control and either non-synonymous sites (for exons), intronic sites, or intergenic sites, bootstrapping (by gene) 1000 times to determine the 5th and 95th percentiles. Intergenic SNPs and fixed differences were assigned to a gene if they fell in the half of the intergenic region closest to the gene. All statistics were performed in R (Version 2. 6. 0) (Ref. [68]). Tests for differences in DAF or MAF used Mann-Whitney U tests. Tests for differences in selective constraint between exon, intron, and intergenic sites within a gene used paired Mann-Whitney U tests.
Malaria causes debilitating ill-health in millions of people and kills about one million people annually, mostly young children. It is caused by a single-cell Plasmodium parasite transmitted to humans via mosquito bites. It is difficult to control this parasite because variable genetic make-up enables it to evade detection by vaccines and because drug resistance has repeatedly evolved. Therefore any progress in our understanding of the evolution and genetic variation of the parasite will be central to controlling the parasite. Genic regions that encode proteins are comparatively easy to characterize, whereas non-genic regions are poorly understood. We compare the genomes of seven distantly-related Plasmodium species and find that some of the non-genic regions are very similar between species. The absence of significant evolutionary differences between these non-genic regions implies that they play an important role in the survival of the organism. We then compare the genomes of thirteen different strains of Plasmodium falciparum. It is currently accepted that several families of antigenic parasite genes evolve rapidly. However, using two methods we demonstrate that many other genes have also undergone adaptive evolution.
genetics and genomics/population genetics
infectious diseases
genetics and genomics/comparative genomics
evolutionary biology/genomics
Long- and Short-Term Selective Forces on Malaria Parasite Genomes