Temporal and spatial dynamics of Plasmodium falciparum clonal lineages in Guyana

Plasmodium parasites, the causal agents of malaria, are eukaryotic organisms that obligately undergo sexual recombination within mosquitoes. In low transmission settings, parasites recombine with themselves, and the clonal lineage is propagated rather than broken up by outcrossing. We investigated whether stochastic/neutral factors drive the persistence and abundance of Plasmodium falciparum clonal lineages in Guyana, a country with relatively low malaria transmission, but the only setting in the Americas in which an important artemisinin resistance mutation (pfk13 C580Y) has been observed. We performed whole genome sequencing on 1,727 Plasmodium falciparum samples collected from infected patients across a five-year period (2016–2021). We characterized the relatedness between each pair of monoclonal infections (n = 1,409) through estimation of identity-by-descent (IBD) and also typed each sample for known or candidate drug resistance mutations. A total of 160 multi-isolate clones (mean IBD ≥ 0.90) were circulating in Guyana during the study period, comprising 13 highly related clusters (mean IBD ≥ 0.40). In the five-year study period, we observed a decrease in frequency of a mutation associated with artemisinin partner drug (piperaquine) resistance (pfcrt C350R) and limited co-occurence of pfcrt C350R with duplications of plasmepsin 2/3, an epistatic interaction associated with piperaquine resistance. We additionally observed 61 nonsynonymous substitutions that increased markedly in frequency over the study period as well as a novel pfk13 mutation (G718S). However, P. falciparum clonal dynamics in Guyana appear to be largely driven by stochastic factors, in contrast to other geographic regions, given that clones carrying drug resistance polymorphisms do not demonstrate enhanced persistence or higher abundance than clones carrying polymorphisms of comparable frequency that are unrelated to resistance. The use of multiple artemisinin combination therapies in Guyana may have contributed to the disappearance of the pfk13 C580Y mutation.


Introduction
Genomic data from pathogens, vectors, and/or human hosts can complement traditional epidemiological data on disease incidence and prevalence to inform decisions regarding control.
In the case of malaria, several distinct use cases for genomic epidemiology have been previously identified [1], including the identification of imported cases and transmission hotspots [2,3], as well as informing strategies for local disease elimination by documenting connectivity among parasite populations mediated by human movement [4].Most importantly, genomic data from malaria parasites can play an important role in surveillance of emerging drug resistance markers [5].Resistance has arisen to every widely deployed antimalarial [6], and molecular surveillance has been endorsed by the WHO as a core intervention for maintaining the efficacy of current malaria drug treatment regimens [7].Genetic surveillance of drug or insecticide resistance is typically conducted using genotyping data from specific polymorphisms associated with resistance [4,8].However, whole genome sequencing (WGS) data and genome-wide genotyping assays can inform understanding of the context for the origin and spread of mutations, especially in cases where compensatory or epistatic mutations are required to generate a high-fitness resistance genotype capable of spreading quickly [9].While measurable phenotypic resistance may be conferred by individual mutations, other genomic changes are often required for those mutations to be evolutionarily successful, with examples in Plasmodium malaria parasites [10][11][12], bacteria [13] and other pathogens [14].
Resistance has been arising in a small number of specific geographic locations to artemisinin, which is administered with one or more partner drugs as artemisinin combination therapy (ACT) as the first line treatment for malaria caused by Plasmodium falciparum in most of the world.
Delayed parasite clearance following ACT treatment was first observed in the Greater Mekong Subregion (GMS) of Southern Asia in early 2000s [15,16].More recently, mutations associated with reduced susceptibility to artemisinin have also been detected in East Africa [17][18][19][20][21] and Papua New Guinea [22].The most important artemisinin resistance mutation, a C to Y substitution at codon 580 (C580Y) in the propeller domain of a kelch-domain-containing protein on chromosome 13 (pfk13) was first observed in the Americas in samples collected in Guyana in 2010, where five out of 94 symptomatic cases were found to carry the pfk13 C580Y mutation [23].
In 2014, a therapeutic efficacy study (TES) from Guyana failed to detect clinical artemisinin resistance [24], but sample size was likely too low to recruit subjects with low-frequency resistance mutations.The pfk13 C580Y mutation was observed in 14 out of 854 clinical samples in a resistance surveillance study conducted in Guyana from 2016-2017, and through whole genome sequencing we determined that all of these samples represented a single clonal parasite lineage, despite being observed in disparate regions of the country [25].This observation of a single clonal background for the pfk13 C580Y mutation in Guyana was unexpected because P. falciparum is a eukaryotic parasite that undergoes sexual recombination in mosquitoes as an obligatory component of its life cycle.However, when a mosquito bites a human host with a monoclonal infection (caused by a single parasite genomic lineage), parasites do not have an opportunity to undergo sexual outcrossing in the mosquito, and instead perform selfing, resulting in the perpetuation of the genomic lineage present in the previous human host.Malaria transmission levels are low in Guyana relative to many settings in sub-Saharan Africa, and therefore most infections are monoclonal, resulting in frequent clonal transmission.Therefore, a null hypothesis to explain the observation of pfk13 C580Y on a single clonal background could simply invoke low transmission in Guyana as a causal mechanism.However, a plausible alternative hypothesis is that the pfk13 C580Y mutation was observed on a single clonal background because that genomic lineage contained important compensatory or epistatic mutations, related to the phenotype of artemisinin resistance directly or resistance to one or more partner drugs commonly administered in ACTs.Historically, resistance to antimalarials has originated de novo in low-transmission settings like Southeast Asia or the Americas and has only later spread to sub-Saharan Africa where malaria is much more common [26], leading to the hypothesis that low sexual outcrossing rates in such settings could facilitate the emergence of highfitness resistance genotypes by preserving key combinations of alleles (in addition to factors such as lower immunity and higher drug pressure).Clonality has been associated with the emergence of pfk13 C580Y P. falciparum in Cambodia and its subsequent spread throughout the GMS [16,26], perhaps facilitated by an additional mutation in this lineage (plasmepsin 2 and/or plasmepsin 3 gene amplification) that confers resistance to an important artemisinin partner drug.
In East Africa, studies from Uganda [19] and Eritrea [27] reported evidence of emergence of resistance through clonal propagation with an increase in prevalence of pfk13 mutations.The official first-line for malaria in Guyana is the ACT artemether-lumefantrine (AL), and no lumefantrine (LMF) resistance mutations are known to be segregating in Guyana P. falciparum populations.However, an important context for malaria transmission in Guyana is among goldminers working in forested regions who are known to frequently self-medicate with the ACT dihydroartemisinin (DHA) -piperaquine (PPQ) -trimethoprim (TMP; DHA + PPQ + TMP; Artecom) tablets [28,29].At least two mutations that are segregating in Guyana confer resistance to piperaquine: a point C350R mutation in the chloroquine resistance transporter (pfcrt) gene that is endemic to the Guiana Shield region and has been increasing in frequency over the last 20 years [28,30,31], and copy number amplification of the plasmepsin 2 (Pfpm2 -PF3D7_1408000) and/or plasmepsin 3 (Pfpm3 -PF3D7_1408100) genes [30].The pfcrt C350R mutation and plasmepsin 2/3 amplifications interact epistatically to yield piperaquine resistance [32], adding credibility to the hypothesis that clonal transmission may be adaptive under DHA-piperaquine pressure.
In the present study we generated whole genome sequencing data from P. falciparum clinical samples collected in Guyana between 2016-2021 to profile the temporal and spatial dynamics of clonal parasite lineages.We identify circulating clonal components (referred to as clones), defined as groups of genomically indistinguishable parasites identified under a graph-based framework [33], and we explore whether limited sexual outcrossing may have been conducive to the de novo origin of the pfk13 C580Y mutation in Guyana.We specifically explore the representation of pfk13 C580Y, pfcrt C350R, and Pfpm2/3 gene amplifications in clonal and unique parasite genomic backgrounds, and their co-occurence in frequency vs. rare clonal lineages.Further, we profile new signatures of selection in the local parasite population using this deep population genomic dataset to determine whether previously uncharacterized mutations may also drive clonal dynamics, or whether the persistence and prevalence of clonal lineages in Guyana are driven by stochastic factors.

Temporal and Spatial Clonal Dynamics in Guyana
We performed selective whole-genome amplification (sWGA) on 1,727 samples collected from Guyana between 2016 and 2021 across three time periods (Fig. 1).A total of 264 genomes (15.3%) did not meet the quality criteria of at least 30% of the genome covered at ≥ 5-fold coverage, resulting in 1,463 samples suitable for analysis.Of this set, 54 samples were classified as multiclonal infections (F ws < 0.7) and were excluded from subsequent analyses.The final dataset for relatedness analysis contained 1,409 monoclonal genomes (Table S9) with an average pairwise IBD across the entire dataset of 0.283 (SD = 1.114 -Fig.S1).The final dataset obtained was composed of 736 genomes from 2016/2017; 130 genomes from 2018/2019 which were collected as part of a therapeutic efficacy study (TES); and finally 523 genomes from patient samples collected in 2019/2021.To explore patterns of relatedness due to shared recent common ancestry, pairwise identity-by-descent (IBD) values were computed between haploid genotypes (Fig. 2).Genome-wide mean IBD estimates across samples revealed patterns of shared ancestry.
Network analysis identified 160 clones (C), which were defined as groups of at least two samples with a mean pairwise IBD ≥ 90%, and 332 singletons.Some of these clones formed larger highly related clusters, defined as a group of multiple clones (n ≥ 3) which displayed a mean IBD ≥ 40% and grouped together in the hierarchically-clustered dendrogram as portrayed in Fig.  4).In terms of clonal persistence, this haplotype was among the top 20% of multi-occurrence clones (n=130).We investigated nonsynonymous (NSY) mutations with a similar allelic frequency (MAF = 0.007 ± 0.05) as pfk13 C580Y, screening 2,360 NSY mutations for their relative clonal size and clonal persistence (Fig. 5).The temporal persistence of the pfk13 C580Y mutation above the mean clonal duration (t = 287.4days, p < 0.211, t C580Y =418.0 days) and clonal size was above average (n C580Y =8.0 p < 0.211, n mean =6.0) but below the 95th percentile of polymorphisms in the same frequency class in each case.Two occurrences of a previously undescribed pfk13 mutation (G718S) were also observed.The two samples were collected the same week in November 2020 in Aranka River in Region 7.They also carried the pfk13 K189T mutation.These samples belonged to a clonal background (C#321) composed of four samples which was first detected in April 2018, but the other members of this clonal background did not have sufficient coverage at this position to permit allele identification.
Decrease in pfcrt C350R frequency across the five year study period.
In the pfcrt gene (PF3D7_0709000), which encodes a transmembrane digestive vacuole protein known to modulate resistance to chloroquine and other drugs [10], Allelic positions 72, 76, 220, 326 (wildtype) and 356 in pfcrt were fixed.The frequency of pfcrt C350R in the dataset was 54.04% (n=709) and was found in 222 clones (Fig. 2).Additionally, six samples harbored a previously undocumented coding polymorphism in pfcrt: D329N.The earliest observation of the D329N mutation was obtained in September 2018 and was sampled in Georgetown as part of the TES.The D329N mutation was found in three clones (C#173, C#9, C#402), which were each composed of two samples.These samples exhibited the pfcrt C350 wildtype allele and were found in different highly related clusters.Between the two study periods, a change in pfcrt C350R frequency was observed.In 2016/17, pfcrt C350R was present in 73.3% of samples (n=478) while in 2020/2021, the frequency of the mutation was 36.2% (n=191).This frequency reversion to the wildtype allele could also be observed within a highly related cluster (Table S7).In 2016/17, highly related cluster 6 displayed one predominant clone (C#134) carrying the pfcrt C350R mutation.The clone was found primarily in Mid Essequibo but also appeared in six other spatial clusters (Fig. 3 & S4, Table S8).In 2020/21, samples in this highly related cluster return to the wildtype (in C#137, C#135, C#136 and C#138).These samples were still widely distributed, with C#137 occuring in nine spatial clusters (Fig. 6).Evidence of multiple events of pfcrt C350R mutation was observed.
The wildtype and pfcrt C350R were both observed in nine clones.For instance, C#45 contained four samples with pfcrt C350R and seven representing the wildtype (Fig. 6).The clone C#268 harboring the pfk13 C580Y did not carry pfcrt C350R.When investigating whether pfcrt C350R had an impact on clonal persistence or clonal size, no significant difference was observed.The average duration of clones carrying the mutation was 268.0 days (p < 0.810), while the average duration of clones representing other nonsynonymous (NSY) mutations of comparable allele frequency (MAF = 0.46 ± 0.05) was 291.0 days (Fig. 5C, n = 683 NSY mutations).The size of clones harboring pfcrt C350R (mean = 6.9, p < 0.526) was also similar to the size of clones representing these comparator mutations (mean size of 6.9 samples) (Fig. 5D).Mutations significantly associated with prolonged clonal duration included one mutation in falcilysin gene (PF3D7_1360800, n= 134.5 days, p < 0.001) as well as two mutations in PF3D7_1133400 (AMA1 -apical membrane antigen 1, n=156.4/150.0days, p < 0.001).

Co-occurrence of plasmepsin 2/3 duplication and pfcrt C350R
The combination of pfcrt C350R and plasmepsin 2/3 copy number amplification has been recently demonstrated to confer piperaquine resistance in the Guiana Shield [30].Pfpm2/3 copy number status was recovered from 62 samples in 2016/17 (8.0%) [30].Although the information was only available for a limited number of samples, a Χ 2 test revealed a significant association between pfpm2/3 copy number and C350R (P < 0.035).We assessed pfpm2/3 copy number for 401 samples in the 2020/21 dataset.A total of 96 (15.2%) samples possessed pfcrt C350R in combination with an increase in pfpm2/3 copy number while 87 samples (21.7%) exhibited pfcrt C350R with a single copy of pfpm2/3.Ninety-six out of 245 (39.2%) wildtype samples presented multiple copies of pfpm2/3.No significant association between pfpm2/3 copy number and pfcrt C350R was observed during the 2020/21 period (Χ 2 = 0.89, P < 0.64).Highly related clusters appeared to carry the pfcrt C350R mutation heterogeneously.Highly related clusters with more than two samples displayed a frequency of pfpm2/3 copy number variation of 37.8% (Fig. 7).Only four out of the 60 clones investigated carried the duplication homogeneously highlighting the genomic lability of this duplication.
While the frequency of pfcrt C350R decreased between 2016 and 2021, two highly related clusters exhibited a frequency increase and were predominantly carrying the mutation (Table S7).

Mutations in drug resistance genes
The rise of drug resistance polymorphisms has the potential to drive clonal dynamics.In mdr1
The selection landscape of these two periods 2016/2017 and 2020/2021 was also investigated using isoRelate [34] to detect genomic regions exhibiting enhanced relatedness, with a false discovery rate of 0.01.The analysis was run on clones sampled across more than three months in

Discussion
In this study, we profiled the dynamics of P. falciparum over a five year study period using the deepest whole genome sequencing dataset yet produced for this parasite species from a single country.In contrast to the GMS and East Africa, where clonal transmission and enhanced population relatedness were directly related to the emergence of mutations conferring resistance to ACTs [27,35], Guyana offers a different perspective on clonal dynamics in the context of drug resistance emergence in a low-transmission setting.Stochastic processes with intermittent recombination appear to be the dominant mechanism driving clonal diversity rather than a selective advantage obtained from particular polymorphisms favoring a specific clonal background.

305Impact of artemisinin on clonal dynamics in Guyana
Resistant lineages can circulate at low frequencies for years before becoming dominant.In this study, a total of 160 clones aggregated into 13 highly related clusters were observed.Two highly related clusters present at the beginning of the study disappeared by 2020, while four highly related clusters increased in frequency (Table 1).Malaria transmission in the Guyana shield is largely driven by mobile populations working in gold mining or other forest-associated professions [36].
Evidence of clonal dispersal among spatial clusters was best represented by highly related cluster 6, which has two clones spreading to seven and nine spatial clusters in a limited amount of time (Fig. 4).In 2016/2017, one clonal component (C#134) dominated and was preferentially found in Mid Essequibo.By 2021, the former clonal component had disappeared and a related clone (C#137) appeared to be circulating predominantly in Potaro.A cautionary note regarding Pf sampling is needed as this dataset was assembled through different sampling schemes performed at different health centers.Parasite origin inferred from patient travel recollection may not be consistently precise and regions with high mining activities might be over represented (Fig. S4).
In 2004, Guyana was the first country on the continent to implement artemether-lumefantrine (COARTEM ® ).Therefore, constant artemisinin pressure and shifting exposure to lumefantrine, piperaquine, and perhaps other partner drugs has imposed heterogeneous selective pressure on P. falciparum lineages.However, only limited evidence of allelic change responding to artemisinin drug pressure has been observed.A pfkic6 (Kelch13 Interacting Candidate; PF3D7_0609700) Q1680K polymorphism increased by 31.65% (Table S3) in the five year study period.The NSY mutation was present in clones which persisted longer than average (Δ=31.8,p < 0.001) and pfkic6 is a gene which could potentially play a role in artemisinin (ART) resistance given its association with the resistance-associated PfK13 protein [37,38].Other polymorphisms that appeared to be favored in the Guyana landscape were associated with potential resistance to artemisinin (Table 2) and their prevalence should be closely monitored, but polymorphisms driving resistance in other regions did not show any obvious signs of selection.

Selection by artemisinin partner drugs
The emergence of drug resistance in the Guyana shield is of concern, considering that resistance to chloroquine and sulfadoxine-pyrimethamine emerged almost simultaneously and in an independent manner in both South America and Southeast Asia [11,39].In Guyana, 54% of gold miners self-medicate to treat fever using Artecom (DHA+PPQ+TMP) tablets before seeking care [40].The association of the pfcrt C350R allele with an amplification of plasmepsin (xpfpm2/3) has been shown to strengthen resistance phenotypes to piperaquine [28,30].In the current study, we observed a reduction in frequency of pfcrt C350R from 73.0% to 24.2% across five years indicating a potential reduction in piperaquine pressure.We can speculate that a change in dominant ACT therapy from DHA+PPQ+TMP to artemether-lumefantrine could have occurred.
Erratic use of DHA+PPQ+TMP during a period of high prevalence of the pfcrt C350R mutation could have contributed to the emergence of the pfk13 C580Y mutation.Subsequent increase in the use of artemether-lumefantrine may have reduced the pressure to maintain pfcrt C350R, and eliminated the pfk13 C580Y mutation primarily through success of the partner drug.Further potential evidence for reduced DHA+PPQ+TMP self-medication in recent years is the reduced prevalence of parasite genomes containing both pfcrt C350R and plasmepsin duplication.
In Southeast Asia, an increase in copy number of the plasmepsin 2 (pfpm2) and/or plasmepsin 3 (pfpm3) genes is associated with piperaquine resistance [35,41].These copy number amplifications have been observed to enhance piperaquine resistance in vitro through epistatic interaction with the pfcrt C350R mutation [30].As observed in French Guiana [30], we found multiple mutational events for pfcrt C350R occurring within a short timespan.Plasmepsin duplication was also highly genomically labile, varying within and among conserved clonal lineages (Fig. 7).The gene amplification appeared more variable compared to the emergence of polymorphism.The frequency of these phenomena unique to this part of the world make it difficult for a clone to thrive to the extent observed in the GMS.Gene copy number may appear as a strategy for regulating expression under environmental stresses [42].In this context, the plasticity of pfpm2/3 might reflect a more rapid adjustment of the parasite responding to heterogeneous drug exposure.For instance, in highly related cluster 4, the dominant highly related cluster circulating in Lower Mazaruni in 2020/21, 14 of45 samples displayed both pfcrt C350R and xpfpm2/3, which might reflect localized recent selection by piperaquine.

Other candidate variants associated with clonal dynamics
Although this study primarily attributes recent spatiotemporal dynamics of parasite clones in Guyana to stochastic processes (e.g., sporadic outcrossing, periods of low-transmission bottleneck) rather than to selection towards the preservation of specific multi-locus haplotypes, we do not suggest that meaningful selective processes are entirely absent.For instance, the clonal background containing pfk13 C580Y was observed in six spatial clusters across 418 days, where the average clone was found in 2.36 spatial clusters and lasted on average 287 days (Fig. 4).It is therefore possible that the pfk13 C580Y mutation improved clone fitness for a period of time.We also noted the previously unobserved pfk13 G718S mutation in C#321, further sign of autochthonous pfk13 polymorphism in Guyana (Fig. 3).Furthermore, we observed persistent and large clones carrying two NSY mutations in AMA1 as well as a NSY mutation (MAF = 0.46 ± 0.05) in falcilysin (PF3D7_1360800) (Table S3-S4).The latter additionally featured among the 61 NSY mutations which increased in frequency between 2016 and 2021.Falcilysin is a metalloprotease believed to be involved in hemoglobin digestion, and has been found to be a target of chloroquine, which inhibits its proteolytic activity [43].Given that degraded products of hemoglobin activate ART [44], it is possible that this polymorphism interferes with parasite clearance.
The outcrossing rate in Guyana appears to maintain sufficient haplotypic diversity in the population to prevent the long-term dominance of specific clones.However, four clones were sampled over four years, indicating the possibility of longer-term clonal persistence in the region.
The selection signal observed at pfcrt was conserved throughout the dataset as previously described in global P. falciparum populations [34] (Fig. 5).These results suggest that selection may yet be influencing clonal dynamics in Guyana, even if the impact of selection is not as stark as in the GMS [9].
Other NSY mutations which increased in frequency tended to be associated with gametocyte maturation, a process which is key to withstanding artemisinin pressure [45] because artemisinin clears only asexual parasites.Moreover, gametocyte production ultimately determines fitness because they are required for transmission.Three polymorphisms were found in transcription factor pfap2-g5 (PF3D7_1139300).Apicomplexan-specific ApiAP2 gene family is a well-known regulator of sexual commitment and gametocyte development [46][47][48].The gene appears as an important mechanism during the maturation of sexual stages through gene repression combined with other chromatin-related proteins [49].Transcription factors (AP2 genes) involved in the gametocyte development have been previously found to display the strongest signatures of selection in French Guiana [50].Seven other genes which increased in frequency are also related to gametocyte development.For instance, PF3D7_0904200 (PH domain-containing protein) transcripts have been shown to be enriched in gametocytes [51] and PF3D7_1474200 was found to be highly expressed in late-stage gametocytes [52].

Relevance of pfK13 C580Y mutation disappearance
Guyana represents the first country where the pfk13 C580Y mutation (or similar ART resistance mutations) have appeared and then subsequently disappeared rather than increase in frequency.
The mutation was restricted to a single clonal background and was last observed in April 2017.
This clonal background lacked the pfcrt C350R mutation, making it likely susceptible to PPQ, which has been subject to fluctuating use through self-medication in the country and might have led to this disappearance in the presence of efficacious artemether-lumefantrine treatment.
Previous therapeutic efficacy studies in the region have hinted at resistance to artemetherlumefantrine [53] and artesunate monotherapy [54] but evidence from TES in Guyana is lacking.
A modeling study exploring factors associated with the spread of pfk13 mutations found that deploying multiple first-line therapies was the best approach to postponing treatment failure [55].
The simultaneous use and potentially shifting balance of at least two ACTs in Guyana might have therefore led to the elimination of the pfk13 C580Y mutation and its clonal background.Clonal turnover in Guyana appears to be different from the patterns observed in other regions like South-East Asia and East Africa.In the GMS, artemisinin was initially used as monotherapy facilitating rapid resistance expansion via hard selective sweep [56].These observations indicate that drug resistance emergence does not result in the same patterns of clonal dynamics in different geographic locations, perhaps due to unique differences in disease epidemiology and drug pressure across settings.Further molecular surveillance of clonal dynamics is warranted in settings where it occurs, given the potential association of clonal transmission with both known and novel mutations associated with drug resistance.

Sample collection and spatial cluster mapping
We evaluated 1,727 clinical samples collected from malaria-diagnosed individuals between 2016 and 2021 who provided informed consent for genetic analysis of their parasite samples.
Samples were collected as dried blood spots on Whatman FTA cards.Samples dating from 2016-2017 (n=837) were collected for a resistance surveillance project [25].Samples dating from 2018-2019 (n=174) were collected in the context of a therapeutic efficacy study.Samples dating from 2020-2021 (n=716) were collected for a separate malaria molecular surveillance study from individuals diagnosed with P. falciparum infection (Fig. 1).Participants provided informed consent in accordance with the ethical regulations of the countries.
To define spatial clusters, we first matched travel history responses to a catalog of malaria survey sites used by the Guyana Ministry of Health (MoH).We then mapped survey sites onto a custom shape file summarizing the country's primary river and road coordinates and onto a raster map of motorized transport resistance [57] available at https://malariaatlas.org/.Sites were clustered based on river/road connectivity in the R package 'riverdist' [58], travel conductance using the R package 'gdistance' [59], and manual assessment of coordinates on river/road and resistance layers in QGIS.Samples were collected in specific recruitment locations and patient travel history was documented.To investigate spatial patterns, 20 spatial clusters were defined following roads and rivers access (Fig. 1a).Patient travel history revealed that a majority of infections were acquired in Lower Mazaruni River in Region 7 (n=434, 36.1%),followed by Potaro River in Region 8 (n=162, 13.5%), as well as along the Cuyuni River (Table S1

Genomic data generation
DNA extraction was performed using two approaches according to year of collection.Samples from 2016-2017-2018-2019 (n=1,011) were extracted from dried blood spots using the QIAamp DNA mini kit according to the manufacturer's instructions (Qiagen, Hilden, Germany).For samples from 2020-2021 (n=716), we performed DNA extraction on all patient samples using a ThermoFisher blood and tissue kit and a ThermoFisher Kingfisher instrument.We performed selective whole genome amplification (sWGA) [60] on all samples to enrich the proportion of parasite DNA relative to host DNA.We performed library construction using a NEBNext kit on the enriched DNA samples and sequenced them on an Illumina NovaSeq instrument using 150 bp paired-end reads.We aligned reads to the P. falciparum 3D7 v.3 reference genome assembly and called variants following the Pf3K consortium best practices (https://www.malariagen.net/projects/pf3k).We used BWA-MEM [61] to align raw reads and remove duplicate reads with Picard tools [62].We called SNPs using GATK v3.5 HaplotypeCaller [63].We performed base quality score and variant quality score recalibration using a set of Mendelian-validated SNPs, and restricted downstream population genomic analyses to SNPs observed in 'accessible' genomic regions determined to be amenable to high quality read alignment and variant calling [64].Individual calls supported by fewer than five reads were removed and any variant within 5 nucleotides of a GATK-identified indel was also excluded.
Samples exhibiting quality monoclonal genome data (>= 5x coverage for >30% of the genome) were included in relatedness analyses.The final dataset to investigate mutation comprised 74,357 SNPs.

Relatedness analysis using identity by descent
We performed analyses of relatedness by estimating pairwise identity by descent (IBD) between all monoclonal patient samples (n=1,409).We estimated IBD using the hmmIBD triplicates for the DNA samples.If the value was between 0.6 and 1.5, the copy number is estimated as 1, whereas if the value was between 1.5 and 2.4, the copy number estimated was 2.
We use the term xpfpm2/3 to designate the amplification of pfpm2 or pfpm3 and 1pfpm2/3 to denote one copy of both genes similarly to [30].
2016/2017 and in 2020/2021 (n 2016/2017 = 55 clones, n 2020/2021 = 49 clones).Selection signals (relatedness peaks) were consistent across different analysis runs with different representative samples of each clonal component (Fig.S7).This allowed investigation of signals of positive selection within "successful" clones to understand whether genes present on genomic segments within these clones were particularly important in Guyana.In 2016/2017, seven segments (159 genes) contained within four strong selection signals on chromosomes 2, 4, 7 and 9 were identified among long-lasting clones (see Fig.6and see supplementary results for details).On chromosome 9 (chr9: 61,342-208,725 and 318,311-432,047), the selection signal found in long-lasting clones in 2016/2017 was observed in both short-and long-lasting clones in 2020/2021.Nine genes (eleven mutations) presented a large increase in mutation frequencies ( FREQUENCY ≥ 28.4%, Table3) in the past five years: PF3D7_0902400 and PF3D7_0902500 (two serine/threonine protein kinase part of the FIKK family gene), PF3D7_0903300 (unknown function), PF3D7_0904200 (PH domain-containing protein), PF3D7_0905500 (unknown function).

Figure 2 -
Figure 2 -The mean IBD between samples highlighting highly related clusters (IBD ≥ 0.4).Different sampling years are indicated as well as the presence/absence of pfcrt C350R.

Figure 3 -
Figure 3 -Clonal temporal dynamics between 2016 and 2021.Clone (clonal component) C#1 in highly related cluster 1 was the largest clone present in the dataset (n=73 samples).C#268 is the clonal background harboring pfk13 C580Y, while C#321 carried the pfk13 G718S.All clones highlighted on the figure are referenced in the text.

Figure 4 -
Figure 4 -Clone persistence and the number of spatial clusters reached (> 1 sample and at least one spatial location, n=130, mean =287 days).Clones are colored by highly related clusters.C#268 carrying the pfk13 C580Y mutation is highlighted right on the persistence 80 th percentile (vertical line).The horizontal and vertical dashed lines represent the 95 th percentiles.All clones highlighted on the figure are ref in the text.

Figure 6 -
Figure 6 -Selection signals from isoRelate in long-lasting clones (sampled over three months) and short-lived clones sampled over two study periods (a-b) 2016/2017 and (c-d) 2020/2021.Dashed lines represent the threshold for the different selection signals investigated using a false discovery rate of 0.01.

Table 2 -Polymorphism which increased in frequency between 2016/2017 and 2020/2021
, and 2 uL DNA template as previously described by[69].The reactions were performed using the following conditions: initial denaturation at 95 °C for 15 minutes followed by 40 cycles at 95 °C for 15 seconds, 58 °C for 20 seconds, and 72 °C for 20 seconds; a melt curve starting at 95 °C for 2 minutes, 68 °C for 2 minutes, followed by increments of 0.2 °C from 68 °C to 85 °C for 0:05 seconds and a final step at 35 °C for 2 minutes.Copy number value was calculated using the 2^-Ct method[69].Means of pfpm2 and Pftub were calculated for 3D7 (a single copy control) using six replicates.Standard deviation should not be more than 25% including all algorithm[65], incorporating all SNPs that were called in ≥ 90% of samples and with minor allele