Identification and Functional Validation of the Novel Antimalarial Resistance Locus PF10_0355 in Plasmodium falciparum

The Plasmodium falciparum parasite's ability to adapt to environmental pressures, such as the human immune system and antimalarial drugs, makes malaria an enduring burden to public health. Understanding the genetic basis of these adaptations is critical to intervening successfully against malaria. To that end, we created a high-density genotyping array that assays over 17,000 single nucleotide polymorphisms (∼1 SNP/kb), and applied it to 57 culture-adapted parasites from three continents. We characterized genome-wide genetic diversity within and between populations and identified numerous loci with signals of natural selection, suggesting their role in recent adaptation. In addition, we performed a genome-wide association study (GWAS), searching for loci correlated with resistance to thirteen antimalarials; we detected both known and novel resistance loci, including a new halofantrine resistance locus, PF10_0355. Through functional testing we demonstrated that PF10_0355 overexpression decreases sensitivity to halofantrine, mefloquine, and lumefantrine, but not to structurally unrelated antimalarials, and that increased gene copy number mediates resistance. Our GWAS and follow-on functional validation demonstrate the potential of genome-wide studies to elucidate functionally important loci in the malaria parasite genome.


SNP Discovery
The SNP discovery methodology was similar to those described in Volkman 2007 [2]. 1X ABI shotgun sequence was obtained for nine geographically diverse parasite isolates that were previously sequenced to 0.25X coverage, bringing total coverage to 1.25X per isolate. These nine isolates include: 7G8, Santa Lucia (El Salvador), V1/S, D10, FCC-2/Hainan, D6, RO-33, Senegal V34.04 and K1. Three of the twelve previously sequenced isolates in Volkman 2007 were excluded from additional sequencing, as they were previously found to be nearly genetically identical, suggesting possible contamination in culture [2]. Reads ends with low quality (PHRED <10) bases were trimmed. Reads less than 100 bases, containing greater than 3% internal N's, or containing a mononucleotide repeat covering greater than 80% of the read were discarded. Reads were aligned to the PlasmoDB version 5 of the 3D7 genome using BLAT37 [3] requiring 95% identity, a minimum score of 100, less than 20% gaps, and coverage of at least half of the read. Only the highest scoring alignment for each read was kept and paired reads which aligned more than 10kb apart or in the wrong orientation were discarded. The Neighborhood Quality Standard (NQS) algorithm was used to distinguish real polymorphisms from sequence errors [4]. We required the SNP to have a minimum quality score of 25, and the five base neighborhood to have a minimum score of 20. We allowed one mismatch and no indels in the neighborhood. We discarded SNPs when another read from the same sample met the NQS criteria at that position but did not have a sequence difference.

Array Development and Assessment
Based on all 111,536 discovered SNPs [2,5,6] in P. falciparum, and given design parameters and unique sequence constraints, we were able to design assays for 74,656 markers. Each of 74,656 SNPs is represented by a probe set of 12 to 84 probes, for a total of 4.4 million genotyping probes on the Affymetrix 49-format array. These were hybridized to 63 unique samples (totaling 81 arrays with replicates). Genotype calls were produced using the BRLMM-P algorithm, a variant of the RLMM algorithm [7], included in Affymetrix Power Tools version 1.8.5, and clustered over all 81 arrays. BRLMM-P was forced into a haploid calling mode by setting assigning all SNPs to the "Y chromosome" and setting all arrays to "male".
The array with sample TM93C1088 is eliminated immediately after clustering (arbitrarily, since the chip claiming to be TM90C6A and the chip claiming to be TM93C1088 are identical). We also remove samples CF04.010 and Senegal Th10.04, which were suspected to be multi-clonal based upon molecular barcode analysis [8]. A halofuginone-resistant version of Dd2, a human-DNA sample, and the P. reichenowi ancestral samples are also removed at this stage, leaving 57 unique samples (totaling 75 arrays with replicates) for analysis. We then calculate a call rate for each SNP and remove 7,778 SNPs that have below an 80% call rate, leaving 66,878 SNPs. Since technical replicates showed 99.9% repeatability between chips, we merged replicate data for each of the 57 samples, producing a no-call when the replicates indicated discordant genotypes.
Concordance against sequencing data was calculated in both major and minor alleles for 17 sequenced reference strains [9]. The following 17 samples were compared against sequencing data for concordance:  [2], removing the three found to be genetically identical, and adding the two strains 3D7 and A4. A total of 18,303 SNPs lacked call overlap between array genotypes and sequencing genotypes in minor alleles and were thus removed, since concordance in both alleles could not be fully calculated. Another 30,993 SNPs were removed due to imperfect concordance, and of these discordant SNPs, most (28,789) exhibited monomorphic behavior on the array, suggesting that much of the discordance may be attributed to either a faulty assay or false discovery. The remaining 17,582 perfectly concordant SNPs constituted the high confidence set of assays used in our analyses.

Copy Number Variation (CNV) Analysis
We examined the ability to detect copy number variants (CNV) using the array by first studying a known CNV using the hybridization intensity signal of the SNP genotyping probes on the array. Kidgell, et al. (2006) [11] reported that the pfmdr1 locus was present in 3-4 copies in the Dd2 strain relative to a collection of other strains. We compared Z-scores of the normalized hybridization intensity of perfect match probes for SNPs in the neighborhood of pfmdr1 for Dd2 and six parasites estimated by Kidgell, et al. to contain only 1 copy of the locus (3D7, 7G8, HB3, D10, D6, K1). For each SNP assay we utilized the average hybridization intensity of all perfect-match probes. Hybridization intensity values were background corrected and normalized to reduce inter-array variation artifacts. SNPs with a hybridization intensity standard deviation equal to or greater than half the magnitude of the average hybridization intensity across all arrays were excluded from analysis. Figure S11 illustrates that probes for many of the SNPs assayed within the pfmdr1 locus exhibit notably higher hybridization intensity values in Dd2 relative to the other parasites, with 13 assays exhibiting average intensities greater than 2 standard deviations higher than observed in the other parasites
The following drugs were obtained from Sigma Aldrich: artemisinin, dihydroartemisinin, chloroquine, mefloquine, and quinine. The following were obtained from AK Scientific: artemether, artesunate, halofantrine, lumefantrine, and piperaquine. The following were obtained from USP: amodiaquine and atovaquone. Each drug was tested in triplicate for each parasite. Additionally, some parasites were tested with multiple biological replicates: 3D7 (nine biological replicates per drug, each in triplicate), Dd2 (three replicates) and RO-33, D10, and 207-89 (two replicates).
SNPs were filtered down to a set that contained at least 5 strains with a minor allele as well as an 80% call rate under every phenotype condition. The final data set includes 7,437 SNPs. This gives us a genome-wide significance threshold of -log 10 (P-value) > 5.17 by Bonferroni correction for multiple testing. For binary phenotype tests (Fisher's exact test, Fisher's permuted, CMH, and HLR), we used IC 50 cutoffs shown in Table S4. For tests requiring defined geographic clusters (CMH, HLR, Fisher's permuted), the three population clusters are defined by PCA, as in the LRH analysis, and the assignments are shown in Table S1.
Pointwise P-values were computed using PLINK [12]. Quantile-quantile plots (qq-plots) were used to examine the resulting P-value distributions for inflating effects due to population structure (Fig. S7). Because most of the genome is assumed to fit the null hypothesis (most of the genome should not be in association with the phenotype), significant, early deviations from this expectation may result in a high false positive rate. The null expectation is plotted as the unity diagonal line in Figure S7. Bonferroni significance is plotted as the dashed line and Benjamini-Hochberg significance is marked with the dotted line. Since most Fisher's results show evidence of inflation, we do not report these results in Figure 2 or Table 1.
Permutations of Fisher's exact test can be used to compute empirical pointwise P-values based on a simulated null distribution. We used PLINK to perform this permutation while respecting the phenotype frequencies present in our three predefined population clusters. The resulting P-value distributions (Fig.  S7) do not show inflation due to population structure, however no significant hits were found for any drug.
Similarly, the Cochran-Mantel-Haenszel (CMH) test can perform population-stratified analyses for association. We used PLINK to compute P-values (Fig. S7), and again, we see appropriate corrections for population structure, but no hits reach genome-wide (Bonferroni) significance.
The Efficient Mixed-Model Association (EMMA) test was specifically designed to handle quantitative trait associations to a data set with complex population structure using a linear mixed model [13]. It calculates a genotype similarity matrix instead of discrete categories and does not require a priori specification of population structure. The resulting P-value distributions demonstrate little remaining effect from population structure (Fig. S8) while retaining power to find a number of associations at genome-wide significance (Fig. S8, 2A, Table 1).
The Haplotype Likelihood Ratio (HLR) test is a multi-marker association test [14]. Unlike a standard, χ 2 -based multi-marker test which looks for differences in haplotype frequencies in cases vs. controls, the HLR test specifically models the likelihood that a single haplotype rose to dominance in cases while all other haplotypes proportionally decreased. It produces a LOD score, which is the maximum likelihood estimate for the haplotype frequencies observed in cases (O 1 …O k ), given the distribution in controls (f 1 …f k ): where j is the haplotype on which the mutation arose, 1-α is the recombination rate, and e j =1 when i=j and e j =0 when i≠j. The test produces maximum likelihood estimates for j and α.
So while a χ 2 -based association finds any significant differences in haplotype frequencies, the HLR test models a specific scenario that is common in rapid selective events. The HLR test does not provide significantly more power than a single marker test in regions of high LD-in extreme cases, these regions may only have two haplotypes, and a multi-marker test will have the same power as a bi-allelic SNP test. But in regions where LD is low in controls and a single, long haplotype is prevalent in cases, the HLR test is highly sensitive. The HLR test is a one-sided test and we ran separate tests for both drug resistance (called "risk") and drug sensitivity ("protect"). Results for drug sensitivity are available in Fig. S10, but are not reported generally as we are more interested in selective events for drug resistance.
We used PLINK to produce sliding window haplotypes across the genome and calculate haplotype frequencies for input to the HLR test. We produced input for all two, four and six-marker windows. The resulting LOD scores did not map well to known distributions, such as the χ 2 1-degree of freedom distribution. We instead converted the pointwise LOD scores to empirical pointwise P-values by performing approximately 370,000 permutations of the null model for each test condition. This allows us to calculate empirical P-values up to a significance of about -log 10 (P-value) = 5.6. Similar to the permuted Fisher's test, we preserved population-specific phenotype frequencies by only allowing permutations within each of our three defined populations. Resulting P-value distributions fit expectations well for the vast majority of test conditions (Fig. S9, S10) and the test demonstrates power to detect a number of loci at genome-wide significance ( Fig. 2A, Table 1). Figure M1. SNPs located in genes (28,576) were more likely to pass concordance filtering than intergenic SNPs (19,999). Figure M2. SNPs that were discovered from only one sequenced strain (35,727 SNPs) show a higher rate of monomorphism on the array than those with higher minor allele counts (12,848 SNPs). Any amount of this discordance may be explained by false discovery from sequencing data. Figure M3. To show the effect of GC composition, we took 16bp of flanking sequence on each side of the SNP to construct a 3D7-based 33mer and calculated the percent GC. The window boundaries for the graph below are chosen as the octiles of the GC distribution. SNP performance appears to worsen at GC levels below 20%, which accounts for roughly half of the SNPs. Figure M4. The effect of unique sequences in flanking regions. Although the initial design of the array excluded probes that had exact matches elsewhere in the genome, many of the remaining SNPs are in neighborhoods that contain 1 or 2 base mismatch similarity to other parts of the genome. We took 16bp of flanking sequence on each side of the SNP to construct a 3D7-based 33mer and aligned it to 3D7 using SSAHA (word length 10, step length 1, max gap 2, max insert 1, min hits 24) [15]. 28,352 SNPs aligned uniquely to their location of origin. The 20,223 SNPs that aligned in multiple locations showed a much higher rate of discordance. Figure M5. A histogram of marker spacing. Most markers are spaced closely with a few large gaps in coverage. The mean spacing of concordant markers is 1316bp with a median spacing of 444bp. Table M1. Statistics on marker spacing (gaps) by chromosome. The table gives the number of gaps, as well as the median, mean gaps, 90% percentile (90% gap), and maximum gap length (max gap) by chromosome (1 -14).