Approval to collect and analyze clinical isolates was granted by the National Ethics and Research Committee of Côte d’Ivoire (CNER-CI) according to the protocols and standard operating procedures of Good Clinical Practices of the ICH harmonized Triplicate Guidelines for Good Clinical Practice and the Helsinki Declaration on human being research. Samples were collected following the acquisition of written informed consent from patients or their legal guardians for study participants under 18 years old.

Abidjan is an urban malaria endemic city on the southern Atlantic coast in Côte d’Ivoire, on the Gulf of Guinea, in West Africa. Malaria transmission occurs often in Abidjan, with recrudescence during the rainy season occurring from July to November.

Blood samples were collected from September 2013 to February 2014 during the rainy and dry seasons in local health facilities located across three localities in Abidjan (Figure 1) including Abobo (5°25’N, 4°1’W), Koumassi (5°18’N, 3°57’W) and Yopougon (5°19’N, 4°4’W). These samples were collected in the context of the Plasmodium Diversity Network in Africa (PDNA) project (206194, 090770). Our study involved 70 samples collected from study participants, including 31 from Abobo, 19 from Koumassi and 20 from Yopougon. Patients were eligible for recruitment if they were aged 2 years and above, presented with an axillary temperature of 37.5 °C or had a history of fever during the last 24 hours. After administrating informed consent forms, 5 mL of venous blood sample was collected in ethylenediaminetetraacetic acid vacutainers (EDTA), and leukocytes were depleted using CF11 cellulose columns and frozen at −20 °C. P. falciparum. DNA was later extracted using the QI Amp blood mini kit (Qiagen, UK), followed by the whole-genome sequencing of P. falciparum.

DNA from 70 clinical samples collected in Abidjan were prepared for library generation and whole-genome sequencing using an Illumina HiSeq platform at the Wellcome Trust Sanger Institute, UK, as part of the MalariaGen P. falciparum Community project [21]. Standard laboratory protocols were used to determine DNA quantity and the proportion of human DNA in each sample. Samples were put forward for whole-genome Illumina paired-end sequencing [14,22]. All the 70 samples passed the whole-genome sequencing quality standards. Sequenced reads were mapped to the P. falciparum 3D7 reference genome using bwa mem, version 0.7.15, with -M parameter to mark shorter split hits as secondary [23]. Standard alignment metrics were generated for each sample using the stats utility from sam tools, version 1.2 [24]. We also used GATK’s CallableLoci (version 3.5) to determine the proportion of genomic positions that could be identified in each sample [25]. Potential SNPs (Single-Nucleotide Polymorphisms) and indels were discovered by running GATK’s HaplotypeCaller (version 3.6) independently for each of the 70 sample-level bam files. This resulted in genotype calls for both SNPs and short indels. SNPs and indels were filtered separately. Each variant was assigned a quality score using GATK’s Variant Quality Score Recalibration (VQSR), version 3.6. The tools, VariantRecalibrator and ApplyRecalibration, were used for this purpose. Regions of the genome that we previously identified as being enriched for errors were masked out. VariantsRecalibrator was employed using the PASS variants from P. falciparum crosses as a training set, 1.0 release. For SNPs, we previously used 15.0 for the training set variants. ApplyRecalibration was then used to assign each variant a quality score named VQSLOD (Variant quality score). Higher values of VQSLOD indicate a higher quality. Variants (both SNPs and indels) with a VQSLOD score under or equal to 0 were filtered out. Variants were excluded from the analysis if they were positioned within subtelomeric regions located within the hypervariable Var, Rifin, and Stevor gene families, or were positioned within repetitive sequences as identified using Tandem Repeat Finder. Data were then filtered out to extract 89578 bi-allelic SNPs in the core genome with a VQSLOD score greater than 6. In addition, we filtered out those SNPs to exclude isolates with missing calls greater than 10% in all positions and SNPs with calls missing in more than 10% of isolates. A total of 89211 SNPs remained after filtering.

Analysis of allele frequency distributions, including within-population Tajima’s D test, was performed using Vcftools and custom R scripts to identify genes under balancing selection [26]. For Tajima’s D test, we extracted bi-allelic SNPs that were segregating within our population. Missing data were excluded via the removal of 8 individual isolates on a gene-by-gene basis due to the observation that most of missing data clustered within a small number of isolates. The allele frequency spectrum for each gene was assessed with at least 3 SNPs using custom R scripts. Typable SNPs in this study were classified as synonymous or non-synonymous based on amino acids changes and compared to the 3D7 reference genome sequence. The ancestral state of SNPs was determined by comparing each sequence to homologous sequences in Plasmodium reichenowi (Pr), a parasite with recent common ancestry.

To detect loci under recent positive selection, we computed the standardized integrated Haplotype Score (| iHS |) for each SNP with no missing data and a minimum minor allele frequency of 0.05 using the REHH R software package. iHS was calculated for each SNP with no missing data and a minor allele frequency greater than 0.05 [27]. The genetic map distance between markers inferred with LDhat 2.2 [28] was measured using a block penalty of 10 million rjMCMC iterations and a burn-in of 100 000 iterations. Selection windows were defined by calculating the distance required for the extended haplotype homozygosity of each SNP to decay to a level of 0.05 in each direction. Overlapping EHH windows from individual high-scoring SNPs (|iHS| >3.29, suggestive line) were combined into continuous windows, and windows supported by only a single SNP position were subsequently discarded. Bonferoni correction was applied for genome-wide significance. Reference and non-reference alleles were described as ancestral and derived alleles, respectively. For significance, the REHH package generated a two-sided p-value where Ф (iHS) represents the Gaussian cumulative distribution function.

Infection complexity was determined using the Fws fixation index [14,29] using the moimix package in R and computed with R scripts. For all bi-allelic-coding SNPs, Fws was calculated using the formula, Fws = 1 − (Hw /Hs), where Hw is the within-individual heterozygosity and Hs is the within-population heterozygosity. For each bi-allelic SNP, heterozygosity was estimated using the formula, H = 1 − (p2 − q2), where p and q are the frequencies of the two alleles (p = 1 − q). At each SNP, p and q were estimated for each sequence as the proportions of sequencing reads that carried each allele in the individual sample. At the population level, the allele frequencies at the SNPs level were estimated as the mean of the allele frequencies in the individual composing the population sample. The minor allele frequency (MAF) was reported as the frequency of the least common allele for that SNP. Within-individual and within-population heterozygosities were computed and assigned to ten equal-sized MAF bins from [0.0–0.05] to [0.45–0.5], and for each bin, the mean within-sample and population heterozygosities were computed.

We conducted PCA and Fst analysis, as implemented using PLINK1.9 and Vcftools, respectively, to assess the structure of study populations. We used Weir and Cockerham’s population genetic differentiation estimator, Fst, to study the sites [30]. For PCA, we applied the Linkage Disequilibrium correction to remove correlated pairs of SNPs and identity by descent (IBD) to identify and remove any closely related samples before computing the principal components (PCs). We further calculated the top 10 eigenvectors from the population genotype. R statistical package was used to analyze our data. For Fst analysis, missing data for some isolates were excluded on a per SNP basis.

There was an excess of rare variants in the population, with the majority (80%) being single isolates (Figure 2A). The coding sequences had higher coverage than the intergenic regions did, probably due to A + T allelic richness. There was an excess of non-synonymous SNPs compared to the number of synonymous SNPs, showing a clear mutation process in the isolates (Figure 2B).

Genes with at least three SNPs had their allele frequency distribution assessed via computing a Tajima’s D statistic based on 16697 SNPs. For the 262 genes analyzed, Tajima’s D values were mostly negative (mean = −0.15), with only 47 genes (18%) having positive Tajima’s D values, indicating an excessively low frequency compared with that which was expected for a mutation-drift equilibrium population (Figures S1 and S2). At the whole-genome level, we identified 18 genes with Tajima’s D values greater than one (Table S1). PfEMP1, SURFIN 4.2, GLURP, and msp7-lik were the most significantly represented genes.

We conducted a scan of the whole genome of P. falciparum using the REHH package in R. Via an assessment of the standardized integrated Haplotype score (| iHS |), we identified a recent positive selection in seven loci that had four or more SNPs with a standardized | iHS | value greater than 3.2 (suggestive line) (Figure 3A) and four loci that had at least one SNP with an | iHS | greater than 5 (Figure 3B and Table 1). We identified windows containing genes that have been under strong and recent positive selection. Indeed, there was a strong signature surrounding one of the major chloroquine resistance gene, Pfcrt on chromosome 7. In addition, a strong signature was also observed around genes encoding erythrocyte membrane protein 1 (PfEMP1), a gene involved in the adhesion of P. falciparum to endothelial cells. Pfdhps, a drug-resistant gene located on chromosome 8, was also found to be under positive selection. No detectable signature was found around multi-drug-resistant and anti-folate drug target genes, mdr1 and dhfr, respectively. Another region with a very high | iHS | value on chromosome 11 had SNPs within the ama1 gene. Windows spanning regions of elevated | iHS | values were also observed on chromosome 10 and 13, covering genes encoding for the merozoite surface protein 7 like antigen, Pftrap (Thrombospodin-related adhesive protein), GLURP (Glutamate-Rich Protein) and msp3 (merozoite surface protein 3).

The within-infection Fws fixation index describes the relationship between the diversity observed within a patient to that of the population using estimates of heterozygosity. It provides a measure of the risk of outcrossing between parasites within an individual to generate new genotypes during recombination in a mosquito host [29]. In our analysis, the Fws scores ranged from 0.50 to 0.99 (mean 0.89; median 0.98), with 64% samples presenting high Fws estimates (i.e., >0.95; Figure S3).

We conducted principal component analysis using 89211 SNPs with no missing data in all 70 clinical samples. The first three principal components (10.8% of the total variation) were plotted and showed no evidence of population structure in most samples from the three populations. Only a few isolates appeared as outliers that were not very divergent (Figure 4).

We measured Fst to assess the genetic differentiation and evaluate the overall effect of population substructure of parasite populations. There was a minimal amount of differentiation between the populations highlighted by a very low value for the mean Fst estimate (mean 0.001) (Figure 5). Only 12 SNPs had Fst values equal or greater to 0.2, and the highest Fst value was 0.7 (Table 2).

Three of the differentiated SNPs were located on chromosome 7. The SNP at position 435368 spanned the region coding for P. falciparum esterase gene Pfpare (P. falciparum prodrug activation and resistance esterase), which is responsible for a resistance to pepstatin esters [31]. Furthermore, an SNP located on chromosome 8 encodes an amino acid within the SURFIN_8.2 gene (PF3D7_0830800), a polymorphic antigen that is expressed on the surface of P. falciparum-infected erythrocytes (IE) and released merozoites [32]. The PF3D7_1035200 gene containing two SNPs with highest genome-wide Fst values were located within a single region on chromosome 10, encoding a conserved protein with an unknown function.