Micro-epidemiological structuring of Plasmodium falciparum parasite populations in regions with varying transmission intensities in Africa

Background: The first models of malaria transmission assumed a completely mixed and homogeneous population of parasites. Recent models include spatial heterogeneity and variably mixed populations. However, there are few empiric estimates of parasite mixing with which to parametize such models. Methods: Here we genotype 276 single nucleotide polymorphisms (SNPs) in 5199 P. falciparum isolates from two Kenyan sites (Kilifi county and Rachuonyo South district) and one Gambian site (Kombo coastal districts) to determine the spatio-temporal extent of parasite mixing, and use Principal Component Analysis (PCA) and linear regression to examine the relationship between genetic relatedness and distance in space and time for parasite pairs. Results: Using 107, 177 and 82 SNPs that were successfully genotyped in 133, 1602, and 1034 parasite isolates from The Gambia, Kilifi and Rachuonyo South district, respectively, we show that there are no discrete geographically restricted parasite sub-populations, but instead we see a diffuse spatio-temporal structure to parasite genotypes. Genetic relatedness of sample pairs is predicted by relatedness in space and time. Conclusions: Our findings suggest that targeted malaria control will benefit the surrounding community, but unfortunately also that emerging drug resistance will spread rapidly through the population.


Introduction
The earliest models of malaria transmission assumed a completely mixed and homogenous parasite population 1,2 . However, malaria transmission is highly heterogeneous, and follows the Pareto principle where 80% of infections occur in only about 20% of the population 3 . Consequently, there is increasing interest in models allowing for spatial heterogeneity and variably mixed populations of parasites 4-7 . There are now several epidemiological studies describing spatial heterogeneity of malaria on varying geographical scales [8][9][10][11][12][13][14][15][16][17][18][19] . This heterogeneity is characterized by infection hotspots which usually persist even after transmission has been reduced in surrounding areas 9,11,20-25 , and thus act as reservoirs of infection 21,26 . Achieving any meaningful reduction in transmission in regions containing malaria hotspots will require a scale up of control activities, including repeated mass administration of Artemisinin Combination Therapy (ACT) drugs, increased coverage of long lasting insecticide treated nets (LLINs) and intensive indoor residual spraying (IRS). These measures are very costly and may not be realistic for universal coverage in most of the resource-poor endemic countries. Thus, targeted control may be more important, and is likely to be required to eliminate malaria 3,21,27,28 .
Mathematical models show that targeting hotspots may reduce transmission in surrounding areas 11,22 . These models, however, assume that hotspots are stable and that mosquito mixing in the community is homogeneous 22 . Studies have shown that certain species of mosquitoes exhibit some level of site fidelity, where they return to the same homesteads to feed 29 . If such behaviour is the norm with very little mixing, then this would greatly reduce the community-wide impact of targeted interventions, and interventions would be beneficial only to individuals within the targeted region. If, however, transmission networks operate freely over large geographical areas, then these interventions would likely have an impact beyond the targeted region. Furthermore, parasite evolution takes place in a micro-epidemiological context and the spread of drug resistance or new antigenic variants through the population will also be critically dependent on the degree of mixing of parasite populations.
Few studies currently provide empiric evidence on the mixing of parasites over space and time, yet this evidence is important as parasite mixing is likely to affect the outcome of targeted control interventions 23 . The community-wide impact of targeted control has not been studied extensively, although early controlled trials showed that bed nets were effective at reducing child morbidity and mortality associated with malaria, in villages or communities randomised to the intervention in The Gambia 30 and Kilifi 31 . More recent studies have shown that the use of bed nets in a village randomized to intervention in Asembo, western Kenya, also protected individuals just outside the intervention village who were themselves not using bed nets 32 . A cluster-randomized controlled trial on the impact of targeting integrated control measures to hotspots showed temporally limited effect on reducing transmission in areas surrounding the targeted hotspots 23 . In order to inform future targeted control strategies more precise empiric data on parasite mixing is required.
We hypothesized that by genotyping parasites with fine-scale temporal and spatial data we would be able to determine fine-scale structure to the population and infer the degree of parasite mixing over small geographical areas which are likely to be the focus of targeted malaria control programs 23,27 . We used SNP genotyping of Plasmodium falciparum field isolates from three African sites and analysed the genetic relatedness among parasites within individual sites, in order to determine the level of parasite mixing on micro-epidemiological scales in each population. Principal Component Analysis (PCA) was used to detect parasite subpopulations in each site, and tests of spatial autocorrelation including Moran's I and spatial scan statistics were used to test for autocorrelation among parasite genotypes. The analyses were carried out at different spatial scales ranging from intensive withinvillage surveillance through to county-wide surveillance.

Materials and methods
Study sites P. falciparum infected blood samples were collected from individuals at three sites in two African countries: Kombo coastal districts of The Gambia on the West African coast; Kilifi, Kenya on the East African coast, and Rachuonyo South District in the Western Kenyan highlands. The Gambia has a subtropical climate with a single rainy season between the months of June and October 33,34 , while Kenya has two rainy seasons, experiencing short rains between October and December and long rains between April and August 35 . In all three sites, P. falciparum is the main causative agent of malaria 22,33,35 and transmission occurs almost exclusively during and immediately after the rainy seasons 34,36 . The common vectors in The Gambia are Anopheles gambiae s.s., Anopheles arabiensis and Anopheles melas 37 , while the common vectors in the Kenyan coast have historically been A. gambiae s.s. and A. funestus, but a recent shift to A. arabiensis and A. merus has been detected along the coast 38 . In Rachuonyo South district, the main vectors transmitting malaria are A. gambiae s.l. and A. funestus 39 . Temporal trends show

Amendments from Version 1
In response to the reviewer's comments, we have included information on the changes in malaria transmission intensity in each study site during the study period. We have also modified Figure 5 to show the 95% confidence intervals around parasites collected one day apart, and added more information on how this figure was generated. We have additionally provided information on the analysis of SNP subsets, specifically those SNPs typed in EBA175 and AMA1. Genes typed included antigen-encoding, housekeeping and hypothetical genes. 52 and 9 SNPs were typed in the antigen-encoding parasite ligands Erythrocyte Binding Antigen 175 (EBA-175) and Apical Membrane Antigen 1 (AMA-1), respectively. In the Kilifi parasite population, between 158 and 226 SNPs were typed in each sample, while in The Gambia and Rachuonyo south populations, 131 and 111 SNPs were typed in 143 and 2744 samples, respectively. Genotyping was done on the Sequenom MassARRAY iPLEX platform, which allows multiplexing of up to 40 SNPs in a single reaction well and differentiates alleles based on variations in their mass 51 . Locus specific PCR and iPLEX extension primers were designed with the sequenom MassAR-RAY designer software (Version 3.1) using 3D7 as the reference genome (PlasmoDB release 9.0) (Dataset 2 52 ). A multiplexed PCR reaction was performed by pooling locus-specific primers, and un-incorporated dNTPs were dephosphorylated enzymatically using shrimp alkaline phosphatase. Extension primers binding immediately adjacent to the SNP site of interest were then extended by a single nucleotide base, using mass-modified dideoxynucleotides. The extended products were resin cleaned to remove excess salts and the mass of the different alleles determined using MALDI-TOF mass spectrometry.
Sample and SNP cut-off selection criteria Genotype data was aggregated to determine genotyping success rates for individual samples and SNPs. Samples where >40% of SNP typing failed were excluded from analysis, and among the remaining samples, SNP typing for which >30% of samples failed were further excluded from analysis. The criteria for successful SNP typing were based on the SNP intensity values (r) and allelic intensity ratios (theta). Alleles were called as successful if they were above an intensity cut-off value ranging between 0.5 and 1.0, set depending on the performance of the individual SNP assay, and were classified as failed if they were below this cut-off. For those SNPs that were above the cut-off, allelic intensity ratios ranging between 0 and 1 were used to classify them as homozygous (single parasite genotype infections) or heterozygous (mixed parasite genotype infections

Statistical analyses
All statistical analyses of genotype data were conducted in R statistical software (version 3.0.2) 53 except for the spatial scan statistics which were computed using SaTScan software (version 9.3) 54 . Analyses were carried out separately for each parasite population, except for the Fixation index (FST) analyses which by definition involve the comparison of populations and so were carried out between samples in the different sites.
In each population, genotype data for all samples was aggregated and analysed collectively. Separate analyses were also carried out for subsets of SNPs typed in EBA 175 (39, 36 and 20 SNPs in The Gambia, Kilifi and Rachuonyo South, respectively) and AMA1 (9 SNPs in The Gambia and 8 SNPs in Kilifi).
Only 3 SNPs were genotyped in Rachuonyo South, so this SNP subset was not analysed separately. In the Kilifi population, we ran additional analyses for samples collected from community surveys (asymptomatic infections) and hospital admissions (symptomatic infections).
Calculating pairwise time, distance and SNP differences. Analyses were carried out separately for each of the three sites. Each parasite was compared to every other parasite in that site (i.e. a pairwise analysis), noting the time, distance and SNP differences between the parasite pair (Dataset 3 -Dataset 5 55-57 ). We took half the lower limit of detection of temporal and spatial differences for parasite pairs collected on the same day and/or at the same location. Parasite pairs collected on the same day were assigned a difference of 0.5 days. For older samples in Kilifi (i.e. collected prior to 2004) where location was known to a 5 km accuracy, pairs collected at the same location were assigned a difference of 2.5km. We had precise geospatial co-ordinates for recent samples in Kilifi (i.e. collected after 2004) as well as all samples from The Gambia and Rachuonyo South, so parasite pairs in these three groups collected from the same location were assigned a difference of 0.02km. SNP differences were computed by comparing genotype data for parasite pairs within each population and counting the number of SNPs between them. Missing SNP data for each parasite was replaced with the major allele in the respective population, after excluding SNP typing where >30% of assays failed as described above.

Population genetics analyses.
Minor allele frequencies were computed for SNPs in each population. Principal components analysis (PCA) was performed using singular value decomposition on a covariance matrix of pairwise SNP differences between parasites in individual populations. To detect inter-population genetic differentiation and within-population genetic diversity, we restricted analysis to 33 SNPs that had been successfully typed in all three populations.

Spatial autocorrelation.
Moran's I was calculated using geographical coordinates to specify location and scores for the first 3 principal components to specify associated attribute values. Moran's I was computed at distance classes of 1 km, 2 km and 5 km, using 100 bootstrap resampling steps to determine statistical significance.
Spatial scan statistics were calculated using SaTScan software and were run separately for each study site. The statistics involved running a purely spatial, retrospective analysis based on a normal probability distribution model using continuous variables (PC scores) and looking for areas with clusters of high PC scores. Latitude and longitude coordinates were used to represent the geographical locations of specific parasites, whereas principal component scores were used to represent individual parasite genotypes. During the analysis, a scanning window that gradually varies in size from including only a single homestead up to 50% of the population moves over the geographical space and at each window size and location, the ratio of parasites with high PCs inside the window versus outside the window is calculated. The window with the highest ratio is noted down as a cluster and its statistical significance is determined after accounting for multiple comparisons using random permutations.

Raster analysis.
To identify possible spatial barriers to parasite movement and mixing over short distances, each study area was divided into pixels of varying sizes which were then scored with 1 or 0, based on whether or not a straight line linking any two parasites crossed their boundaries. These pixels were then used as independent variables in a multivariable linear regression analysis that had the number of SNP differences as the dependent variable. Significance of the coefficient estimates were determined using non-parametric bootstrapping with 100 resampling steps.
To test for correlations between transmission intensity and population genetics at fine scale, each pixel was assigned the mean of the PC scores and either Malaria Positive Fraction (for Kilifi data) or asymptomatic parasite prevalence by PCR (for Rachuonyo) for all samples found within that pixel. The correlation between PC score and MPF or between PC score and parasite prevalence was tested by Spearman's rank ordered correlation coefficient.

Results
Study populations 5199 P. falciparum parasite isolates were collected from the Kombo coastal districts in The Gambia, and Kilifi County and Rachuonyo South district in Kenya ( Figure 1) between 1998 and 2011. 107, 177 and 82 SNPs were successfully genotyped in 133, 1602, and 1034 parasite isolates from The Gambia, Kilifi and Rachuonyo South district, respectively (Table 1). 26, 57 and 49 SNPs were present at frequencies of 5% and above in The Gambia, Kilifi and Rachuonyo, respectively. In each of the populations, there was a positive correlation between SNP assay performance and parasite density.
In all study sites, separate analyses of EBA175 and AMA1 did not reveal qualitatively different results from the pooled analyses Analysis of within-population genetic diversity (π), based on a set of 33 SNPs that had been typed in samples from all three populations, showed that parasites in Rachuonyo South had the highest genetic diversity with an average of 3.384 (95% CI: 3.380 -3.388) SNP differences per parasite pair. Those  .3%, PC3=4%) of the variability seen in Rachuonyo South. We were unable to resolve parasite populations into distinct sub-populations using principal component analysis ( Figure 2 and Figure 3, Supplementary Figure 4).

Global and local spatial autocorrelation analysis
Having not seen sub-populations by PCA alone, we then included spatial analyses to test for spatial structure to the principal component values. Moran's I analysis for spatial autocorrelation showed slight positive correlations for parasites that were statistically significant for at least one principal component at 2 km and below in The Gambia, 5 km and below in Kilifi, and 1 km and below in Rachuonyo South ( Figure 4).
Spatial scan statistics using SaTScan identified statistically significant (p≤0.01) clusters of distinct parasite sub-populations of different sizes in Kilifi and Rachuonyo South. In Kilifi, one cluster with a radius of 1.54 km (p=0.01) containing 15 parasites was detected, while in Rachuonyo South, a smaller cluster of genetically distinct parasites was detected with a radius of 0.5 km (p=0.001) containing 14 parasites. No clusters were detected in The Gambian population, indicating that parasites did not group into distinct sub-populations in this study site.

Spatio-temporal variations in genetic differences between parasite isolates
We examined the effect of distance and time separating parasite pairs on genetic relatedness to determine the spatial extent and rate of parasite mixing. We used linear regression models where the number of SNP differences between parasite pairs was an outcome predicted by the distance between parasite pairs and the time between parasite pairs. Time was not included for the Rachuonyo South population as the samples were collected in a single cross-sectional survey taken over a few days. Across all three datasets, distance was independently associated with increasing variation in genotype, i.e. the further apart in space any two parasites were, the greater the number of SNP differences between them. In The Gambia and Kilifi populations, time was also shown to be associated with increasing variation in genotype, with parasite pairs collected further apart in time having greater number of genetic differences. Additionally, in The Gambia and Kilifi populations, time interacted antagonistically with distance to attenuate the effect of distance on genotype relatedness ( Figure 5). This means that the genetic differences between any two parasites increased with distance, but at a decreasing rate when time between these samples increased. We observed that in The Gambian population, parasites acquired SNP differences over distance at a slower rate than in the Kilifi and Rachuonyo populations.
Bootstrapping the analyses (to take into account the linked nature of pairwise observations) gave statistically significant effects of distance, time and the interaction between distance and time ( Table 2).

Identification of geographical barriers to parasite movement
We conducted raster analysis by pixels to examine a) the spatial relationship between distinct parasite genotypes as represented by the principal component analysis and either malaria positive fraction (MPF) data (in Kilifi) or PCR positive data (in Rachuonyo South) and b) the presence of possible spatial barriers to parasite movement that would act as factors. "The range of MPF and parasite prevalence per pixel varied depending on the size of the pixels analysed. In the Kilifi population, MPF ranged from 0 -100% (interquartile range (IQR) = 20%) for the 0.5km pixels; 0 -100% (IQR = 14%) for the 1.0km pixels; 20 -83% (IQR = 7%) for the 2km pixels; and 33 -63% (IQR = 4.7%) for the 4km pixels. In the Rachuonyo South population, PCR positive prevalence varied from 0 -75% (IQR = 19.4%) for the 0.5km pixels; 0 -47% (IQR = 17.4%) for the 1.0km pixels; 3.5 -35.8% (IQR = 14.3%) for the 2km pixels; and 6.2 -33.4% (IQR = 7.8%) for the 4km pixels." The analysis of principal components did not show any consistent or statistically strong associations with markers of transmission intensity (i.e. malaria positive fraction and prevalence of asymptomatic parasitaemia by PCR) (Supplementary Figure 5).
Bootstrapping the multivariable linear regression analysis of pairwise comparisons of samples for SNP differences using 189, 703 and 340 pixels for The Gambia, Kilifi and Rachuonyo South, respectively, showed that the majority of pixels were not significant influences on SNP differences (Supplementary Figure 6). The few pixels that were significant (p<0.05) were non-significant after applying Bonferroni correction to account for multiple testing. Furthermore the distribution of p values was uniform for each dataset (mean p value ~0.5 in each population).

Discussion
As malaria transmission declines, targeted control at the microepidemiological scale is likely to be important in eliminating malaria in any remaining transmission foci. The effectiveness of such targeted measures will depend on the extent of parasite mixing in and around these foci 23 . In the current analysis, we did not identify any population structure by simple inspection of the Principal components derived from SNP genotyping in The Gambia, Kilifi and Rachuonyo South (Figure 2 and Figure 3), indicative of a parasite population that is well mixed. However we did not conclude that there was no structure to the population, only that we could not identify it in the absence of spatial data.
We therefore went on to analyse the genotype data using spatiotemporal data, and identified spatial autocorrelation using Moran's I in all three populations, with statistical significance (p<0.01) for the first principal component in The Gambia and Kilifi and the third principal component in Rachuonyo South (Figure 4). Overall, the consistent pattern observed in the Moran's I analyses was that of spatial auto-correlation at close proximity (i.e. at a range of a few km), and little or no auto-correlation at larger distances. The auto-correlation was modest in effect size but statistically significant with p values ranging from 0.01 to 0.001 at < 1 km. However, using scan statistics we identified only two specific clusters of distinct parasite sub-populations based on PC scores, one in Kilifi and another in Rachuonyo South. The limited evidence of specific local clusters of parasite sub-populations in the face of evidence of spatial auto-correlation over the whole study site implies that there is a high degree of mixing among parasites within the study sites, leading to limited clustering of parasites into genetically distinct sub-populations.
We further looked at the effect of time, distance and time-distance interaction on the variation in SNP differences between parasite pairs within individual study sites. Since the number of days differed for almost all parasite pairs, dummy data were included in the regression analysis to enable the generation of time-distance interaction graphs. For each study site, a distance range of 1 -10km (with an interval of 0.1km between adjacent distances) was used. Temporal distance with 14 and 10 day intervals were assigned to parasite pairs in Kilifi and the Gambia, respectively, whereas time was not considered for the Rachuonyo South population. Constant SNP differences of 14, 10 and 8 were used for parasite pairs in Kilifi, Rachuonyo South and the Gambia, respectively. We found that time and distance were independently associated with increasing variation between parasite genotypes (i.e. the further apart in time or space two parasites were, the greater the genetic differences observed between them). However, in the case of The Gambia and Kilifi populations where we had longitudinal data, time was shown to interact antagonistically with distance, with an increase in time reducing the variations in genetic differences between parasites as distance between the parasites increased ( Figure 5). This implies that distance between samples was no longer predictive of genetic variation when there were longer time periods between samples, indicating that, given enough time, even parasites that are separated by large distances would get a chance to interact and recombine, especially if they are not geographically isolated. The number of SNP differences were seen to plateau at approximately 1km in the Gambia, 3km in Kilifi and 10km in Rachuonyo South. This may be attributed to the characteristics of the local parasite population, which in turn may be explained by the distribution of human settlement in the areas sampled, for example in the Gambia, homesteads tend to be clustered together in distinct, autonomous villages whereas in Rachuonyo South there is a denser and more uniform pattern of human settlement over the study area, enabling the interaction of parasites over a much larger distance.
Lack of genetic structuring of parasite populations observed in this study is indicative of a population that is well mixed. This observation of a highly mixing parasite population is in agreement with results of similar studies of parasite genetic diversity and population differentiation using microsatellites 58,60,61 , immune selected genes 62,63 and SNPs 64 . These studies were carried out in parasite populations from different geographical regions representing a diverse range of transmission intensities, from the highest in Africa and oceania, intermediate in Southeast Asia, to the lowest in south America. However, other studies have shown population structure when looking at the same population 50,65-67 , although these analyses were carried out on larger geographical scales than those analysed here and mostly involved analyses at provincial, country or continental levels. Population structure was most evident in regions with low transmissiion intensities such as south America or southeast Asia 58 , and less evident in Africa where transmission intensity is much higher 61 .
On an international level, for example, some studies have been able to distinguish between Senegalese and Thai parasite isolates using a 24-SNP barcode 68 , and another study using 4 SNPs out of a set of 384 SNPs was able to resolve East and West African parasites 50 , showing that parasite populations can be resolved on a large geographical scale. A study in Senegal was also able to identify population structure among parasites using a 24 SNP barcode, despite a high level of similarity among the parasites analysed 69 . It is possible that more detailed genotyping using a larger number of markers, for instance by whole genome sequencing, would start to identify mutations that are private to particular sub-populations at a finer geographical scale, although the degree of mixing observed here suggests that discrete populations are unlikely.
We identified spatial autocorrelation among parasites in the different study areas. However, most of these correlations were found over short distances, pointing to the existence of parasite subpopulations over small spatial scales. This indicates the presence of clusters of genetically distinct parasites at micro-epidemiological scales within the study sites. Previous studies have identified parasite sub-populations based on clustering of serological responses to the important antigen Plasmodium falciparum Erythrocyte Membrane Protein 1 (PfEMP1) in children in Kilifi 70 , supporting our observations of parasite sub-populations at this site. In Papua New Guinea, sub-populations of parasites have also been identified at a micro-epidemiological scale using PfEMP1 71 , indicating that this may be a good marker for population differentiation at the micro-epidemiological level.
Studies on hotspots of symptomatic malaria infection have identified hotspots or clusters of infections down to the level of individual homesteads in Kilifi 9 . The lack of consistent correlations between parasite genotypes and infection prevalence shown through raster analysis of pixels in this study (Supplementary Figure 5) indicate that infections within higher incidence areas are likely not caused by distinct parasite sub-populations. Instead, such infections are likely caused by parasites that are well mixed within the general population. Our inability to detect barriers to parasite movement over short distances indicates that parasites move freely within the study areas, and the spatial extent of such parasites may be limited only by the ecology and dispersal range of mosquito vectors. Furthermore, recent examination of the epidemiology of hotspots shows that they occur at the full range of spatial scales, with a pattern of spatial auto-correlation that does not show a discontinuity at any scale (i.e. a smooth semi-variogram) 9 . This further argues against the existence of discrete "units" of transmission with sub-populations of parasites.
This has implications for public health interventions that may target transmission hotspots. If hotspots consist of distinct parasite populations that do not mix with parasite populations in the wider parasite community, the impact of hotspot-targeted interventions beyond the hotspot boundaries can be expected to be limited. If parasites mix freely, as suggested by our data, the impact of hotspot-targeted interventions may affect community-wide malaria transmission. This assumes that hotspots can be detected, are stable in time 20 and the spread of parasite populations indeed primarily occurs from hotspots to the surrounding community 23 .
This study had some limitations. First, the number of SNPs typed was relatively small, and this would have limited our power to detect genetic structuring among the highly similar parasite populations, especially in The Gambia. Detecting structuring in highly similar parasite populations may require either a much larger panel of SNPs or the use of more informative SNPs, as shown in the study by Campino et al, 2011 50 . Advances in sequencing technologies have increased the use of whole genome sequence data in the analysis of P. falciparum parasite population genetics, and this has led to the identification of hundreds of thousands of SNPs, most of which are present at very low frequencies especially in African parasite populations 72 . Additional analyses will require the use of whole genome sequence data to identify rare variants and distinguish between closely related parasites, thus allowing parasite population structure to be analysed at fine spatial scales. However, despite the small SNP panel used in this study, we were still able to detect population structuring on a micro-epidemiological scale. Our analysis suggests that this structure was a uniform spatial and temporal auto-correlation rather than driven by discrete clusters of parasites at specific locations. Despite the limitations of our SNP typing and sample size we can therefore conclude that any specific clustering is less prominent as a feature than the auto-correlations in space and time that we can detect.
A second limitation is that we conducted our study in only two sites in Kenya, and one site in the Gambia. It may be premature to generalize our results more widely and an analysis of more sites will be required to make confident generalizations. On the other hand the three sites selected do demonstrate differing transmission intensities typical of many endemic Sub Saharan African countries, and this was reflected in the level of genetic diversity observed in the populations. Furthermore, our findings are consistent across all three sites. Nevertheless, patterns of parasite mixing may differ between populations based on distinctive features such as geographic isolation and patterns of human movement. Further data are required to make more general conclusions. Furthermore, as transmission continues to decline and malaria programmes gradually shift their focus from control to elimination, the analysis of parasite gene flow between different transmission foci, e.g. Kilifi and Rachuonyo South, will become increasingly important in informing the mitigation measures needed to prevent importation of parasites as a result of human movement and migration. These analyses were not carried out in the current study since the numbers of common SNPs between the two Kenyan sites was low, and we only had parasites from one timepoint in Rachuonyo South, hence we were unable to conduct an informative analysis of gene flow between sites.
Finally, we used genetic data to show that there is high parasite movement and mixing within individual study sites. Additional analyses using gene flow models, e.g. as implemented in Migrate-N software, can be used to further validate our hypothesis of rapid gene flow and to confirm whether the parasites are part of a panmictic population or whether there exists underlying population structure, as well as to determine directionality of parasite movement between different populations, assuming that such populations can be identified within the region.
In conclusion, we have shown that Plasmodium falciparum parasite populations mix evenly within The Gambia, Kilifi and Rachuonyo South and there appear to be no detectable geographical barriers to parasite movement over short distances within these sites. That said, autocorrelations of genotype were detected at the micro-epidemiological level. We would conclude that control strategies that efficiently target hotspots will likely benefit the wider community outside the hotspots at the District/County level (we are however unable to comment on larger geographical scales), although this is likely to be affected by factors such as the underlying transmission level, heterogeneity of transmission, and patterns of human movement 23 . On the other hand, based on the high level of parasite mixing observed at each study site, we would predict that ineffective application of control interventions such as mass drug administration that result in residual foci of transmission would lead to rapid re-infection of the wider community, and also that parasites acquiring mutations conferring drug resistance or immunological escape would spread rapidly at the micro-epidemiological level. This underscores the need for effective and sustained control until malaria elimination is achieved.

Competing interests
No competing interests were disclosed. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supplementary figure 1. Time-distance interaction curves showing the effect of distance on changes in genetic variation between
Gambian P. falciparum parasite pairs over time. The analyses were carried out for a) pooled SNPs, b) EBA-175, c) AMA1 and d) 'other' SNPs (all SNPs excluding EBA175 and AMA1). Dashed lines represent time intervals separating parasite pairs at 1 day (red), 1 month (green), 6 months (blue) and 1 year (purple). The interaction term was log-transformed prior to running the analysis.
Click here to access the data. 1. This study analyzed large sample sets of malaria parasites taken from the western and eastern coasts of Africa (The Gambia and Kenya) and genotyped at 276 SNPs. For two of the sample sets, parasites were collected at different time points, allowing identification of population changes over time and space. Overall, the analysis was sound and results were well explained. The authors also notified the limitations of the study. For example, inclusion of additional parasite samples between these western and eastern sites, and use of more SNP markers would validate whether the conclusions drawn here represent the whole African continent.

Comments:
The assumption for comparing the temporally collected samples is that malaria case numbers have been reduced, which might lead to genetic isolation and structuring of parasite populations. It would be great if malaria epidemiology at the beginning and end of sample collection in the sites where samples were collected is clearly stated. It is possible that despite the overall reduction in malaria cases, some of the sites may represent hotspots where malaria epidemiology remained 1.
malaria cases, some of the sites may represent hotspots where malaria epidemiology remained more or less unchanged over the time. As a result, this would make the parasite populations and genetics relatively stable over the time.
The inclusion of numerous SNPs for this type of analysis is a nice practice. However, the authors may want to separate those that are clearly under selection (such as EBA175 and AMA1), since these mutations are subject to strong immune selection and will have different evolutionary trajectories as compared to more neutral SNPs.
More detailed comparison of the two Kenyan sites might be interesting to see whether gene flow between these sites exists, given that these sites are relatively closely located, yet separated by potential gene flow barriers (such as the rift valley).

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate? I cannot comment. A qualified statistician is required. We are grateful for this review and the helpful comments and suggestions that have been made. We have included a point-by-point response (in bold) to the issues raised.
The assumption for comparing the temporally collected samples is that malaria case numbers have been reduced, which might lead to genetic isolation and structuring of parasite populations. It would be great if malaria epidemiology at the beginning and end of sample collection in the sites where samples were collected is clearly stated. It is possible that despite the overall reduction in malaria cases, some of the sites may represent hotspots where malaria epidemiology remained more or less unchanged over the time. As a result, this would make the parasite populations and malaria cases, some of the sites may represent hotspots where malaria epidemiology remained more or less unchanged over the time. As a result, this would make the parasite populations and genetics relatively stable over the time.
The following statement has been added to show the changing epidemiology of malaria during the study period: "Over the study period, malaria transmission, as measured by malaria slide positivity rate, fell from 56% in 1998 to 7% in 2009 in Kilifi , and rose slightly in Fajara and Brikama in the Gambia ." There was no data showing temporal variation in malaria transmission in Rachuonyo South because the samples were collected in a single cross-sectional survey.
The inclusion of numerous SNPs for this type of analysis is a nice practice. However, the authors may want to separate those that are clearly under selection (such as EBA175 and AMA1), since these mutations are subject to strong immune selection and will have different evolutionary trajectories as compared to more neutral SNPs. More detailed comparison of the two Kenyan sites might be interesting to see whether gene flow between these sites exists, given that these sites are relatively closely located, yet separated by potential gene flow barriers (such as the rift valley).
The current study focused on parasite movement and mixing within small, geographically-defined areas, and hence concentrated on analysing parasite genetics within individual sites. Furthermore, we did not have sufficient number of SNPs typed in these two populations to carry out a meaningful comparison. That said, your comment is important, and has been noted as a recommendation for future work, in the statement below: "Furthermore, as transmission continues to decline and malaria programmes gradually shift their focus from control to elimination, the analysis of parasite gene flow between different transmission foci, e.g. Kilifi and Rachuonyo South will become increasingly important in informing the mitigation measures needed to prevent importation of parasites as a result of human movement and migration. These analyses were not carried out in the current study since the numbers of common SNPs between the two Kenyan 1 2 out in the current study since the numbers of common SNPs between the two Kenyan sites was low, and we only had parasites from one timepoint in Rachuonyo South, hence we were unable to conduct an informative analysis of gene flow between sites" No competing interests were disclosed. 1.

2.
"On the other hand, based on the high level of parasite mixing observed at each study site, we would predict that ineffective application of control interventions such as mass drug administration that result in residual foci of transmission would lead to rapid re-infection of the wider community, and also that parasites acquiring mutations conferring drug resistance or immunological escape would spread rapidly at the micro-epidemiological level. This underscores the need for effective and sustained control until malaria elimination is achieved." No competing interests were disclosed. Competing Interests: 11 April 2017 Referee Report doi:10.21956/wellcomeopenres.11628.r21345

Cristian Koepfli
University of California, Irvine, Irvine, CA, USA This is a relevant study, assessing the ability to identify small-scale foci of transmission and parasite gene flow to surrounding areas based in SNP-typing. While it is overall clearly presented and well written, more detail, in particular in the results section, would help to better understand the data, and to assess its power.

Specific comments:
Abstract: Please state how many samples and SNP markers were included in the final analysis.
I wonder whether "relatedness in space and time" is the correct term, or "distance in space and time" would be more appropriate.

Results:
In the part on "Identification of geographical barriers to parasite movement" it would be useful to include the range of prevalence or MPF per pixel analyzed. The second paragraph of this part is difficult to follow, as the term 'cluster' is used consistently, without further indication on what the clusters represent. It would help to include a sentence describing that spatial clusters were analyzed based on the PCA values of all isolates found within the cluster. Thus, clusters of isolates differing from all other isolates were identified. The same is the case in the discussion. What sizes were the clusters identified, and how many haplotypes were included per cluster? Figure 5: Given that for almost every pair of samples the number of days differs, how were the days for the different curves calculated? I assume each color represents a range, yet only a point estimate is given.
Also, please indicate in brackets for each curve the number of samples included. For example, how many samples were available for the 1-day and 31-days analysis in The Gambia? Could the apparent reduction in SNP difference at 10 km be a chance finding due to limited sample size?
Including the number of SNPs analyzed in each population would further help to interpret the data. 3. 1.

2.
Including the number of SNPs analyzed in each population would further help to interpret the data. E.g. it is interesting that in Rachuonyo South the proportion of different SNPs is approx. twice as high as in the other sites, yet this is only evident when Figure 5 is compared to Table 1.
Would it be possible to include confidence intervals for the 1-day curves in the figure? This would help to understand the power of the data. For example, the statement in the abstract "Genetic relatedness of sample pairs is predicted by relatedness in space and time" suggests that genetic relatedness can be inferred, once the distance by space and time is known. This is however difficult to assess without more detail on the variance of the data.
In We would like to sincerely thank the reviewer for taking time to review this work. We appreciate the comments and issues raised, and have addressed them on a point-by-point basis as indicated (in bold).

Abstract:
Please state how many samples and SNP markers were included in the final analysis.
The following statement has been inserted in the abstract to show the number of samples and SNPs used in the final analysis. "Using 107, 177 and 82 SNPs that were successfully genotyped in 133, 1602, and 1034 parasite isolates from The Gambia, Kilifi and Rachuonyo South district, respectively, we show that there are no discrete geographically restricted parasite sub-populations." I wonder whether "relatedness in space and time" is the correct term, or "distance in space and I wonder whether "relatedness in space and time" is the correct term, or "distance in space and time" would be more appropriate.
"Relatedness in space and time" has been replaced with "distance in space and time".

Results:
In the part on "Identification of geographical barriers to parasite movement" it would be useful to include the range of prevalence or MPF per pixel analyzed. The second paragraph of this part is difficult to follow, as the term 'cluster' is used consistently, without further indication on what the clusters represent. It would help to include a sentence describing that spatial clusters were analyzed based on the PCA values of all isolates found within the cluster. Thus, clusters of isolates differing from all other isolates were identified. The same is the case in the discussion. What sizes were the clusters identified, and how many haplotypes were included per cluster?
The following statement has been added to show the MPF and parasite prevalence range analysed per pixel. The interquartile range of both measures at each pixel size has also been added.
"The range of MPF and parasite prevalence per pixel varied depending on the size of the pixels analysed. In the Kilifi population, MPF ranged from 0 -100% (interquartile range (IQR) = 20%) for the 0.5km pixels; 0 -100% (IQR = 14%) for the 1.0km pixels; 20 -83% (IQR = 7%) for the 2km pixels; and 33 -63% (IQR = 4.7%) for the 4km pixels. In the Rachuonyo South population, PCR positive prevalence varied from 0 -75% (IQR = 19.4%) for the 0.5km pixels; 0 -47% (IQR = 17.4%) for the 1.0km pixels; 3.5 -35.8% (IQR = 14.3%) for the 2km pixels; and 6.2 -33.4% (IQR = 7.8%) for the 4km pixels." The second paragraph under 'spatial autocorrelation' heading in the methods section has been expanded to indicate what clusters represent and how they are identified, and now reads as follows: "Spatial scan statistics were calculated using SaTScan software and were run separately for each study site. The statistics involved running a purely spatial, retrospective analysis based on a normal probability distribution model using continuous variables (PC scores) and looking for areas with clusters of high PC scores. Latitude and longitude coordinates were used to represent the geographical locations of specific parasites, whereas principal component scores were used to represent individual parasite genotypes. During the analysis, a scanning window that gradually varies in size from including only a single homestead up to 50% of the population moves over the geographical space and at each window size and location, the ratio of parasites with high PCs inside the window versus outside the window is calculated. The window with the highest ratio is noted down as a cluster and its statistical significance is determined after accounting for multiple comparisons using random permutations." The number of samples contained in each significant cluster has been added in the results section and now reads as follows: "In Kilifi, one cluster with a radius of 1.54 km (p=0.01) containing 15 parasites was "In Kilifi, one cluster with a radius of 1.54 km (p=0.01) containing 15 parasites was detected, while in Rachuonyo South, a smaller cluster of genetically distinct parasites was detected with a radius of 0.5 km (p=0.001) containing 14 parasites".
And the following section has been added in the discussion section to make it clear that the clusters were based on PC scores: "However, using scan statistics we identified only two specific clusters of distinct parasite sub-populations based on PC scores, one in Kilifi and another in Rachuonyo South." Figure 5: Given that for almost every pair of samples the number of days differs, how were the days for the different curves calculated? I assume each color represents a range, yet only a point estimate is given.
The following statement has been added in the results section to show how the graphs were generated: "Since the number of days differed for almost all parasite pairs, dummy data were included in the regression analysis to enable the generation of time-distance interaction graphs. For each study site, a distance range of 1 -10km (with an interval of 0.1km between adjacent distances) was used. Temporal distance with 14 and 10 day intervals were assigned to parasite pairs in Kilifi and the Gambia, respectively, whereas time was not considered for the Rachuonyo South population. Constant SNP differences of 14, 10 and 8 were used for parasite pairs in Kilifi, Rachuonyo South and the Gambia, respectively." Also, please indicate in brackets for each curve the number of samples included. For example, how many samples were available for the 1-day and 31-days analysis in The Gambia? Could the apparent reduction in SNP difference at 10 km be a chance finding due to limited sample size?
The number of samples used to draw each curve has been noted in the Figure legend. We agree with the reviewer that the decrease past 10km is likely due to limited sample size and is there probably not significant.
Including the number of SNPs analyzed in each population would further help to interpret the data. E.g. it is interesting that in Rachuonyo South the proportion of different SNPs is approx. twice as high as in the other sites, yet this is only evident when Figure 5 is compared to Table 1.

The number of SNPs analysed in each population has been added to the figure legend.
Would it be possible to include confidence intervals for the 1-day curves in the figure? This would help to understand the power of the data. For example, the statement in the abstract "Genetic relatedness of sample pairs is predicted by relatedness in space and time" suggests that genetic relatedness can be inferred, once the distance by space and time is known. This is however difficult to assess without more detail on the variance of the data. difficult to assess without more detail on the variance of the data.

Figure 5 has been regenerated to show the confidence intervals for the 1-day curves.
In Table 2, what is the unit of the results showed? I assume it is SNP-difference/day (or SNP-difference/km), with days and distance log-transformed. Please state if/how data was transformed.
The following statement has been added to Table 2 to show the units of measurement of the results.
"Values represent the change in the number of SNP differences between parasite pairs per day (time), per kilometre (distance) and per day/kilometre (time-distance interaction). Time, distance and the product of time and distance (time-distance interaction) were log transformed prior to running the regression analyses."

Discussion:
Paragraph 3 of the discussion could be expanded. At what spatial scales was population structure found in previous studies (as compared to the approx. 50 km range of the present study)? Have any of these studies included relatedness? This information would help to assess the feasibility to identify foci of higher transmission, and to estimate the level of gene flow to surrounding areas in different transmission settings.
The paragraph has been expanded and now reads as follows: "Lack of genetic structuring of parasite populations observed in this study is indicative of a population that is well mixed. This observation of a highly mixing parasite population is in agreement with results of similar studies of parasite genetic diversity and population differentiation using microsatellites , immune selected genes and SNPs . These studies were carried out in parasite populations from different geographical regions representing a diverse range of transmission intensities from highest in Africa and oceania, intermediate in Southeast Asia, and lowest in south America. However, other studies have shown population structure when looking at the same population , although these analyses were carried out on larger geographical scales than those analysed here and mostly involved analyses at provincial, country or continental levels. Population structure was most evident in regions with low transmissiion intensities such as south America or southeast Asia, and less evident in Africa where transmission intensity is much higher." The number of SNP differences plateaus at approx. 1 km in The Gambia, 3 km in Kilifi, and increases up to 10 km in Rachuonyo South. Are there possible explanations for these differences due to the characteristics of the local parasite populations?
Although the exact reason for the plateau observed is currently unknown, the following statement has been added to postulate possible reasons for the observation: We would like to sincerely thank the reviewer for taking time to review our work. We appreciate your comments, and have included a point-by-point response to them as follows (our responses are in bold).
How can heterozygous genotypes be detected in haploid populations of the parasite?
The use of the words "homozygous" and "heterozygous" in the context of haploid organisms has been clarified to mean single parasite genotype infections and mixed parasite genotype infections, respectively.
As noted by the authors themselves, using 276 SNPs is rather limited. Genetic studies dealing with human populations at nowadays routinely rely on 500000 SNPs or more. One main feature of such studies is that microgeographical structures are deteted mostly from low frequency variants and rare variants, which of course are undetectable when using a limited set of markers. Moreover, these low frequency and rare variants are supposed to be highly relevant for phenotypic expression, in particular disease susceptibility and are largely responsible for recent and localized evolution in human populations. (see for example Leslie et al. (2015) ). It is most probable that these patterns exist in parasite populations too. The authors should discuss this point more, since it is probably one of the main avenues of future researches in microbiology.
The following additional information has been added in the discussion section to address this point: "Advances in sequencing technologies have increased the use of whole genome sequence data in the analysis of parasite population genetics, and this has P. falciparum led to the identification of hundreds of thousands of SNPs, most of which are present at very low frequencies especially in African parasite populations . Additional analyses will require the use of whole genome sequence data to identify rare variants and distinguish between closely related parasites, thus allowing parasite population structure to be analysed at fine spatial scales." No competing interests were disclosed. Competing Interests: 1