Impact of intensive control on malaria population genomics in Eastern Myanmar

TL;DR

Malaria control efforts in Eastern Myanmar have significantly reduced transmission and altered parasite genetic diversity.

Contribution

This study provides genomic evidence of malaria control impacts, including parasite population structure changes and localized transmission patterns.

Findings

Malaria control reduced transmission by 97% and decreased parasite genetic diversity.

Localized transmission of unique parasite genotypes was observed within 20 km.

Mass drug administration caused parasite founder effects and reduced effective population size.

Abstract

The malaria elimination program in Kayin State (Myanmar) utilizes malaria posts
for rapid detection and treatment together with mass drug administration (MDA) in high
transmission villages, and has reduced transmission by 97%. We examined the impact of
control on parasite genomic parameters, using 2270 genome sequenced Plasmodium falciparum
infections from 283 malaria posts, sampled over 58-months (2015 – 2020). Parasites
were genetically depauperate: 1726 single-genotype infections comprised 166 unique genomes
(≥90% IBD), while nine families (≥45% IBD) accounted for 62.5% of parasites
sampled. Parasite effective population size decreased over the study period, but there was
minimal change in artemisinin resistance alleles until 2020, when just one predominant
genotype (carrying kelch13-R561H) remained. We observed sustained localized transmission
of unique parasite genotypes revealing transmission chains: this resulted in positive
correlations in parasite relatedness for ≤20 km. MDA resulted in parasite founder
effects, providing genomic evidence for the efficacy of this control tool. These results
reveal changes in population structure driven by control, and rapid shifts in allele
frequency in a parasite population close to elimination.

Figures

C](#F1)). We removed 343 samples with > 20% genotypes missing. We then filtered the SNP calls following a “stringent” filtering method ^[34](#R34)^, to generate a final list of 447,435 high-quality, biallelic, core-genome located (defined in ^[33](#R33)^) SNPs. To analyze complexity of infection and population structure, we further removed SNPs that were genotyped in less than 50% of samples or with minor allele frequency (MAF) < 0.05.

C](#F3), Supplementary Figure 11) ^[6](#R6)^; (ii) the size of clonal clusters was not significantly different when comparing *kelch13* wildtype and mutant parasites ([Figure 3C](#F3), Supplementary Figure 12). The small effective population size of malaria parasite populations may contribute to the patterns observed, because selection is inefficient when population sizes are small and genetic drift is enhanced ^[29](#R29)^. The initial effective population size of malaria parasites in the Kayin State dataset was much smaller (Ne = 11.5 to 26.6) compared to SMRU clinics (15.5 – infinite) ([Fig

A](#F4)), where IBD cluster 1 carrying *kelch13*-R561H is at high frequency. Despite predominating in northern Hpapun Township since 2017, the *kelch13*-R561H didn’t spread into other areas of Kayin State. The rapid frequency increase of IBD cluster 1 in 2020 is consistent with bottlenecks and genetic drift in a *P. falciparum* population nearing elimination. Further sampling will determine whether this parasite genotype spreads further in Kayin State and elsewhere in Myanmar and Thailand. In contrast, Wasakul et al ^[22](#R22)^ describe a classic outbreak driven by a selective sweep in Laos,

A](#F6)). The relatedness between pre and post MDA parasites from these 3 HCPC units was significantly lower than observed between malaria parasites collected during the sampling time period from HCPC regions where MDA was not used ([Figure 6B](#F6)). Hence, MDA impacted parasite relatedness, consistent with efficient control of pre-MDA parasite genotypes, and post-MDA recolonization with unrelated genotypes.

Full text

Introduction

Regions of low malaria transmission intensity predominate in Southeast (SE) Asia and South America and are becoming increasingly common in Africa ^1^. A central challenge for malaria control is to develop efficient approaches to eliminate malaria from such regions. Rapid selection of drug-resistant parasites is a central concern for intensive malaria control programs. For example, kelch13 mutations conferring artemisinin resistance increased from 0–90% frequency in parasites collected from patients visiting Shoklo Malaria Research Unit (SMRU) clinics on the Thailand Myanmar border between 2003 – 2014 ^2^. Use of mass drug administration (MDA) is controversial for malaria treatment due to concerns about resistance. Prior use of chloroquine treated salt is thought to have accelerated selection of chloroquine resistance in the last century ^3^. However, MDA is effective for treating submicroscopic malaria infections which comprise the majority of infections in many SE Asian locations^4^, and are missed by passive malaria surveillance. Submicroscopic infections may be cured effectively by MDA because there are few parasites per patient, and so treatment is more likely to be completely successful than in high-parasitemia infections. It has therefore been argued that MDA is not likely to promote resistance spread in low prevalence regions like SE Asia ^5^. A previous paper examined the epidemiology of kelch13 haplotypes in Kayin State and showed limited changes of kelch13 artemisinin resistance alleles between 2013–2019 ^6^. Spread of resistance was much less rapid than occurred in the SMRU clinics (2000–2014); this is consistent with combined use of near-exhaustive coverage of communities with malaria posts and MDA imposing limited selection for ART-resistance.

Here, we examine genomic epidemiology of P. falciparum samples collected between November 2015 and August 2020 during malaria elimination efforts in Kayin State. Between 2014–2020, the Malaria Elimination Task Force (METF) targeted four townships (Myawaddy, Kawkareik, Hlaingbwe, and Hpapun) in Kayin State, eastern Myanmar, for malaria elimination ^7,8^. This was done using a combination of interventions: (i) 1475 village malaria posts (MPs) were opened, providing rapid diagnosis and malaria treatment (artemether–lumefantrine plus single low-dose primaquine); (ii) Mass Drug administration (MDA) (dihydroartemisinin–piperaquine [DHA-piperaquine] plus single dose primaquine once per month for 3 consecutive months) was used for 69 “hotspot” villages where malaria remained prevalent (>40% malaria and

20% Plasmodium falciparum) ^8^. The combination of these two approaches reduced P. falciparum cases by 97% from an incidence of 39 cases per 1000 person-years (May 2014–April 2015) ^8^ to 1 case per 1000 person-years (May 2019 to April 2020) ^9^.

Our central goal was to use genome sequence data to understand parasite transmission, population genomics and resistance evolution, and to use these data to inform future control efforts in elimination settings of SE Asia. Our key questions were to: (i) determine the distribution and stability of malaria populations within the Kayin State target region; (ii) to evaluate evidence for long-distance gene flow between Kayin State, Myanmar and other eastern SE Asian countries (Cambodia, Vietnam, Laos); (iii) document the origins of kelch13 resistance alleles; (iv) to determine the impact of MDA on parasite population structure; and (v) to evaluate appropriate genomic metrics for assessing transmission intensity.

Results

Extreme clonal expansion and inbreeding in a region under massive drug

selection

A total of 5014 dried blood spot (DBS) samples with geographic references from 413 MPs were collected between November 2015 and August 2020 as part of the METF malaria elimination program in eastern Myanmar, led by the SMRU. (Figure 1, Supplementary Figures 1 & 2). We processed and analyzed genome-wide sequencing data of 2270 DBS samples collected from 283 MPs (Table S1, Figure 1). After filtration of low-coverage samples and low-quality genotypes, the final Kayin State dataset contains 1927 P. falciparum samples with a set of 25,461 high-quality SNPs. 89.6% (1726/1927) of the samples were from single-genotype infections (Fws ≥0.90, Supplementary Figure 4).

To identify highly related individuals, we estimated pair-wise genetic relatedness (r, proportion of genome that was IBD) and used these to cluster samples with ≥90% of the genome IBD (r ≥0.90) (Figure 2). From the 1726 single-genotype Kayin samples, we identified 93 IBD clusters with unique genomes (2 to 229 samples per cluster), and 73 singletons, giving a total of 166 unique genomes. The Kayin population has an R_G_ of 0.10 (166 unique genomes from 1726 single-genotype infected samples), containing 14 large clonal expansion clusters (>30 samples per cluster, Figure 2). In contrast, the R_G_ ratios are much higher for the other SE Asian populations, with 0.72 (589/814) for SMRU clinics, 0.67 (123/184) for other Myanmar regions (not include Kayin), 0.54 (398/731) for Cambodia, 0.59 (102/173) for Viet Nam and 0.89 (80/90) for Laos. R_G_ ratios range from 0.94 to 0.98 for African parasite populations (Figure 2B).

Samples from the most common 10 IBD clusters account for 51.62% of the population, showing low parasite diversity in Kayin State. 151 out of 166 (91.0%) unique genomes were genetically related (r ≥0.25) to at least one other unique genome, indicating high levels of inbreeding in the Kayin State parasite population (Figure 3A). We identified 9 closely related families, that account for 110/166 unique genomes. Each individual in a family has a r ≥0.45 with at least another family member. Samples from these families account for 65.2% (1126/1726) of all the single-genotype infections in Kayin State. To further estimate genealogical relationships of individuals inside each family, we analyzed the distribution of chromosomal IBD segments between closely related family members (Figure 3B, Supplementary Figure 5). For two of the families (family 1 and family 7), we were able to identify both parents and F1 progeny based on their chromosomal recombination patterns, indicating extremely small parasite population size in Kayin State.

Localized transmission and regional stability of haplotypes

Genetic relatedness in space. Spatial groups of IBD clusters reflect direct or indirect transmission chains of malaria parasite clones, so are particularly informative for understanding transmission dynamics. The IBD clusters and closely related families show localized spatial distribution (Figure 4A, Supplementary Figure 6). For example, 86% samples from the largest IBD cluster (carrying kelch13 R561H) were collected to the north of Hpapun Township. The second largest IBD cluster (carrying kelch13 F446I) was found mainly in the center of Hpapun Township, while the third largest IBD cluster (carrying kelch13 P441L) was found in the west of the same township. Similarly, different IBD clusters carrying kelch13 wildtype alleles, show localized distribution (Figure 4A). Spatial correlograms confirm that parasite relatedness is positively correlated at distances ≤20 km (Figure 4B). We also observed significant negative correlations in relatedness between 27.5–90km. The negative correlations further indicate the spatial pattern of parasite relatedness. Not only are parasites that are geographically proximate more likely to be related – those that are geographically distal (27.5–90km) are less likely to be related.

Genetic relatedness in time. The length of time in which clonal lineages are sampled provides an indication of the frequency of outbreeding within malaria parasite populations. We detected clonal IBD clusters (r ≥ 0.90, contain ≥2 samples) that were sampled across the 56-month study period, as well as new genome haplotypes generated through recombination (Supplementary Figure 7). Out of the 93 IBD clusters, 9 lasted ≥ 36 months (3 years of sampling), while 34 lasted ≤ 6 months. The mean duration of IBD clusters in Kayin population was 13.8 (1 se = 1.4) months. Furthermore, the mean sampling duration of closely related parasite family members (r ≥0.45) was 48.7 (1 se = 3.0) months (Supplementary Figure 7), consistent with a low frequency of outbreeding in the population. We further analyzed the two parasite families for which both parents and progeny were identified. There were only 12 recombination events among the 109 samples (14 unique genomes, 41 months of family duration) from family 1, and 6 among the 31 samples (8 unique genomes, 48 months of family duration) from family 7 (Table S1, Supplementary Figure 5), consistent with low parasite outbreeding frequency.

We used a temporal correlogram to examine the number of days over which correlations in relatedness were observed (Figure 4B). This revealed positive correlations in relatedness between parasites sampled ≤ 170 days apart. To evaluate the relative impact of space and time on parasite relatedness, we divided the control region into 50 regions using hierarchical clustering on principal components (HCPC) (Supplementary Figure 8); 29 of these HCPC regions contained from 10–161 parasites. We examined relatedness in parasites sampled within and between HCPC regions for parasites sampled parasites sampled 1–12, 13–24, 25–36 and 37–48 months apart (Figure 4C). We observed significantly greater relatedness among parasites sampled from the same HCPC unit, relative to those sampled from different HCPC units. This remained significant for parasites sampled up to ≤36 months apart, demonstrating spatial stability of parasite populations within Hpapun Township (Figure 4C).

Long distance connectivity of parasite populations in SE Asia

We used parasite relatedness to measure connectivity within west SE Asian populations and between west and east SE Asia (Supplementary Figure 9). We detected a high level of gene flow between parasite populations from SMRU clinics and Kayin State: 38.3% of Kayin samples had >25% IBD (r > 0.25) with at least one sample from SMRU clinics, and 18.1% had >35% genome IBD (r

0.35). We identified two subpopulations, corresponding to west SE Asia (Kayin State, SMRU clinics and other Myanmar regions) and east SE Asia (Cambodia, Viet Nam and Laos) based on their genetic similarity and population structure (Supplementary Figure 10). Hence, we found no evidence for clonal transmission or recent recombination (r > 0.15) between west and east SE Asia.

The low connectivity between west and east SE Asia is further confirmed by the distribution of pfcrt mutations conferring piperaquine (PPQ) resistance. 52.19% of east SE Asian (Cambodia, Viet Nam and Laos) parasites carried PPQ-resistant pfcrt alleles between 2015 and 2018 (MalariaGEN Pf7 ^13^). These included T93S (23.55%), I218F (11.87%)^14,15^, H97Y (5.37%), F145I (8.20%), G353V (2.125)^16^ and G367C (1.08%)^17^. In contrast, these pfcrt mutations were absent from the Kayin dataset (Table S1) and from other west SE Asia regions (other Myanmar regions and west Thailand).

Genomic measures of parasite population size

We evaluated three genetic metrics (proportion of multiple-genotype infections, R_G_, and Ne) in both Kayin State and SMRU parasite populations, for assessing how control efforts impact parasite population size (Figure 5). The malaria incidence decreased significantly in both regions over studying time. The incidence decreased from 273.9 cases per person-year in 2001 to 22.4 in 2011 for SMRU clinics ^18^, and from 58.8 in 2016 to 1.0 in 2020 for Hpapun Township in the northern part of Kayin State, where >96% sequenced samples were collected ^8,9^ (Table S3).

Proportion of multi-infections. The proportion of multiple clone infections (10.4%, 201/1927) in Kayin samples was low compared to other P. falciparum populations (Table S3) and remained low throughout the year (range: 3.7–15.1%). The proportion of multiple clone infections decreased significantly in SMRU clinics, from 34.3% in 2001 to 4.3% in 2014 (p-value = 6.76e-05, R^2^ = 0.78), consistent with a prior analysis ^18^. However, this statistic showed no significant decline in Kayin (p-value = 0.61, R^2^ = 0.10, Figure 5B).

R_G_. R_G_ ratio is expected to be negatively correlated with the level of clonal expansion and positively correlated with transmission intensity ^18^. R_G_ ratios were lower in Kayin (0.10 to 0.30) between 2016 and 2020 than in SMRU clinics (range: 0.50–1.00) between 2001–2014. The R_G_ ratio decreased over time in SMRU clinics (p-value = 1.85e-03, R^2^ = 0.78), but not in Kayin (p-value = 0.54, R^2^ = 0.14).

Ne. We computed single-sample estimates of effective population size (Ne) using unique genomes from each population for each sampling year. Ne estimates were lower in Kayin (11.5 to 26.6) compared to SMRU clinics (15.5 to infinite). We detected significant reductions in Ne in both Kayin (p-value = 3.69e-03, R^2^ = 0.96) and SMRU clinics (p-value = 0.03, R^2^ = 0.40).

The evolution of pfkelch13 alleles

Impact of malaria elimination efforts on drug resistance. 61.32% of samples from Kayin State carried nonsynonymous SNP mutations in kelch13 (Supplementary Figure 11). The major mutant alleles were P441L (15.19%), F446I (15.01%), R561H (14.02%), and G449A (7.87%). Only 2.28% of Kayin samples carry C580Y. In comparison, in the adjacent SMRU clinics, C580Y was the dominant kelch13 mutant allele, reaching 71.05% allele frequency in 2014 (Table S2). There were 47 IBD clusters (r ≥ 0.90) carrying mutant kelch13 alleles and 46 IBD clusters carrying wild-type kelch13. We compared the size of IBD clusters carrying mutant kelch13 alleles with those carrying wild type kelch13 and found no significant difference (Figure 3C, Supplementary Figure 12). These results suggest that artemisinin selection was not the main driver for clonal expansion in Kayin State.

Clonal expansion of parasite carrying kelch13-R561H in 2020. Despite strong drug selection, frequencies of mutant kelch13 alleles remained stable between 2016 and 2019 ^6^ (Figure 3C, Supplementary Figure 12). However, in 2020 one of the kelch13 alleles - R561H - reached 74.2%. This clonal expansion results from near elimination of parasites from most areas in Kayin, with the exception of northern Hpapun Township where parasites carrying kelch13-R561H predominate (Figure 4D). 54.8% (40/73) of samples collected between January and August 2020 before the COVID-19 pandemic lockdowns were from one single malaria post (LH-0266B) (Table S1).

The change in kelch13 allele frequencies were reflected by changes in diversity in this gene and its flanking regions. We measured expected heterozygosity (He) to quantify diversity. He in kelch13 remained high between 2016–2019 (He = 0.78) but dropped to 0.41 in 2020. We observed parallel reductions in flanking region diversity with a drop from 0.48 to 0.19 between 2019 and 2020 (Supplementary Figure 11D).

Origins of kelch13 alleles. We reconstructed the haplotypes surrounding the kelch13 gene (100kb upstream and 100 kb downstream). We found a wide variety of kelch13 genetic backgrounds, with one or more unique haplotypes per resistance allele (Supplementary Figure 13 & 14). Two P441L, one F446I, one G449A and one R561H haplotypes had shared ancestry between Kayin State and SMRU clinics or other Myanmar regions. However, none of these alleles had high frequency in SMRU clinics or other Myanmar regions. Two F446I and one C580Y haplotypes were uniquely observed in Kayin State, consistent with local origin. For two Kelch13 resistance alleles (G449A, R561H), the wildtype kelch13 haplotypes on which these resistance mutations arose were sampled in the early 2000s in SMRU clinics. Of the three C580Y haplotypes identified in Kayin, two were also found in SMRU clinics, while one was unique to Kayin State. The two C580Y haplotypes shared with SMRU clinics were only found to the south of Hpapun Township and 60–120km north of the SMRU clinics (Supplementary Figure 15). None of these haplotypes shared IBD with east SE Asia C580Y haplotypes (Supplementary Figure 16).

Impact of mass drug usage on parasite population structure

We predicted that MDA would reduce relatedness of pre and post MDA parasite populations due to clearance of the local parasite population and replacement with new parasite genotypes post-MDA, resulting in founder effects. There were three HCPC units in which sufficient parasites (n ≥ 20) were sequenced both pre and post MDA (Figure 6A). The relatedness between pre and post MDA parasites from these 3 HCPC units was significantly lower than observed between malaria parasites collected during the sampling time period from HCPC regions where MDA was not used (Figure 6B). Hence, MDA impacted parasite relatedness, consistent with efficient control of pre-MDA parasite genotypes, and post-MDA recolonization with unrelated genotypes.

Discussion

The METF elimination efforts, combining community malaria posts and MDA, significantly decreased malaria case numbers in the target area (near the Myanmar-Thailand border) between 2014–2020 ^8,9^. We described the key results from genomic surveillance during these elimination efforts, including parasite transmission patterns, population diversity and genomics, evolution of drug resistance, and genomic impacts of MDA using over 2000 whole genome sequenced P. falciparum samples collected between Nov 2015 – Aug 2020.

Spatial and temporal structure of malaria populations

Parasite sequences from the METF study region reveal extremely high levels of inbreeding and low levels of genetic variation. We found 166 genotype clusters (≥90% of the genome IBD) among the 1726 single clone samples sequenced. Hence only 10% of parasite genomes sampled are unique (R_G_ = 0.1), and most infections show a clonal structure. In contrast, genetic richness is >0.94 in African populations sampled, and ranges from 0.54 to 0.89 in other SE Asian populations examined ^19^. Furthermore, 110 of the 166 unique genomes are distributed among 9 different families (r ≥ 0.45). Parasites within these families typically carry one or two kelch13 alleles, that mostly likely inherited from the parents. This is clearly the case in the two families with both parents identified - family 1 (F446I and wild-type kelch13) and family 7 (C580Y and wild-type kelch13). Hence recombination is rare in these populations and most infections are clonally related.

We observed strong spatial structure in the parasite population. This is particularly clear from the distribution of unique parasite genotypes. Such parasites result from self-fertilization of male and female gametes of the same genotype and allow spatial tracing of transmission networks. That these transmission networks are clustered in space is clear visually (Figure 4A), and statistically evident from autocorrelation analyses (Figure 4B), which reveal positive correlations in relatedness between genotypes for up to 20km. The local distribution of unique genotypes indicates (i) local transmission, and very few long-distance transmission events and (ii) reintroduction of circulating genotypes from asymptomatic carriers.

Unique parasite genotypes were long lived in this low transmission setting, with some IBD clusters observed over the complete 56-month study period. This is clear from (i) the temporal autocorrelation, which reveals positive correlation in relatedness between parasites collected 170 days apart; (ii) from the elevated relatedness observed in parasites collected from the same HCPC regions but up 3 years apart. The strong spatial and temporal sub-structure of parasite populations in Kayin State is comparable to that observed in Cambodia ^20^ and Guyana ^21^.

A clonal expansion of parasites carrying kelch13-R561H in

2020

As elimination approaches, genetic drift is expected to play an increasingly important role and expansions of parasite lineages may occur ^22^. We observed a sudden increase in the frequency of IBD cluster 1 parasites carrying kelch13-R561H in 2020. In this case, the increase of kelch13-R561H frequency resulted from elimination of malaria from all regions of Kayin other than northern Hpapun Township (Figure 4A), where IBD cluster 1 carrying kelch13-R561H is at high frequency. Despite predominating in northern Hpapun Township since 2017, the kelch13-R561H didn’t spread into other areas of Kayin State. The rapid frequency increase of IBD cluster 1 in 2020 is consistent with bottlenecks and genetic drift in a P. falciparum population nearing elimination. Further sampling will determine whether this parasite genotype spreads further in Kayin State and elsewhere in Myanmar and Thailand. In contrast, Wasakul et al ^22^ describe a classic outbreak driven by a selective sweep in Laos, where a lineage carrying kelch13-R539H (named LAA1) rose from a low frequency to replace the previously dominant KEL1/PLA1 (C580Y) population.

MDA impacts parasite population structure

Two control measures were used in Kayin State by METF: malaria posts (early diagnosis and community case management), and regional MDA in malaria hotspots ^8^. The combination of these approaches significantly decreased malaria incidence ^8,9^. Encouragingly, there was minimal evidence of selection for drug-resistance parasites through these elimination efforts which agrees with observations by Imwong et al ^23^, McLean et al ^24^ and Thu et al ^6^. At the genomic level, these combined control measures reduced parasite effective population size (Figure 5).

This study provided an opportunity to examine the impact of one component of this elimination strategy - MDA - on parasite population structure. MDA is expected to generate bottlenecks between pre and post MDA malaria populations, because reservoirs of asymptomatic malaria are removed. We therefore expect to see founder effects resulting from newly colonizing parasites and large divergence between pre and post MDA populations, when compared to populations with no MDA. Our power to detect an impact was limited because most parasites sequenced were collected after MDA: only three HCPC units had sufficient numbers of infection sampled both pre and post MDA to examine this hypothesis. Nevertheless, as predicted, we saw lower relatedness between pre and post treatment parasites in these three HCPC regions than in control regions where MDA was not implemented. These results provide genomic evidence for the effectiveness of MDA in removing local parasite populations through effective clearance of both asymptomatic and symptomatic parasite infections.

Genomic metrics for assessing transmission intensity

Genetic metrics, such as proportion of multiple clone infections can provide useful metrics for assessing control efforts ^18,25,26^. Such metrics are particularly useful in low transmission regions, where high prevalence of low-density asymptomatic infections complicates assessment of transmission using standard epidemiological methods. However, genetic metrics perform poorly when transmission levels are extremely low. In Senegal ^25,26^, complexity of infection provided an unreliable metric for evaluating transmission when transmission is incidence < 100 cases per 1000 person-years. Consistent with this, we observed that both R_G_ and proportion of Multiple infections worked poorly in Kayin where incidence ranged from 1 – 39 cases per 1000 person years. However, we found that another metric, effective population size (N_e_) calculated from LD between unlinked markers, showed significant decline in both the Kayin State parasite population and from SMRU clinics. We suggest that N_e_ may be a useful genomic indicator of transmission dynamics, particularly in parasite populations in which transmission has been reduced to near elimination levels. N_e_ is typically used in conservation biology to assess viability in endangered populations of animals and plants. Our results suggest that this metric may also have utility for assessing whether parasite populations are approaching local extinction.

No evidence for long-distance gene flow between Kayin State and Eastern SE Asian

countries

The current frontline treatments for P. falciparum malaria parasites have been failing in east SE Asia ^23,27,28^, due to the spread of the multidrug-resistance parasites carrying kelch13-C580Y mutation and plasmepsin 2 amplifications, named KEL1/PLA1. The KEL1/PLA1 lineage was first detected in Cambodia as the DHA-piperaquine was heavily used. When Cambodia withdraw DHA-piperaquine and adopted to artesunate–mefloquine, KEL1/PLA1 subgroups with acquired pfcrt mutations conferring piperaquine resistance rapidly spread to other ESEA countries, such as Laos and Vietnam ^27^. There is a concern that this parasite lineage will further spread to west SE Asia, which has the majority of malaria cases in SE Asia ^24^ and where DHA-piperaquine is the frontline treatment for P. falciparum. Two lines of evidence suggest minimal geneflow between east and west SE Asia. First, we did not detect pfcrt mutations associated with piperaquine resistance on Thailand-Myanmar border or in the Kayin State sampling sites. Second, examination of genome-wide IBD sharing among 3718 infections (1458 unique genomes) revealed no recent recombination or clonal transmission between east and west SE Asia (Supplementary Figure 9).

Origins of kelch13 resistance alleles

kelch13 mutations conferring artemisinin resistance are established in both east and west SE Asia. C580Y is the major mutation in regions other than northern Myanmar, where F446I predominates ^23^. In contrast, the dominant kelch13 mutations in the Kayin State include P441L, F446I, R561H, and G449A, depending on the location (Figure 4). While the majority of infections from nearby SMRU clinics carry C580Y (71.05% in 2014), the C580Y frequency in Kayin State was only 2.28%. The most widespread F446I haplotypes in Kayin State originated independently from the dominant F446I haplotype in northern Myanmar.

What factors lead to the patterns of artemisinin resistance evolution seenin Kayin? Longitudinal studies in both Cambodia and from SMRU clinics have revealed that multiple independent kelch13 mutations emerged and spread initially (soft selective sweeps). Single kelch13 genotypes (typically kelch13-C580Y) eventually outcompete other lineages leading to hard selective sweeps ^2,20,23,27^. In contrast, we found limited evidence that strong drug selection drove drug resistance evolution in Kayin State: (i) we found no significant increase in kelch13 mutant allele frequencies before 2020 (Figure 3C, Supplementary Figure 11) ^6^; (ii) the size of clonal clusters was not significantly different when comparing kelch13 wildtype and mutant parasites (Figure 3C, Supplementary Figure 12). The small effective population size of malaria parasite populations may contribute to the patterns observed, because selection is inefficient when population sizes are small and genetic drift is enhanced ^29^. The initial effective population size of malaria parasites in the Kayin State dataset was much smaller (Ne = 11.5 to 26.6) compared to SMRU clinics (15.5 – infinite) (Figure 2, Figure 5). Other factors that may also limit the impact of drug selection. Human population movements were more limited in Kayin compared to nearby SMRU clinics, especially in Northern Hpapun Township where human movement is limited by difficult terrain, the heavily militarized landscape, and a lack of year-round roads ^7,30^, which can hinder transmission of resistance alleles. Similarly, low levels of recombination in Kayin State limits the rate of formation of new multi-locus parasite genotypes. The small parasite population size, limited population movement, and minimal recombination enhance the role of genetic drift rather than selection in determining drug resistance evolution in the Kayin State region. Our results, and those from other studies ^23^ illustrate how genetic drift can result in rapid stochastic changes in parasite population genomics and drug resistance status in small parasite populations close to elimination.

This study had several limitations: (i) We analyzed malaria genomes collected from 2015 onwards. However, control efforts began earlier than this in 2014. Hence, we were unable to examine malaria population structure and diversity prior to initiation of control efforts. (ii) Use of finger prick blood samples and whole genome amplification resulted in bias towards sequencing high parasitemia infections. (iii) we were unable to score copy number variants, in genes such as Plasmepsin II/III, associated with piperaquine resistance from whole genome amplified DNA. However, the sustained sampling of a high proportion of blood spots collected over a 5-year period provides a unique dataset for examining impact of malaria control efforts on parasite population structure and resistance evolution.

Methods

Study area and sample origins

The samples for this analysis were collected during routine diagnosis and treatment efforts in Kayin State as part of the METF malaria elimination effort led by the Shoklo Malaria Research Unit (SMRU, based on the Thailand-Myanmar border) (Figure 1, Supplementary Figure 1& 2). This METF project was established in 2014 and utilized two primary P. falciparum-focused interventions: the establishment of a large network of community-based malaria diagnosis and treatment posts (MPs), and targeted MDA in communities determined to have a high prevalence of asymptomatic P. falciparum infections. The MPs were stocked with filter papers (Whatman 3mm blotting paper) and were asked to collect dried blood spots (DBSs) from finger prick blood samples for patients with rapid diagnostic test (RDT) confirmed P. falciparum infection. Each DBS sample is linked to the MP from which it originated, and all MPs have geographic references (latitude and longitude). 5014 DBS samples were collected between November 2015 and August 2020 (Table S1). The DBS samples were then transported to SMRU and subsequently shipped to the Texas Biomedical Research Institute (in the U.S.A.) for molecular analyses.

Sequencing library preparation

We extracted DNA from the dried blood spots and enriched parasite genomes using selective whole genome amplification (sWGA) to following Li et al ^31^ and Oyola et al ^32^. We extracted and purified genomic DNA using QIAamp^®^ 96 DNA Blood Kit or QIAamp DNA Mini Kit, following the instruction manual for DNA isolation from dried blood spots. The DNA was then eluted in 100ul of 10mM Tris-HCl (pH 8.0–8.5) buffer. We used real-time quantitative PCR (qPCR) to estimate the numbers of parasite genomes in each DNA sample as described in Li et al ^31^.

For samples with more than 200 copies of parasite genome per ul, we used selective whole genome amplification (sWGA) to enrich parasite DNA. sWGA reactions were performed following Oyola et al ^32^. Each 25 μl reaction contained at least 1000 copies of Plasmodium DNA, 1× BSA (New England Biolabs), 1 mM dNTPs (New England Biolabs), 3.5 μM of each amplification primer, 1× Phi29 reaction buffer (New England Biolabs), and 15 units of Phi29 polymerase (New England Biolabs). We used a PCR machine (SimpliAmp, Applied Biosystems) programmed to run a “stepdown” protocol: 35 °C for 10 min, 34 °C for 10 min, 33 °C for 10 min, 32 °C for 10 min, 31 °C for 10 min, 30 °C for 6 h then heating at 65 °C for 10 min to inactivate the enzymes prior to cooling to 4 °C. Sample were cleaned with AMPure XP Beads (Beckman Coulter), at a 1:1 ratio.

We used the Quant-iT^™^ PicoGreen^®^ Assay (Invitrogen) to determine the total amount of sWGA product, and quantified the proportion of Plasmodium DNA by qPCR. Only sWGA products with more than 50% DNA from Plasmodium were used for further library preparation and Illumina sequencing. We constructed PCR-free next generation sequencing libraries using 300ng sWGA product following the KAPA HyperPlus Kit protocol. All libraries were sequenced to an average coverage of 60× using Illumina Hiseq X or Novaseq sequencers.

Whole-genome sequencing data generation

We individually mapped whole-genome sequencing reads for each library against the P. falciparum 3D7 reference genome (PlasmoDB, release 46) and human GRCh38 reference genome, using the alignment algorithm BWA mem (http://bio-bwa.sourceforge.net/) under the default parameters. The resulting alignments were then converted to SAM format, sorted to BAM format, and deduplicated using picard tools v2.0.1 (http://broadinstitute.github.io/picard/). Reads mapping to the human genome were discarded before genotyping.

We used Genome Analysis Toolkit GATK v3.7 (https://software.broadinstitute.org/gatk/) to recalibrate the base quality score based on a set of verified known variants ^33^. We called variants for each sample using HaplotypeCaller, and calls from every 100 samples were merged using CombineGVCFs with default parameters. Variants were further called at all sample-level using GenotypeGVCFs, with parameters: --max_alternate_alleles 6 --variant_index_type LINEAR --variant_index_parameter 128000 --sample_ploidy 2 -nt 20.

The recalibrated variant quality scores (VQSR) were calculated by comparing the raw variant distribution with the known and verified Plasmodium variant dataset. SNPs and indes with VQSR less than 1 or located outside of the core genome regions (21 Mb, defined in ^33^) were removed from further analysis. Samples with less than 50% of the core genome callable were also excluded from further analysis. Only biallelic SNPs that pass all the quality filter were used, unless otherwise specified. The final variants in VCF format were annotated at functional effect to genes and proteins using snpEff v4.3 (https://pcingola.github.io/SnpEff/) with 3D7 (PlasmoDB, release 46) as reference.

We initially identified 1,302,006 single-nucleotide polymorphisms (SNPs) and 703,138 indels (Figure 1C). We removed 343 samples with > 20% genotypes missing. We then filtered the SNP calls following a “stringent” filtering method ^34^, to generate a final list of 447,435 high-quality, biallelic, core-genome located (defined in ^33^) SNPs. To analyze complexity of infection and population structure, we further removed SNPs that were genotyped in less than 50% of samples or with minor allele frequency (MAF) < 0.05.

Complexity of infection

We measured multiplicity of P. falciparum infections using the within-infection Fws fixation index ^35^. Samples with FWS > 0.9 were assumed to come from single-genotype infections for samples from Kayin State. Allele frequencies across the genome were plotted and manually inspected to detect further possible complex infections.

Relationships among parasite genotypes

We used relatedness - r, defined as the fraction of the genome that is identical-by-descent (IBD) between a pair of individuals ^36,37^ - to estimate parasite relationships. Based on the distribution of relatedness among F1 progeny from malaria parasite genetic crosses (Supplementary Figure 3), we assume that parasites are genetically related if ≥ 25% of their genome is identical (r ≥ 0.25); parasites are closely related (such as siblings or parent and progeny) if their relatedness is greater than 45% (r ≥ 0.45). We considered samples to be clonal if their relatedness is over 90% (r ≥ 0.90). We visualized relatedness among samples using the R package pheatmap and the Cytoscape software. We also examined the recombination patterns between closely related parasites and plotted shared IBD regions between estimated parents and progeny using karyoploteR.

Surveillance of kelch13 haplotypes

We extracted SNPs distributed within 100 kb upstream and 100 kb downstream of the kelch13 gene. We measured expected heterozygosity (He) at the kelch13 locus by treating kelch13 as a single locus with multiple alleles. We also measured He over the 200kb kelch13 haplotype region. To compare the relationships between different kelch13 haplotypes, we measured pairwise IBD sharing among all kelch13 haplotypes. We assume that haplotypes with IBD sharing ≥ 0.90 originated from the same mutational event; that when 0.35 ≤ IBD < 0.90, there was a one least recombination event to break the original haplotype; and when IBD < 0.35, these haplotypes have emerged independently.

Comparisons of malaria parasite populations

We compared the Kayin State parasite population with other world-wide malaria parasite populations (Figure 1, Table S2). The SMRU clinics are located around Mae Sot, in Tak Province along the international Thailand-Myanmar border. We used “other Myanmar” to represent sampling sites in Myanmar but not from Kayin State. West SE Asia population includes samples from Kayin State, SMRU clinics and other Myanmar regions, while east SE Asia population includes Cambodia, Viet Nam, and Laos.

We merged raw SNP genotypes from the Kayin dataset with those from MalariaGEN P. falciparum Community Project ^19^ (release 6). We performed “stringent” filtration as described above, and selected loci with minor allele frequency > 0.05. We calculated genetic richness $[eqn]$ ^38,39^ to quantify the richness of clonal parasites in each population, where $[eqn]$ is the number of unique genomes, and $[eqn]$ is the total number of single genotype infected samples. For samples with relatedness > 0.9, only one representative sample per population with the highest genotype rate was selected and used for further analysis (Table S2). We pruned SNPs for linkage disequilibrium (LD) and generated a pairwise genetic distance matrix using PLINK with default parameters. We conducted principal component analyses (PCA) and ADMIXTURE analyses based on the pruned genotypes and distance matrix. We measured the proportion of pairs IBD across the genome within populations following the scripts in isoRelate ^40^. We estimated effective population size (Ne) based on patterns of LD at unlinked loci, using methods implemented in NeEstimator v2.0 ^41^.

Statistical analysis

All statistical analysis was performed using R version 4.1.3. For pairwise comparisons between groups, we used Welch Two Sample T-test. We measured correlations between parasite genetic relatedness and geographic distance or time using the Mantel statistic using the mantel function in the vegan package. p<0.05 was considered statistically significant. We used hierarchical clustering on principal components (HCPC) following scripts in FactoMineR ^42^ to divide the 283 malaria posts with samples sequenced into 50 HCPC regions based on latitude and longitude. We then compared parasite relatedness within and between HCPC regions for parasites collected in the same year, between 1–2, 2–3 and 3–4 years apart. We compared relatedness between parasites collected from HCPC regions 6 months before and 6 months after MDA. As controls, we examined relatedness of parasites collected from HCPC regions where MDA was not used during the same time windows.

Supplementary Material

Additional Information

Supplementary Information is available for this paper.

Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

• 2.Supplementarytables1306.11.2025.xlsx

• SupplementaryFigures.docx