Plasmodium falciparum parasite population structure and gene flow associated to anti-malarial drugs resistance in Cambodia

Plasmodium falciparum isolates and samples size

A set of 533 blood samples collected in 2010–2011 from P. falciparum malaria patients was analysed. These samples originate from 16 health centres (11 health centres and five reference hospitals) located in 10 provinces in Cambodia. Isolates were grouped in four regions: western, southern, eastern and northern Cambodia (Additional file 1). Four control DNA samples were used to validate barcode detection (3D7, Dd2, HB3 and RO33).

DNA extraction and PCR amplification

The genomic DNA was extracted from 200 µl of blood using the DNA Mini blood kit (Qiagen, Germany) according to manufacturer’s instructions. DNA extracts were stored at −20 °C until use. DNA from reference strains 3D7, Dd2, HB3 and RO33, provided by Malaria Research and Reference Reagent Resource Center-MR4, were used as controls. Primary PCR was carried out in 25 µL of final volume with 5 µL of DNA, 0.25 µM of each corresponding primers (Additional file 2), 0.2 mM of each deoxynucleoside triphosphate (dNTP) (Solis Biodyne), 1× of reaction Buffer, 2.5 mM of MgCl2, 1.25 U FirePol^® Taq DNA Polymerase (Solis Biodyne), with the following conditions: 94 °C for 15 min, then 30 cycles of a three step program (94 °C for 30 s, 52–55 °C for 1 min and 72 °C for 1 min) and final extension at 72 °C for 10 min to reach the corresponding target between 164 and 385 bp (Additional file 3).

The nested PCR was performed in 25 μl containing 0.5 μM of each primer (Additional file 2), 0.2 mM of each deoxynucleoside triphosphate (dNTP) (Solis Biodyne), 1× of reaction Buffer, 2.5 mM of MgCl₂, 2.5 U Taq polymerase (FirePol^® DNA Polymerase, Solis Biodyne). 5 µl of the primary PCR reaction were used as the template. PCR conditions were: 94 °C for 15 min, then 40 cycles of 94 °C for 30 s, annealing temperature between 55 and 60 °C for 1 min and 72 ^°C for 1 min. A final extension at 72 °C for 10 min was performed to obtain the corresponding fragments between 100 and 200 bp. PCR of valid SNPs were performed in four multiplexed reactions (Additional file 3).

Ligation and detection assays

Nested PCR products were pooled together in two sets according to microsphere combinations. One microlitre of the pooled PCR products were used for the ligase detection reaction (LDR). The LDR was based on two allele-specific primers and one locus-specific probes (Additional file 4). The allele-specific primers were composed of two parts: the 5-prime part hybridizing with the MagPlex-Tag probe and the 3-primer part hybridizing with the PCR product. 33 different MagPlex-Tags were used to detect 40 alleles (Additional file 4) corresponding to the 20 loci which were successfully amplified by PCR (Additional file 3). Ligation was performed after hybridization of the locus-specific primer. Several MagPlex anti-Tag probes were used twice. Locus-specific probes were 5′ phosphorylated and 3′ biotinylated. LDRs were performed in a final volume of 15 µL holding in 1× of Taq Ligase buffer, 10 nM of each LDR (allele- and locus-specific primers), 4 U of Taq DNA ligase (Genesearch) and 1 µL of pooled Nested PCR. Thermocycling conditions were carried out by denaturation of the double stranded DNA at 95 °C for 1 min, followed by 32 cycles at 95 °C for 15 s and hybridization at 58.0–60 °C (Additional file 3) for 2 min. Quality control was performed by using DNA from reference strains provided by MR4. Two multiplexed reactions were used to characterize final valid SNPs (Additional file 3).

Hybridization and labeling of magnetic beads

A 5 µL fraction of the LDR product was poured into 60 µL of hybridization solution TMAC buffer (3× of tetramethylammonium chloride [TMAC] (Sigma-Aldrich), 3 mM of EDTA (Gibco), 50 mM Tris–HCl, pH 8.0 (Sigma-Aldrich), 0.1 % sodium dodecyl sulfate) and 1000 beads of each MagPlex-Tag microspheres used in the multiplex LDR, as described above. Beads quantification was performed as previously described [13]. Mixtures were heated to 95 °C for 1 min 30 s and incubated at 37 °C for 35 min to allow hybridization between SNPs-specific LDR products (Tag-probe) and bead-labeled anti-TAG probes. Then, 6 µL of 1:50 dilution of streptavidin-R-phycoerythrin (Invitrogen) in TMAC buffer was added to the post-LDR mixture and incubated at 37 °C for 20 min in 96-well plate (Eppendorf). PCR and LDR reactions were conducted in 96-well plate. The fluorescence of each allele-specific LDR products was measured on a MagPix instrument with xPonent 4.2 software (LUMINEX).

The measurement of the signal for an allele was decomposed into the signal intensity without noise and the background noise. Negative samples show reduced signal-to-noise ratio and positive samples show increased signal-to-noise ratio. The identification of negative and positive samples was based on a classification method which minimizes the variance associated to the two series of measures. This algorithm was analogous to the k-mean algorithm where k = 2. A test was used to address each measurement to the negative or positive value of the allele. Negative and positive results for the two alleles were combined to assess barcode value at considered position. Double negative results were considered as positive for the third allele at BC07 barcode position (BC07_ALT_G allele). DNA of the four reference strains 3D7, Dd2, HB3 and RO33 was introduced as the positive control on each 96-well plates. A set of 282 samples were successfully genotyped out of the 533 blood samples initially selected for the analysis, using an 11 positions barcode.

Data and statistical analysis

Comparison with the available genomic resources was performed after calling of the mutations from BAM files deposited in the ENA database for 167 samples isolated between 2007 and 2011. Correspondence of these isolates and the earlier defined subpopulations in Cambodia [8] was provided by O. Miotto. Alleles for the 11 barcode positions and k13 locus were recovered from the VCF files.

Statistical analysis was performed on the 282 samples using R software [14]. For determining the dependency of the alleles of some selected genes on the locations, firstly the association between alleles and health centres was visually addressed using correspondence analysis (Fig. 1). Analysis was performed for the 11 barcode positions. Space distribution of health centres was questioned using Between-Class analysis. The significance of differences in allele distribution was tested using Chi squared tests for independence. Loci presenting a p value <0.05 were considered as exhibiting a significantly different distribution among centres. To identify localities that account most for the SNP allele dependency on health centres, a threshold value of 1 for the Chi squared test statistics components was used (one component for each health centre, Additional file 1). Weblogos [15] were used to highlight conserved alleles among health centres.

Average F_ST (fixation index) was calculated to measure the extent of genetic differentiation within health centres. It was based on 11-SNPs barcode of 282 samples. An average allele diversity value was first calculated per sampling area (Hs) from the 11-SNPs barcode. It was then compared with the allele diversity measured for the 282 samples (Ht) to obtain the average F_ST value, ranging from 0 (no differentiation/high diversity) to 1 (complete differentiation/low diversity: subpopulations fixed for different alleles).

To define the population structure, hierarchical clustering was performed on 282 samples described by the 11-SNPs barcode. The pairwise distances between the samples were estimated as the proportion of base substitutions between them over the barcode. Ward’s minimum variance method was used to build the dendrogram (Additional file 5A). Random sampling was performed in order to obtain robust results: the clustering approach was implemented on 10,000 subsets of 230 samples each, randomly selected out of the 282 samples. Based on these 10,000 clustering results, pairwise distances between samples were calculated as the percentage of clustering results in which two samples are in the same cluster. This distance matrix was then used to build a dendrogram for all 282 samples. The number k of final clusters (conserved groups) was selected based on the dendrogram structure (Additional file 5B). The value for k = 9, 10 and 14 were producing relevant clusters, but on comparing the barcodes with the previously defined population structure [8], it was observed that the clusters produced on increasing k after 9 are not associated to a specific resistant population (KH2, KH3 and KH4) and instead are related to admixed population KHA and the mixture of other resistant subpopulations. These unspecific clusters were localized close to centre of the country. Hence, the population structure was represented by nine conserved groups (G1 to G9). Weblogo was used to highlight conserved alleles among health centres in G1 to G9 groups. The groups are geolocalized by the geographical centroids on the 2D map. The coordinates for the geographical centroids are calculated as the average of the coordinates (measured in the coordinate space of the 2D map) of the health centres included in the group. To determine the groups significantly related to geographical locations, the average distance to geographical centroid within a group was compared to the average distance to the geographical centroid when the 282 samples were randomly assigned to the nine groups. The average distance for each group was compared to the average distances calculated for 5000 random re-sampling of the 282 samples within the nine conserved groups. The p value was calculated as the proportion of average values below the average distance of the health centres to the geographical centroid for each group.

Barcode of alleles from 11 variable sites

PCR fragments were successfully amplified for 20 SNPs. There were no PCR amplification for the other four SNPs and were rejected. Assays for only 13 SNPs provide interpretable LUMINEX signal (Additional file 3) and others were excluded from the assay. Out of the 13 SNPs, the locus #7 was abandoned due to non-reproducibility and non-accuracy of the detection on control DNAs and locus #5 was rejected because genotyping analysis revealed that this locus was monomorphic. Finally, 11 SNPs were validated for barcoding (BC01 to BC11). Four multiplex PCRs and two multiplex LDRs were set up for the LUMINEX detection according to their annealing temperature.

Since the work of Daniels et al. [11], the genome version and annotation has been improved. Databases revealed that 23 SNPs amongst the 24 were located in a coding region and they are equally distributed between synonymous and non-synonymous mutations (Additional file 6). Five SNPs were located in subtelomeric regions. Genomic analysis revealed that non detection in the locus #11 was due to the presence of two nearly identical copies of the rifin gene in which the SNP is located (Additional file 6). PlasmoDB v11.1 suggested that locus #15 corresponding locus was tri-allelic. LUMINEX data treatment was adapted for this locus. Locus #24 (BC11) was validated for LUMINEX genotyping despite its low variation in Cambodian parasite population (Additional file 6). Initial analysis was performed on 533 samples, 79 were resulting from mixed infection, 183 of the samples present no significant signal with LUMINEX for at least one barcode position and 50 samples could not be amplified using PCR for at least one locus. Among the 251 rejected samples, 59 samples show more than one type of errors (Additional file 7). Finally, 282 samples among 533 blood samples were successfully genotyped at 11 SNP loci.

Allele distribution associated to health centres

Correspondence analysis was performed using barcode of all the samples (Fig. 1a). Axis 1, 2 and 3 were explaining 21, 16 and 12 % of the information respectively. For each SNP, REF (reference 3D7) and ALT (non-reference) alleles were in opposite quadrants except for BC07 barcode position which was tri-allelic. The BC11_ALT_T allele was located at the centre of the representation as it was present in nearly all samples. Despite the low number of isolates with corresponding BC11_REF_G allele, correspondence analysis showed association of this allele with eastern and southern Cambodia. Association of alleles with health centres was questioned using Between-Class analysis (Fig. 1b). Matching of relative position of health centres in the Between-Class analysis with their geographic position suggests that some alleles show association with samples geographic origin. A strict opposition between eastern and western Cambodia was observed, which can be due to specific distribution of BC03 and BC05 alleles in samples from these areas (Fig. 1). Northern and southern localities present a similar distribution pattern in other projections (Additional file 8). Present analysis suggests that allele frequencies are in agreement with the geographic location of health centres. Pursat is not in a correct position in both correspondence and Between-class analysis. It is located in the western part of Cambodia and it clusters with eastern localities. The discrepancy with its geographic localization could be due to the BC06_ALT_T and BC07 allele frequencies.

Allele frequency gradient between localities

Uneven distribution of alleles was confirmed by Chi squared analysis (p value <0.05) and the allele frequencies were represented using Weblogo (Fig. 2). Barcode BC11 was excluded from this analysis because of its low variation. Most important allele enrichments were highlighted in red and blue colors for ALT and REF alleles respectively, using the Chi squared values (Additional file 1). Allele enrichment was not restricted to a single health centre but is often present in geographically close localities. Western Cambodia, including Battambang, Pailin and Pursat provinces displayed significant enrichment of ALT alleles for BC01, BC02, BC03, BC07 and BC08 and the REF allele for BC04, BC05 and BC06 (Additional file 8). The region of Kampot (Koh Slar and Chhouk health centres) showed strong enrichment of BC02_REF_T and BC04_ALT_T with a quasi-absence of opposite alleles (Additional files 1, 9B, D). For some alleles, the gradient from ALT to REF significant allele frequency was emphasized by the presence of health centres located between ALT and REF significant geographic area where no significant enrichment could be specified for any of the alleles. For example, western Cambodia appeared as the starting point for the diffusion of BC02_ALT_C and northern Cambodia was associated with BC02_REF_T (Additional file 9B). The region of Battambang shows no significant allele frequency for this locus. An example for the West-East axis is observed for BC10 locus where the region of Pailin was associated with BC10_REF_G alleles whereas BC10_ALT_A allele was found in eastern Cambodia. Accordingly, Battambang and Pursat health centres located between these two areas show no significant bias in allele frequency (Additional file 9J). Therefore, association between barcode alleles and localities could suggest the presence of specific subpopulations with fixed allele in restricted geographic distribution and overlap between these subpopulations or even gene flow.

Presence of fixed alleles at the border of Cambodia

The presence of subpopulations was confirmed using an average F_ST value calculated per health centres. High F_ST values are observed at the localities near the borders of Cambodia (Fig. 3), including Keov Seima (eastern Cambodia). Tasanh and Sampov Loun health centres in western Cambodia are associated with high F_ST values and accordingly BC02_ALT_C allele was observed to be fixed in Tasanh region. Similarly, BC04_ALT_T and BC09_ALT_C might have contributed to high F_ST values in northern localities. The fixation of BC04_ALT_T allele was also observed in Kampot province (Chhouk HC).

F_ST analysis and gradients of allele frequencies (Additional file 9A–J) over the country suggest gene flow in a centripetal orientation. According to high F_ST values, the five locations Anlong Veng, Keov Seima, Sampoev Loun, Tasanh and Trapaing Prasat might be associated with parasite subpopulations. Crossing of subpopulations could be responsible for allele diffusion over the country. Especially in western Cambodian sites, where the low F_ST values could result from overlap between subpopulations. This likely reflects gene flow driving the homogenization of the population.

Identification of emerging subpopulations in Cambodia

Results presented in sections above suggest that subpopulations were restricted to small geographic areas. Unsupervised clustering runs based on different random subset of the 282 isolates suggested the existence of 9 robust clusters (referred as G1–G9, size of the groups, n = 18–44) representative of the parasite subpopulations. The relationship between groups and health centres was established based on distance of samples to their geographical centroid. None of the groups had samples restricted to a single health centre, and most of the geographical centroids are focused in north-west area of the country (Fig. 4).

The three groups G1, G3 and G7 were significantly associated to specific geographic area (p value <0.05). The samples in these groups were mostly isolated in western Cambodia, but also include samples from the north or from the south of the country (Additional file 10). G1 includes two samples from eastern Cambodia. G3 samples were originating from western and southern Cambodia only. G7 had two samples from the southern and one from northern Cambodia. Relationship with previously described Cambodian parasite subpopulations shows that the three groups could be associated with KH2 and KH3 subpopulations. Accordingly, samples that have been probed, carried C580Y k13 mutation. Weblogos were added to the analysis to illustrate the frequency of alleles at the 11 barcode position among conserved clusters (Fig. 4). In accordance with the results presented above (Figs. 1, 2), G1 shows conserved allele positions: BC01_T, BC04_A, BC05_G, BC08_A and BC10_G. This genotype was very close to the two barcodes associated with Pailin and Ou Chra health centres. Group G3 weblogo was more reminiscent of Promoy HC barcode (Additional file 11). The barcode analysis based on 11 SNPs was efficient to describe conserved subpopulations that emerged recently in western Cambodia concomitantly with artemisinin resistance.

The three groups G2, G4 and G8, are localized in the area between north-western region and the centre of the country (Fig. 4; Additional file 10). The average distance of these groups to the geographical centroid presents a p value between 0.1 and 0.3. The two groups, G2 (n = 32) and G4 (n = 27), have samples from various localities. Most of the samples in G2 are originating from localities in southern Cambodia and include 10 samples with C580Y k13 allele. Three barcodes in this group are found in the admixed KHA subpopulation described earlier [8]. The samples in the group G4 are mostly originating from localities in western Cambodia and include eight samples with C580Y, two samples with R539T and one sample with N458Y k13 allele. Two barcodes in this group are identified in the previously defined parasite subpopulations, one barcode is found in KH4 subpopulation (also carrying Y493H allele) and the other is found in KHA subpopulation (also carrying R539T allele). The samples in the group G8 (n = 44) are mostly originating from the localities in the north (Trapaing Prasat and Anlong Veng health centre) of the country and include four samples with C580Y and six samples with R539T k13 allele. In this group the 11 barcode loci are conserved in most of the samples and some samples have variation at BC01, BC03 and BC08 locus. Three barcodes in this group are identified in the previously defined KH3 (shown to carry R539T alleles), KH1 (ancestral population) and KHA (admixed population with C580Y alleles) subpopulations.

The three other groups G5, G6 and G9 are localized close to the centre of the country and show no significant geographical centroid p values (0.99, 0.99 and 0.44, respectively). The samples in these groups are originating from localities from all over the country. In the group G5 (n = 40) only four samples are originating from the localities in the south of Cambodia. This group includes four samples with C580Y, one sample with R539T, one sample with P553L and one sample with Y493H k13 alleles. Four barcodes of this group are identified in the KH subpopulations, two barcodes in KH1 and two barcodes in KHA. The samples in the group G6 (n = 33) are mostly coming from the localities in the southern and eastern regions of Cambodia. This group includes six samples with C580Y, one sample with I543T and one sample with V568G k13 alleles. Only one barcode of this group is identified in the KH3 subpopulation. The samples in the group G9 (n = 27) are mostly originating from the localities in eastern and western regions and 4/5 tested sample are negative for k13 allele. In this group, two barcodes are identified in the KH4 subpopulation, two barcodes in the KHA subpopulation, one barcode in KH1 subpopulation and one barcode in KH3 subpopulation. The relationship between barcodes matching the KH4 subpopulation and Y493H allele was confirmed for seven isolates of this group.

The map (Fig. 4) represents a gradient of distribution of relevant subpopulation based on barcode description from north-west to the centre of Cambodia emphasizing gene flow in that orientation. The barcodes of the groups including samples from the north-western localities are mostly associated to the KH2 and KH3 subpopulations and most of the samples carry C580Y k13 allele only. Moving towards the centre it is observed that the barcodes are associated more with the KHA, KH1 and a specific KH3 subpopulations in the north (G8). The samples are mostly carrying C580Y and R539T k13 alleles and also a rare N458Y allele in one of the samples. The groups in the centre of the country are including more barcodes associated to the KHA and KH1 subpopulations and some barcodes matching the KH3 and KH4 subpopulations. The samples are shown to carry C580Y, R539T, Y493H and the three rare mutations P553L, I543T and V568G (Additional file 10). This could suggest the localization of the admixed populations with high diversity towards the centre of the country.

Mefloquine resistance is strongly associated to northern Cambodia

Mutations in k13 gene associated with resistance to artemisinin were determined in 98 patients, as described earlier [6]. From these patients, 70 % of the samples were positive for one of the k13 resistant alleles (C580Y, R539T, Y493H, I543T, P553L, V568G & N458Y). Artemisinin resistance was more frequent in western and northern Cambodia (Chi squared test p < 0.01). The mutant alleles Y493H, I543T, P553L, V568G, and N458Y were found once in the 282 isolates. The C580Y allele was the most prevalent (54/68 positive patients) and was found to be present in all the conserved groups. Thirty-seven different barcodes were found among these 54 samples. No association was found between the C580Y allele and 11-SNPs barcode. The R539T was the second most frequent allele (9 isolates over 68 positive patients) with six isolates belonging to G8, two to G4 and one from G5. Four isolates were from northern Cambodia, three from western and two from southern Cambodia. Barcodes of these nine samples have BC01_REF_C, BC03_ALT_A, BC04_ALT_T, BC05_REF_G, BC09_ALT_C and BC11_ALT_T in common. All these alleles were significantly associated to Trapaing Prasat health centre (Fig. 2; Additional file 11).

In vitro IC₅₀ susceptibilities to chloroquine (n = 109), mefloquine (n = 111), piperaquine (n = 103) and quinine (n = 107), were assessed in isolates with a parasitaemia >0.1 % [16]. Samples were distributed among all geographical locations and clustering groups. Piperaquine showed no geographical bias. The susceptibility for chloroquine and mefloquine were lower in eastern Cambodia (Additional file 12). High mefloquine IC₅₀ values were found in isolates from Promoy, Takavit and Trapaing Prasat health centres (Fig. 5). Mefloquine resistant parasites in the region between Promoy and Takavit were mostly carrying C580Y allele. R539T mutant parasites had significantly high mefloquine IC₅₀ values (Fig. 5) suggesting two geographic loci for mefloquine resistance in Cambodia, one associated to C580Y allele and one associated to R539T allele. Large proportion of G8 samples were carrying R539T alleles and most originating from Trapaing Prasat HC in the north. The samples in G4 and G8 groups show high mefloquine IC₅₀ values (Additional file 13). The F_ST values are shown to be high for Trapaing Prasat HC (Fig. 2). These results suggest the presence of a recently emerging P. falciparum subpopulation in northern Cambodia.

Implementation of 11-SNPs barcode for mid-throughput analysis

A PCR-LDR-FMA technique for the Mid-throughput detection of a barcode for P. falciparum Cambodian isolates was optimized. This strategy included multiplex PCR and ligase detection reactions prior to hybridization with magnetic microspheres (MagPlex-Tag probes). A specific algorithm was developed for signal discrimination between the 23 alleles corresponding to 11 SNPs which enabled us to characterize 282 samples. The choice of the loci was a critical step. It was based on both possibility to run detection at the corresponding locus and allele frequency. In fact, the PCR-LDR-FMA required the design of several primers at each locus. This was the main source of rejection of SNPs. Indeed, this task was hampered by the AT richness of P. falciparum genome. The frequency of alleles in the present dataset (Additional file 6) was sometimes different from the expected non-reference allele frequency (NRAF). One explanation could be that the study of Daniels et al. included only Thai samples [11]. Discrepancy was also observed with the MalariaGEN database. Genome annotation of P. falciparum is available at PlasmoDB. Bioanalysis revealed that some of the selected SNPs were present in genes located in subtelomeric regions and encoded surface antigens (Additional file 6). These regions are known to be highly variable and to encode surface antigens. Variation in these regions might not reflect population evolution, but more parasite-host interactions such as antigenic variation.

11-SNPs barcode successfully identifies parasite subpopulations

The 11-SNPs barcode analysis at a country-wide scale confirmed that P. falciparum population is fragmented into subpopulations. Evidences were given both by individual SNPs and by global barcode analysis. The present study covers all regions where P. falciparum is endemic in Cambodia. The uneven distribution of most alleles confirmed the presence of subpopulations with restricted geographic distribution, leading us to define nine groups. Isolates from clustering groups G1, G3 and G7 were mostly located in western Cambodia. They were associated with the C580Y k13 mutant allele and correspond to the KH2 and KH3 subpopulations described earlier [6, 8]. These three groups might have restricted geographic origin as barcode allele frequencies from health centres (Fig. 2) and clustering groups (Fig. 4) were similar (Additional file 11). Subpopulation in these localities have been shown to emerge recently and correlates with artemisinin resistance specific to C580Y allele. It should be interesting to compare the barcodes in the two subpopulations (WKH-F01 and WKH-F03) described recently in the same localities in western Cambodia [9] from parasites isolated later than the present study (2011–2013).

Among the three less significant groups defined by the barcodes, G2 and G4 could not be associated to any of the specific subpopulations described earlier. The reason for G2 could be the origin of samples from the area which was not covered by genomic analysis. In fact, high F_ST value was observed in Chhouk suggesting possible emergence of subpopulation in southern Cambodia. Only five barcodes from these two groups were matching with the genotypes of defined subpopulation origin, one from KH4 and 4 from the admixed subpopulation KHA. The group G8 was associated with northern Cambodia and most of the samples were carrying the R539T k13 allele. Allele frequencies of this group were very close to that of the Trapaing Prasat HC, as illustrated by the Weblogos (Additional file 11). The high F_ST value of northern localities is in agreement with the hypothesis of gene flow associated with R539T mutation from the north. The prevalence of R539T k13 allele associated with increased ring stage survival and delayed parasite clearance rates near eastern Thailand and northern Cambodia border has been shown in recent studies [17, 18].

The 11-SNPs barcode questions the origin of subpopulations in Cambodia. Subpopulations are expected to emerge at different places and at different times over the country. Introduction of ACT is clearly associated with emergence of western parasite subpopulations. The origin of other subpopulations is not known. Barcode analysis shows that the genetic drift induced by the emergence of subpopulation disappears rapidly over time, most likely through the presence of admixed populations. Unfortunately, the present 11-SNPs barcode seems less relevant for the characterization of populations with high heterogeneity such as core population KH1 and admixed population KHA.

Relationship between subpopulations and drug susceptibilities

Groups G1, G3 and G7 were associated with resistance to artemisinin. Resistant parasites in these groups were carrying the C580Y allele only (Additional file 10). The diffusion of the C580Y allele in 2010–2011 over Cambodia was described previously [6]. The present study provides evidence that the C580Y allele can be found in all groups. In the group G8, though the C580Y allele is present, most of the samples are associated to the R539T resistant allele. Three groups contains more than one k13 allele: G4 (N458Y, 1; R539T, 2; C580Y, 8), G5 (Y493H, 1; P553L, 1; R539T, 1; C580Y, 4), G6 (I543T, 1; V568G, 1; C580Y, 6) and G8 (R539T, 6; C580Y, 4). The R539T allele was strongly associated to a subpopulation originating from northern Cambodia and represented by G8 group. The two groups G4 and G8 were associated to high mefloquine IC₅₀ values (Additional file 13). Interestingly, G8 includes highest number of samples (21/28) with associated to high mefloquine IC₅₀ values (IC₅₀ > 30 nM). Parasites were genetically close according to their barcode. Allele frequency in G8 was close from those found in Trapaing Prasat samples, suggesting that it is the region of origin of these parasites. This result is in agreement with the study by Chaorattanakawee et al., which shows the increase in occurrence of R539T allele from 2009 to 2013 and association with increased mefloquine IC₅₀ value for R539T allele (ex vivo drug susceptibility test) [19]. Despite the observed association between drug resistance (artemisinin and mefloquine) and subpopulations, the present study provides no evidence that drug pressure is responsible for emergence of subpopulations in western and northern Cambodia.

Genetic exchange between parasite subpopulations

Clustering approach provided evidence that subpopulations have emerged in different parts of the country. They might have emerged from the ancestral KH1 population described by Miotto et al. [8]. Significance of geographical centroids shows that the most recently emerged subpopulations could be well localized. The decrease of the significance could be related to the diffusion of alleles, which in that case follows the west to east major axis (Fig. 4). Nevertheless, other gene flow axis might also be present over Cambodia. This was confirmed by high F_ST values at the periphery of the country (Fig. 3). Regions with high F_ST were strongly related to forest areas which are mostly distributed at the border of Cambodia. Parasite subpopulations might have encountered independent drift of mutation and selection. Parasites moving out of their region of origin will progressively mix their genetic background with other parasites. This hypothesis is supported by the presence of G2, G4 and G8 groups containing barcodes associated to admixed population KHA. Genetic exchange between subpopulations takes place by mating. It might be supported by human socio-economical migrations. The presence of this ongoing gene flow might have supported the eastward dispersal of artemisinin resistance k13 alleles after introducing ACT in the country. Currently, environmental factors such as deforestation, development of communication axes and global welfare are changing rapidly. It will be interesting to develop a specific barcode analysis to follow the evolution of these subpopulations in this new socio-economical context.