Rapid identification of genes controlling virulence and immunity in malaria parasites

Identifying the genetic determinants of phenotypes that impact disease severity is of fundamental importance for the design of new interventions against malaria. Here we present a rapid genome-wide approach capable of identifying multiple genetic drivers of medically relevant phenotypes within malaria parasites via a single experiment at single gene or allele resolution. In a proof of principle study, we found that a previously undescribed single nucleotide polymorphism in the binding domain of the erythrocyte binding like protein (EBL) conferred a dramatic change in red blood cell invasion in mutant rodent malaria parasites Plasmodium yoelii. In the same experiment, we implicated merozoite surface protein 1 (MSP1) and other polymorphic proteins, as the major targets of strain-specific immunity. Using allelic replacement, we provide functional validation of the substitution in the EBL gene controlling the growth rate in the blood stages of the parasites.


Introduction
Malaria parasite strains are genotypically polymorphic, leading to a diversity of phenotypic characteristics that impact on disease severity. Discovering the genetic basis for such phenotypic traits can inform the design of new drugs and vaccines. Both association mapping and linkage analyses approaches have been adopted to understand the genetic mechanisms behind various phenotypes of malaria parasites [1][2][3][4][5] and with the application of whole genome sequencing (WGS), the resolution of these methodologies has been dramatically improved, allowing the discovery of selective sweeps as they arise in the field [6]. However, both approaches suffer from drawbacks when working with malaria parasites: linkage mapping requires the cloning of individual recombinant offspring, a process that is both laborious and time-consuming, and association studies require the collection of a large number of individual parasites (usually in the thousands) from diverse geographical origins and over periods of several months or years to produce enough resolution for the detection of selective sweeps.
Linkage Group Selection (LGS), like linkage mapping, relies on the generation of genetic crosses, but bypasses the need for extracting and phenotyping individual recombinant clones. Instead, it relies on quantitative molecular markers to measure allele frequencies in the recombinant progeny and identify loci under selection [7,8]. This approach bears similarity to Bulked Segregant Analysis (BSA) [9], a technique developed to study disease resistance in plants. In BSA, individuals from a population are segregated based upon their phenotype (e.g. disease resistance), following which the frequencies of genetic markers in each population are analysed, identifying loci at which different alleles are found for the differently phenotypes populations. Segregating individuals by phenotype, while relatively straight forward for large organisms such as plants, is not feasible for unicellular pathogens such as malaria parasites. Instead, in LGS, the segregating population is grown both in the presence or absence of a selection pressure (e.g. drug treatment, immune pressure, etc.). Selection removes susceptible individuals in the selected "pool", while leaving both susceptible and resistant individuals in the unselected "pool". In its original implementation, LGS was successfully applied in studying strain-specific immunity (SSI) [10,11], drug resistance [7,12] and growth rate [8] in malaria and SSI in Eimeria tenella [13].
LGS is essentially identical to the extreme QTL approach (xQTL) that was independently developed by yeast researchers based on BSA [14].
In both the original implementations of BSA and LGS a limiting factor is the availability of molecular markers differentiating the two populations. One step in increasing the number of molecular markers was through the use of array hybridisation that allowed the identification of thousands of SNPs as molecular markers in Arabidopsis thaliana [15]. BSA (still using preselected pools) was also combined with tiling microarray hybridisation and used probe intensities to detect a gene underlying xylose utilisation in yeast [16]. The xQTL method increased the power and rapidity of the approach by making use of available yeast microarray data as well as Next Generation Sequencing (NGS) of DNA hybridised to microarray probes to identify a large number of markers across the genome, this time comparing selected and unselected populations, rather then generating pools based on phenotype [14,17]. In the absence of microarray databases, an alternative approach was to use NGS short reads to identify genomewide SNPs between two parents and then use these SNPs as molecular markers to identify target genes in the selected progeny population compared against the unselected population, as done to study chloroquine resistance in malaria [18].
In this study, we apply an improved LGS approach for the identification of genes controlling two independent and naturally occurring phenotypic differences between two strains of the rodent malaria parasite Plasmodium yoelii; growth rate, and strain-specific immunity. A mathematical model, built upon methodological improvements in the analysis of genetically crossed populations [19,20], was developed to analyze the data. This modified LGS approach relies on the generation and selection of at least two independent crosses between the strains. The progeny from both crosses pre-and post-selection are then subjected to high-throughput WGS, and SNP marker movement analyzed using best fitting modeling (Fig 1).
Our novel statistical framework both accounts for the influence of clonal growth in the cross population, and allows for a locally variable recombination rate in the parasite population, unlike previous analyses applied to comparable data [21]. Applying this framework to crosses between two strains of P. yoelii that induce SSI, and which differ in their growth rates, we were able to identify three genomic regions and alleles controlling both phenotypes, demonstrating that the approach can be used to analyze multiple complex phenotypes concomitantly with high genomic resolution within a short space of time.

Results
Characterization of strain specific immunity and growth rate phenotypic differences between CU and 17X1.1pp The difference in blood-stage parasite growth rate between the two clones was followed in vivo for nine days in CBA mice. A likelihood ratio test using general linear mixed models indicated a more pronounced growth rate for 17X1.1pp compared to CU clone by time interaction term, L = 88.60, df = 21, p<0.0001, Fig 2A). To verify that the two malaria clones could also be used to generate protective SSI, groups of mice were immunized with 17X1.1pp, CU or mock immunized, prior to challenge with a mixture of the two clones (S1 Fig). The relative proportions of the two clones were measured on day four of the infection by real time quantitative PCR (Q-RT-PCR) targeting the polymorphic msp1 locus [22]. A strong, statistically significant SSI was induced by both parasite strains in CBA mice (Fig 2B).

Identification of high-confidence SNPs
Two kinds of selection pressure were applied in this study: growth rate driven selection and SSI. Two independent genetic crosses between 17X1.1pp and CU were produced, and both these crosses were subjected to immune selection (in which the progeny were grown in mice The process starts with the identification of distinct selectable phenotypes in cloned strains of the pathogen population (in this case malaria parasites) and their sequencing, usually from the vertebrate blood stage. A genetic cross between two cloned strains is subsequently produced, in this case inside the mosquito vector. The cross progeny is then grown with and without selection pressure(s). Selection pressure will remove those progeny individuals carrying allele(s) associated with sensitivity to the selection pressure(s), while allowing progeny individuals with the resistant allele(s) to made immune to either of the two parental clones), and grown in non-immune mice. Progeny were harvested from mice four days after challenge, at which point strain-specific immune selection in the immunized mice, and selection of faster growing parasites in the non-immune mice had occurred. Using deep sequencing by Illumina technology, a total of 29,053 high confidence genome-wide SNPs that distinguish the parental strains were produced by read mapping with custom-made Python scripts. SNP frequencies from these loci from each population were filtered using a likelihood ratio test to remove sites where alleles had been erroneously mapped to the wrong genome location.
Identification of clonality within the data A hidden Markov model was applied to the data to identify allele frequency changes ( Table 1) that were likely to have arisen from the clonal growth of individuals within the cross population or possible incorrect assembly of the reference genome, as described in the Materials and Methods section and in more detail in the supplementary mathematical methods (S1 Appendix). In a genetic cross population, an especially high fitness clone generated by random recombination events can grow to substantial frequency, this being manifested as sudden jumps in allele frequency occurring at the recombination points in this individual [23]. Jumps of this type were primarily identified in the 17X-immunized population, where the increased virulence of the 17X strain had less of an effect in driving alleles to high frequency, and in the first replica experiment; the data in the first experiment seemed to have been more affected by clonal growth in the population. The consistency of identified jumps between treatment conditions reflects the common origin of the differently treated populations; the jump at the end of chromosome XIV inferred in both replicas may be artefactual.
survive. DNA is then extracted from the whole, uncloned progeny for sequencing. SNPs distinguishing both parents are used to measure allele frequencies in the selected and unselected progenies. A mathematical model is then applied to identify and define loci under selection. Regions in these loci are then analyzed in detail to identify potential target polymorphisms underlying the phenotype(s) under investigation. Targeted capillary sequencing can be employed to verify or further characterize polymorphisms. Finally, where applicable, allele replacement experiments can be carried out to confirm the effect of target polymorphisms. https://doi.org/10.1371/journal.ppat.1006447.g001

Identification of loci under selection
Based upon an analytical evolutionary model describing patterns of allele frequencies following selection, a maximum likelihood approach was used to define confidence intervals for the positions of alleles under selection in each of the genetic cross populations. In the absence of selection acting for a variant in a region of the genome, the allele frequencies in that region are Table 1. Sudden changes in allele frequency identified using a jump-diffusion model. Details are given for loci at which a sudden jump in frequency was inferred with probability at least 1%. The latter value is the inferred probability that the change in allele frequency at a given locus arose from a jump to a random position between 0 and 1, as opposed to arising from a small change to the frequency at the previous locus. Data are shown for the naïve and 17-X immunized experiments; no jumps of this significance were inferred for the CU-immunized experiment.

Condition
Replica expected to be locally constant. In common with a previous approach to identifying selected alleles [21], a search was therefore made for regions of the genome in which allele frequencies varied substantially according to their position in the genome. Next, wherever deviations of this form were consistently identified in both replica experiments a model of selection was applied to the data, inferring for each set of replica data the position in that region of the genome that was most likely to be under selection; this model was based upon expected changes in allele frequency under a constant local rate of recombination and is described further in the Methods section. Regions of the genome in which this inference of selection produced consistent results across replica datasets were then identified ( Table 2). Of a total of 11 genomic regions suggesting evidence of non-neutrality, six showed sufficient evidence of consistent selection. For each of these regions of the genome, a more sophisticated evolutionary model, accounting for variation in the local recombination rate, was then applied to the data, refining the position of the putatively selected allele. At this point, a putative selected allele in chromosome IV was removed from consideration, leaving five cases of potential alleles under selection in three regions of the genomes; confidence intervals for the positions of the selected loci are given in Table 3. Optimal positions of variant loci derived from each replicate are detailed in S1 Table; results of the variable recombination rate model are shown in S2 Table, with inferred recombination rates in S3 Table. Of the final three putative loci, two were detected under multiple experimental conditions (Fig 3). When considering the combined largest intervals, a selective sweep was inferred at position 1,436-1,529 kb on Chromosome (Chr) XIII in replicate crosses grown in both nonimmunized mice and 17X1.1pp-immunized mice, resulting from selection against CU-specific Table 2. Identification of candidate regions by non-neutrality score and SD model selected allele location. The non-neutrality score for region in replica r is denoted S r . The optimal driver location in the same region is given by i Ã r . Where a chromosome is divided into parts, by potential jump alleles, the resulting genomic regions are denoted by their chromosome number, a subscript indicating which part of the genome was under consideration. Identified candidate regions were defined as those at which selection was identified at positions within 200kb in both replicates, and are here highlighted in bold type.

Condition
Chromosome Score S  alleles at the target locus. A second sweep was inferred at position 1,229-1,364 kb on Chr VIII, detected in the parasite crosses grown in both CU and 17X1.1pp immunized mice, though not in the non-immunized mice. Here, selection pressure acted against different alleles according to the strain against which mice were immunized. The third sweep was detected at a locus between positions 725-814 kb on Chr VII. This event was only detected in mice replicates immunized with the 17X1.1pp strain, albeit that a consistent change in allele frequencies was also observed between replicas grown under these conditions (Fig 3B). The remaining loci (on Chrs VIII and XIII) were not consistently detected between replicates (S1 Table) and were thus considered to be non-significant.

Potential target genes within the three main loci under selection
All the genes in the combined conservative intervals of the three main loci under selection are listed in S4-S6 Tables, along with annotation pertaining to function, structure, orthology with P. falciparum genes and Non-synonymous/Synonomous SNP (NS/S) ratio in the P. falciparum orthologue, which is calculated by the PlasmoDB website (6.2) based on SNP data from 202 individual strains. These include both laboratory strains and field isolates obtained from six collections (see Methods for more details). The locus associated with SSI on Chr VIII contains 41 genes. We considered the presence of either transmembrane (TM) domains or a signal peptides as necessary features of potential antigen-encoding genes. Only 16 genes met these criteria. Functional annotation indicated 10 likely candidates among these; eight genes described as "conserved Plasmodium proteins", and two encoding RhopH2 and merozoite surface protein 1 (MSP1). Of these genes, the P. falciparum orthologue of msp1 had the highest NS/S SNP ratio (8.43). MSP1 is a well characterized major antigen of malaria parasites that has formed the basis of several vaccine studies [24] and has been previously linked to SSI in Plasmodium chabaudi [10][11][12].
The locus under selection on Chr VII consists of 21 genes. Only seven contained TM domains and/or a signal peptide motif. Based on functional annotation, four of these could be potential targets for SSI. One of these genes, PY17X_0721800, encodes an apical membrane protein orthologous to Pf34 in P. falciparum. This protein has recently been described as a surface antigen that can elicit an immune response [25]. Three conserved proteins of unknown function (PY17X_0720100, PY17X_0721500 and PY17X_0721600) also displayed potential signatures as target antigens.
The growth rate associated selected locus on Chr XIII contains 29 genes. In this case, the presence of TM domains or signal peptide motifs were not considered informative criteria. Only eight genes contained NS SNPs between the parental strains 17X1.1pp and CU according to the WGS data. Among these was a duffy binding protein, Pyebl. Pyebl, is a gene that has been previously implicated in growth rate differences between strains of P. yoelii [8,26]. A single NS SNP was predicted from the WGS data in this gene. Due to the very high likelihood of its involvement based on previous work, this gene was considered for further analysis.
Characterization of EBL as the major driver of growth rate differences through allelic replacement Examining the Pyebl gene, Sanger capillary sequencing re-confirmed the existence in 17X1.1pp of an amino acid substitution (Cys >Tyr) at position 351 within region 2 of the encoded protein. When aligned against other P. yoelii strains and other Plasmodium species, this cysteine residue is highly conserved, and the substitution observed in 17X1.1pp was novel ( Fig 4A). Crucially, no other polymorphisms were detected in the coding sequence of the gene, including in region 6, the location of the SNP previously implicated in parasite virulence in other strains of P. yoelii [8]. Structural modeling of the EBL protein in both wild-type and 17x1.1pp (C351Y) mutants predicted the abolition of a a disulphide bond between C351 and C420 in the mutant parasites that alters the tertiary structure of the receptor binding region of the ligand in these parasites (Fig 4B and 4C).
The functional role of this polymorphism was verified by experimental means. In order to study the functional consequences of the polymorphism, the Pyebl alleles of slow growing CU and faster growing 17X1.1pp clones were replaced with the alternative allele (i.e. CU-EBL-351C>Y and 7x1.1pp-EBL-351Y>C), as well as with the homologous allele (i.e. CU-EBL-351C>C and 17x1.1pp-EBL-351Y>Y). The latter served as a control for the actual allelic swap, as the insertion of the plasmid for allelic substitution could potentially affect parasite fitness independently of the allele being inserted. To establish whether the C351Y substitution affected EBL localization, as was shown for the previously described region 6 mutation, Immunoflurescence Analysis (IFA) was performed. This revealed that, unlike the known mutation in region 6 [8], the EBL proteins of 17X1.1pp and CU were both found to be located in the micronemes (Fig 5 and S2 Fig).
Transgenic clones were grown in mice for 10 days alongside wild-type clones. Pair-wise comparisons between transgenic clones with the parental allele against transgenic clones with the alternative allele (that is CU-EBL-351C>C vs CU-EBL-351C>Y and 17x1.1pp-EBL-351Y>Y vs 17x1.1pp-EBL-351Y>C) showed that allele substitution could switch growth phenotypes in both strains (Fig 6A and 6B). This confirmed the role of the C351Y mutation as underlying the observed growth rate difference.
RNA-seq analysis revealed that transfected EBL gene alleles were expressed normally, (S3 Fig), thus indicating a structural effect of the polymorphism on parasite fitness, rather than an alteration in protein expression.

Discussion
The development of LGS has facilitated functional genomic analysis of malaria parasites over the past decade. In particular, it has simplified and accelerated the detection of loci underlying selectable phenotypes such as drug resistance, SSI and growth rate [7,8,10]. Here we present a radically modified LGS approach that utilizes deep, quantitative WGS of parasite progenies and the respective parental populations, multiple crossing and mathematical modeling to identify loci under selection at ultra-high resolution. This enables the accurate definition of loci under selection and the identification of multiple genes driving selectable phenotypes within a very short space of time. This modified approach allows the simultaneous detection of genes or alleles underlying multiple phenotypes, including those with a multigenic basis.
Applying this modified LGS approach to study SSI and growth rate in P. yoelii, we identified three loci under selection that contained three strong candidate genes controlling both phenotypes. Two loci were implicated in SSI; the first time LGS has identified multigenic drivers of phenotypic differences in malaria parasites in a single experimental set-up. The strong locus under selection in Chr VIII, associated with the gene encoding MSP1, is consistent with existing knowledge of malaria immunity. The Chr VII locus, which includes the orthologue of Pf34 as well as other potential unannotated antigens, underscores the power for hypothesis generation and gene detection of the LGS approach using multiple crosses.
Our approach also provided a genetic rationale for the difference in growth rate of the parental clones CU and 17X1.1pp. Phenotypically, this occurs due to the ability of 17X1.1pp to invade both reticulocytes and normocytes, while CU is restricted to reticulocytes [22]. Previously, differences in growth rates between strains of P. yoelii have been linked to a polymorphism in Region 6 of the Pyebl gene that alters its trafficking so that the protein locates in the dense granules rather than the micronemes [8,26]. In the case of 17x1.1pp however, direct sequencing of the Pyebl gene revealed a previously unknown SNP in region 2, the predicted receptor-binding region of the protein, with no polymorphism in region 6. Consistent with this, the EBL protein of 17X1.1pp was shown to be located in the micronemes, indicating that protein trafficking was unaffected by the region 2 substitution. Allelic replacement of the parasite strains with the alternative allele resulted in a switching of the growth rate to that of the other clone, thus confirming the role of the substitution.
Region 2 of the Pyebl orthologues of P. falciparum and Plasmodium vivax [27][28][29] are known to interact with receptors on the red blood cell (RBC) surface. Furthermore, the substitution falls within the central portion of the region, which has been previously described as being the principal site of receptor recognition in P. vivax [29]. Wild-type strains of P. yoelii (such as CU) preferentially invade reticulocytes but not mature RBCs, whereas highly virulent strains are known to invade a broader repertoire of RBCs [30]. Further structural and functional studies are required to elucidate how the polymorphism described here enables mutant parasites to invade a larger repertoire of erythrocytes than wild type parasites. We show that the cysteine residue at position 351 in EBL forms a disulphide bond with a cysteine at position 420, and that this is abolished following the C351Y substitution, altering the tertiary structure of the binding region. This leads to the possibility that such an alteration of the shape of the binding domain may enable the ligand to bind to a larger repertoire of receptors.

Rapid identification of functional genes in malaria parasites
LGS with multiple crosses offers a powerful and rapid methodology for identifying genes or non-coding regions controlling important phenotypes in malaria parasites and, potentially, in other apicomplexan parasites. Through bypassing the need to clone and type hundreds of individual progeny, and by harnessing the power of genetics, genomics and mathematical modeling, genes can be linked to phenotypes with high precision in a matter of a few months, rather than years. Here we have demonstrated the ability of LGS to identify multiple genetic polymorphisms underlying two independent phenotypic differences between a pair of malaria parasite strains; growth rate and SSI. This methodology has the potential power to identify the genetic components controlling a broad range of selectable phenotypes, and can be applied to studies of drug resistance, transmissibility, virulence, host preference, etc., in a range of apicomplexan parasites that are amenable to genetic crossing.
The applicability of the approach to human malaria species has been recently demonstrated: the original LGS approach was successfully applied to study P. falciparum immune evasion in mosquitoes in vivo [31], while we recently tested its applicability in vitro to detect loci under selection following antifolate drug treatment and in vitro growth rate competition. With the advent of humanized mice that are able to support the complete malaria life cycle, the generation of new genetic crosses between strains of human malaria has become more feasible, as recently demonstrated [32]. With the ability to maintain these crosses without the need of simian hosts, application of a broader range of selection pressures (excluding, for now, selection mediated by the presence of a complete immune response) is now more feasible in vivo, thus extending the application of the LGS approach to medically relevant malaria species.

Materials and methods
Full and complete details of the mathematical methods are given in S1 Appendix, Supplementary Mathematical Methods.

Parasites, mice and mosquitoes
Plasmodium yoelii CU (with slow growth rate phenotype) and 17X1.1pp (with intermediate growth rate phenotype) strains [33] were maintained in CBA mice (SLC Inc., Shizuoka, Japan) housed at 23˚C and fed on maintenance diet with 0.05% para-aminobenzoic acid (PABA)-supplemented water to assist with parasite growth. Anopheles stephensi mosquitoes were housed in a temperature and humidity controlled insectary at 24˚C and 70% humidity, adult flies were maintained on 10% glucose solution supplemented with 0.05% PABA.

Testing parasite strains for growth rate and SSI
Plasmodium yoelii parasite strains were typed for growth rate in groups of mice following the intravenous inoculation of 1 × 10 6 iRBCs of either CU, 17X1.1pp or transfected clones per mouse and measuring parasitaemia over 8-9 days. In order to verify the existence of SSI between the CU and 17X1.1pp strains, groups of five mice were inoculated intravenously with 1 × 10 6 iRBCs of either CU or 17X1.1pp parasite strains. After four days, mice were treated with mefloquine (20mg/kg/per day, orally) for four days to remove infections. Three weeks post immunization, mice were then challenged intravenously with 1 × 10 6 iRBCs of a mixed infection of 17X1.1pp and CU parasites. A group of five naïve control mice was simultaneously infected with the same material. After four days of growth 10 μl of blood were sampled from each mouse and DNA extracted.
Strain proportions were then measured by Quantitative Real Time PCR using primers designed to amplify the msp1 gene [34]. All measurements were plotted and standard errors calculated using the Graphpad Prism software (v6.01) (http://www.graphpad.com/scientificsoftware/prism/). Wilcoxon rank sum tests with continuity corrections were used to measure the SSI effect, and were performed in R [35]. Linear mixed model analyses and likelihood ratio tests to test parasite strain differences in growth rate were performed on log-transformed parasitaemia by choosing parasitaemia and strain as fixed factors and mouse nested in strain as a random factor, as described previously [22]. Pair-wise comparisons of samples for the transfection experiments were performed using multiple 2-way ANOVA tests and corrected with a Tukey's post-test in Graphpad Prism software (v6.01).
Preparation of genetic cross. Plasmodium yoelii. CU and 17X1.1pp parasite clones were initially grown separately in donor mice. These parasite clones were then harvested from the donors, accurately mixed to produce an inoculum in a proportion of 1:1 and inoculated intravenously at 1 × 10 6 infected red blood cells (iRBCs) per mouse into a group of CBA mice. Three days after inoculation, the presence of gametocytes of both sexes was confirmed microscopically and mice were anesthetized and placed on a mosquito cage containing 400 female A. stephensi mosquitoes six to eight days post emergence. Mosquitoes were then allowed to feed on the mice without interruption. Seven days after the blood meal, 10 female mosquitoes from this cage were dissected to examine for the presence of oocysts in mosquito midguts. Seventeen days after the initial blood meal, the mosquitoes were dissected, and the salivary glands (containing sporozoites) were removed. The glands were placed in 0.2-0.4 mL volumes of 1:1 foetal bovine serum/Ringer's solution (2.7 mM potassium chloride, 1.8 mM calcium chloride, 154 mM sodium chloride) and gently disrupted to release sporozoites. The suspensions were injected intravenously into groups of CBA mice in 0.1 mL aliquots to obtain blood stage P. yoelii CU17X1.1pp cross progeny. Three days after inoculation with sporozoites, blood stage P. yoelii CU17X1.1pp cross progeny parasitized-RBC (pRBC) were harvested.
Two independent genetic crosses between CU and 17X1.1pp were produced. In the first cross, 150 mosquitoes were allowed to feed on mice inoculated 3 days previously with a 50:50 mixture of the two parental strains. Seven days later, a sub sample of mosquitoes (n = 25) were dissected for oocyst detection. In this case, 90% of mosquitoes were infected, with an average burden of 87 oocysts per mosquito. Given that 50% of the oocysts are expected to be the products of selfing (i.e. CU male gametes fertilizing CU female gametes, and 17X1.1pp male gametes fertilizing 17x1.1pp female gametes), and that the remaining 50% of oocysts resulting from cross-strain fertilization would each produce four recombinant progeny types, we estimate that this cross resulted in the inoculation of 15,660 recombinant progeny types to recipient mice on day 21 post-mosquito feed, when 100 mosquitoes were dissected and the sporozoites removed from the salivary glands for inoculation. For the second cross, which followed the same protocol, 60% of mosquitoes were infected with an average oocyst burden of 77 oocysts per mosquito, leading to an estimated 9240 recombinants in the cross inoculation. Selection of uncloned cross progeny for linkage group selection analysis. For immune selection, mice immunized with blood stage parasites of either P. yoelii CU or 17X1.1pp through exposure and drug cure (as above) were inoculated intravenously with 1 × 10 6 parasitized-RBC (pRBC) of the uncloned cross progeny, as described above. The resulting infections were followed by microscopic examination of thin blood smears stained with Giemsa's solution.
DNA and RNA isolation. Parental strains and growth rate-or immune-selected recombinant parasites were grown in naïve mice. Parasite-infected blood was passed through a single CF11 cellulose column to deplete host leukocytes, and the genomic DNA (gDNA) was isolated from the saponin-lysed parasite pellet using DNAzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. For RNA isolation, a schizont-enriched fraction was collected on a 50% Nycodenz solution (Sigma Aldrich) and total RNA was then isolated using TRIzol (Invitrogen).
Whole genome re-sequencing and mapping. Plasmodium yoelii genomic DNA was sequenced using paired end Illumina reads (100 bp), which are available at the European Nucleotide Archive (ENA: PRJEB15102). The paired-end Illumina data were first qualitytrimmed using Trimmomatic [36]. Illumina sequencing adaptors were then removed from the sequences. Following that, trailing bases from both the 5' and 3' ends with less than Q20 were trimmed. Lastly, reads with an average base quality of less than Q20 within a window size of four bases were discarded. Only read pairs where both reads were retained after trimming were used for mapping with BWA [37] version 0.6.1 using standard options onto the publicly available genome of P. yoelii 17X strain (May 2013 release; ftp://ftp.sanger.ac.uk/pub/ pathogens/Plasmodium/yoelii17X/version_2/May_2013/). The SAM alignment files were converted to BAM using Samtools [38]. Duplicated reads were marked and removed using Picard (http://picard.sourceforge.net).
SNP calling. The Python script [18] used to determine SNP functions as a wrapper for SAMtools mpileup and SNP calls based on mapping quality and Phred base quality scores. In this experiment the values were set at 30 for mapping quality and 20 for base quality. Also, since the P. yoelii genome is haploid and the parental strains are clonal, only SNPs where the proportion of the major non-reference allele was more than 80% were retained, to exclude possible sequencing errors or genuine but uninformative SNPs. The script produces a tab-delimited, human readable table that shows the total number of reads for each of the four possible nucleotides at each SNP. SNPS were called on both parental strains. CU SNPs were then filtered against the 17X1.1pp SNPs to remove any shared SNP calls. The remaining CU SNPs were then used as reference positions to measure the number of reads for each nucleotide in the genetic crosses produced in this study through another Python script [18]. This script produced a final table consisting of read counts for each nucleotide of the original CU SNPs in every sample.
Mathematical methods for the identification of loci under growth rate and immune selection SNP frequencies were processed to filter potential misalignment events. We note that, during the cross, a set of individual recombinant genomes are generated. Considering the individual genome g, we define the function a g (i) as being equal to 1 if the genome has the CU allele at locus i, and equal to 0 if the genome has the 17X1.1pp allele at this locus. In any subsequent population of N individuals, the allele frequency q(i) at locus i can then be expressed as for some set of values n g , where n g is the number of copies of genome g in the population. To filter the allele frequencies, we note that each function a g (i) changes only at recombination points in the genome g. As such, q(i) should change relatively smoothly with respect to i.
Using an adapted version of code developed for the inference of subclones in populations [39], we therefore modeled the reported frequencies q(i) as being (beta-binomially distributed) emissions from an underlying diffusion process (denoted by x(i)) along each chromosome, plus uniformly distributed errors, using a hidden Markov model to infer the variance of the diffusion process, the emission parameters, and an error rate. A likelihood ratio test was then applied to identify reported frequencies that were inconsistent with having been emitted from the inferred frequency x(i) at locus i relative to having been emitted from an inferred global frequency distribution fitted using the Mathematica package via Gaussian kernel estimation to the complete set of values {x(i)}; this test filters out reported frequencies potentially arising from elsewhere in the genome.
Next, the above logic was extended to filter out clonal growth. In the event that a specific genome g is highly beneficial, this genome may grow rapidly in the population, such that n g becomes large. Under such circumstances the allele frequency q(i) gains a step-like quality, mirroring the pattern of a g (i). Such steps may potentially mimic selection valleys, confounding any analysis. As such, a jump-diffusion variant of the above hidden Markov model was applied, in which the allele frequency can change either through a diffusion process or via sudden jumps in allele frequency, modeled as random emissions from a uniform distribution on the interval [0, 1]. For each interval (i,i + 1) the probability that a jump in allele frequency had occurred was estimated. Where potential jumps were identified, the allele frequency data were split, such that analyses of the allele frequencies did not span sets of alleles containing such jumps. The resulting segments of genome were then analyzed under the assumption that they were free of allele frequency change due to clonal behavior.
Inference of the presence of selected alleles was performed using a series of methods. In the absence of selection in a chromosome, the allele frequency is likely to remain relatively constant across each chromosome. A 'non-neutrality' likelihood ratio test was applied to each contiguous section of genome, calculating the likelihood difference between a model of constant frequency x(i) and the variable frequency function x(i) inferred using the jump-diffusion model. Next, an inference was made of the position of the allele potentially under selection in each region. Under the assumptions that selection acts for an allele at locus i, and that the rate of recombination is constant within a region of the genome, previous work on the evolution of cross populations [19,20] can be extended to show that the allele frequencies within that region of the genome at the time of sequencing are given by for each locus j not equal to i, where X is the CU allele frequency at the time of the cross, ρ is the local recombination rate, Δ ij is the distance between the loci i and j, x is an allele frequency, and Δ x describes the effect of selection acting upon alleles in other regions of the genome. A likelihood-based inference was used to identify the locus at which selection was most likely to act. In regions for which the 'non-neutrality' test produced a positive result for data from both replica crosses, and for which both the inferred locus under selection, and the direction of selection acting at that locus were consistent between replicas, an inference of selection was made.
For regions in which an inference of selection was made, an extended version of the above model was applied, in which the assumption of locally constant recombination rate was relaxed. Successive models, including an increasing number of step-wise changes in the recombination rate, were applied, using the Bayesian Information Criterion [40] for model selection. A model of selection at two loci within a region of the genome was also examined. Given an inference of selection, a likelihood-based model was used to derive confidence intervals for the position of the locus under selection.

Information on genes in identified loci under selection
For each combined conservative interval of relevant loci under selection, genes were listed based on the annotation available in version 6.2 of PlasmoDB and verified against the current annotation (release 26). For each gene, information on predicted transmembrane domains, signal peptides and P. falciparum orthologues. For the P. falciparum orthologues, the NS/S SNP ratios were obtained from PlasmoDB, based on the count of synonymous and non-synonymous SNPs found in 202 individual strains collected from 6 data sets stored on the website. More details on the data sets can be found at the following link: https://goo.gl/lUwKn1.

Plasmid construction to modify P. yoelii ebl gene locus
All primer sequences are given in Supplementary S7 Table. Plasmids were constructed using MultiSite Gateway cloning system (Invitrogen).
PCR amplification and sequencing of the Pyebl gene. The Pyebl gene was PCR-amplified from gDNA using KOD Plus Neo DNA polymerase (Toyobo, Japan) with specific primers designed based on the ebl sequence in PlasmoDB (PY17X_1337400). Pyebl sequences of CU and 17X1.1pp strains were determined by direct sequencing using an ABI PRISM 310 genetic analyzer (Applied Biosystems) from PCR-amplified products. Sequences were aligned using online sequence alignment software Clustal Omega (https://www.ebi.ac.uk/Tools/msa/ clustalo/) provided by EMBL-EBI.

Phenotype analysis
To assess the course of infection of wild type and transgenic parasite lines, 1 × 10 6 pRBCs were injected intravenously into five 8-week old female CBA mice for each parasite line. Since the 17X1.1p and CU-recipient strains were transfected on separate occasions, the transgenic lines were tested separately. Thin blood smears were made daily, stained with Giemsa's solution, and parasitaemias were examined microscopically.
RNA-seq. Whole blood from mice infected with P. yoelii on day 5 post-infection were host WBC depleted and saponin lysed to obtain the parasite pellet. Total RNA was extracted using TRIzol reagent. Strand-specific RNA sequencing was performed from total RNA using TruSeq Stranded mRNA Sample Prep Kit LT according to manufacturer's instructions. Libraries were sequenced on an Illumina HiSeq 2000 with paired-end 100bp read chemistry and are publicly available at the European Nucleotide Archive (ENA: PRJEB15102). RNA-seq reads were mapped onto P. yoelii 17X version 2 from GeneDB (http://www.genedb.org) using TopHat 2.0.13 [43] and visualized using Artemis genome visualization tool [44].
The prediction of disulfide bonds in our homology models were performed using DISUL-FIND (http://disulfind.dsi.unifi.it) [52][53][54][55]. Our analysis showed high probability of disulfide bond formation by this Cys351 residue. Confirming that C351 is a potential residue for forming a disulfide bond, the energy minimized stable homology models were subjected to Disulfide bond visualization to check whether the Cys351 is involved in any disulfide bond formation with any other Cys and what is the effect of the C351Y substitution.
The homology models along with their disulfide bonds were visualized (Fig 4B and 4C) and the images were obtained using the "Disulfide by Design 2.0" server (http://cptweb.cpt.wayne. edu) [56].