Routine collections of dried blood spots (DBS) were conducted at selected sentinel site health centers (2018-2020) in the malaria-endemic provinces Binh Thuan, Binh Phuoc, Gia Lai, Khanh Hoa, Kon Tum, Lam Dong, Ninh Thuan, and Quang Tri to monitor uncomplicated malaria (P. falciparum and P. vivax) in Vietnam. In each province at least 1 sentinel site was chosen supporting several surrounding communities. In provinces with many malaria cases (such as Gia Lai, Khanh Hoa, Lam Dong), samples were collected at 2 sentinel sites, often one in a district hospital and one in a commune. Distribution of this collection depended heavily on the caseload in each province, with some provinces reporting a very low number of malaria cases (Table 2). Individuals between 7–60 years old living around the sentinel site clinic, visiting the clinic with suspected malaria were invited to take part in the surveillance. In addition, 11 samples from previous routine collections in 2015–2017 from Dak Nong (n = 3), Dak Lak (n = 3) and Khanh Hoa (n = 4) were collected for whole-genome sequencing (WGS) analysis for the barcode design, as there was no prior genomic data available from these sites, and one additional sample from Binh Phuoc included. DBS were collected from confirmed malaria cases (positive by microscopy and/or RDT (Rapid Diagnostic Test)) when informed consent had been given. Details of the P. vivax cases are described elsewhere (Kattenberg et al., 2022).

A total of 1010 DBS samples were available from sentinel site collections (Figure 1) and DNA was extracted using one 5 mm punch (2019 samples) or three 5 mm punches (all other DBS samples) of DBS using the QIAamp DNA Blood Mini Kit (Qiagen) with a final elution in 200 μL elution buffer from the kit. VarATS and 18S-qPCR (Hofmann et al., 2015) were performed to confirm diagnosis (Ct < 42) and 407 samples were included for library preparation (Figure 1).

To design a SNP barcode with in-country resolution in Vietnam, raw WGS data (FASTQ files) of P. falciparum isolates from Vietnam (2015-2017) generated in-house (n = 11) was combined with a VCF file previously generated using Vietnamese and Cambodian P. falciparum genomes (n = 262) from our lab that were externally sequenced as part of MalariaGEN SpotMalaria Project (https://www.malariagen.net/projects/SpotMalaria as described in (16). Newly generated sequencing data for the 11 isolates was deposited in SRA under project number PRJNA1156587, accession numbers are listed in Supplementary Table S1.

In house generated FASTQ files were aligned to the 3D7 reference genome version 41 from PlasmoDB [https://plasmodb.org/plasmo/app/downloads/release-41/Pfalciparum3D7/] and variants were called using GATK4 HaplotypeCaller and merged using bcftools (Danecek et al., 2021; Poplin et al., 2017). Variants were filtered on DP > 50, keeping only SNPs with at least 1 alternate allele (i.e., removing SNPs that were homozygous reference in all samples or missing in all samples from the SpotMalaria dataset), resulting in 53,362 high quality SNPs.

Three samples were identified as outliers in principal component analysis (PCA) and removed from further analysis (SPT25139, SPT25147 and SPT25151). The remaining samples were used in a pipeline to select SNPs for the within-country barcode, in a similar procedure as the selection of SNPs for the P. vivax barcode for Vietnam (Kattenberg et al., 2022). Briefly, LD-pruning was performed iteratively using 500 bp windows, removing SNPs with pairwise LD > 0.2 via scikit-allel (Miles et al., 2024). Remaining SNPs were filtered for minor allele frequency (MAF) > 10%, and selectively neutral SNPs were retained with |Tajima’s D| > 0.5 (calculated in 500 bp windows), resulting in 638 SNPs. Subsequently, the contributions of the SNPs to geospatial genetic clusters were determined using discriminant analysis of principal components (DAPC) with province populations using adegenet package in R (Jombart et al., 2010). We conducted a parameter sweep for the DAPC varying the number of discriminant analysis (DA) clusters between 30-40, number of principal components between 25-40 and number of axes in the DA between 4-6 resulting in a list of 638 potential SNPs. Per chromosome, we selected 4-7 SNPs with loading contributions in the DAPC > 0.003 with the highest count in the DAPC parameter sweep to add to the design.

The Illumina Concierge team (Illumina, San Diego, United States) used DesignStudio software with the Pf3D7 reference genome to design amplicons for the new barcode and other selected markers. Amplicon design was successful for 46 SNPs. Based on recent literature research, we identified several new genes of potential interest for artemisinin resistance. Therefore, compared to our previous Pf Peru AmpliSeq panel, we targeted the additional genes PfKIC7 (PF3D7_0813000) and PfEps15/Formin 2 (PF3D7_1025000), both associated with in vitro resistance to ART (Henrici et al., 2020; Birnbaum et al., 2020; Kattenberg et al., 2023a). In addition, the targets Pfmdt (PF3D7_0516500) and PftetQ (PF3D7_1235400) were added, where CNVs in both genes have been associated with doxycycline resistance (Briolant et al., 2010; Gaillard et al., 2015). We no longer targeted the cytB gene and 23S rRNA genes, as the mitochondrial and apicoplast templates are more predominant than nuclear targets and are overrepresented in the resulting sequences (Kattenberg et al., 2023a). In the amplicon design, full resistance associated genes were targeted (Table 1), as well as six KEL1-lineage markers (incl. RAD5 and arps10), variable regions D1 and D2 in PfamaI, and Pfhrp2 and Pfhrp3. The final designed panel for the Pf AmpliSeq v2 Vietnam included 265 amplicons, ranging in size from 115–379 bp (Supplementary Table S2).

Library preparation was done for each sample and control using custom primers (Pf AmpliSeq v2 Vietnam) and an AmpliSeq Library PLUS kit (Illumina), as previously described in a step-by-step laboratory protocol (Kattenberg et al., 2023b). Briefly, libraries were indexed (Illumina), purified using AMPure XP magnetic beads (Beckman Coulter) and quantified using Qubit v3 High sensitivity DNA kit (Invitrogen). Libraries were diluted to 2 nM with low Tris-EDTA buffer, and pooled (pooling 96–144 libraries per sequencing run). Denatured library pool (7 pM) was loaded on a MiSeq system (Illumina) for 2 × 300 paired-end sequencing (Miseq Reagent Kit v3, Illumina) with 5% PhiX spike-in (Illumina).

A set of DNA controls provided through BEI Resources was used to validate the detection resistant markers in PfK13 (MRA-1241, MRA-1251, MRA-1255). The specificity of the assay was validated using extracted human malaria-negative blood, 3D7 reference wild type, D10 drug sensitive lab strain, Dd2 lab strain and mock Pf-Pv co-infections. Sensitivity of the assay was validated using a 3D7 dilution from 2,450 to 1,000 parasites/µL, followed by 1:10 dilutions up to 1 parasite/µL. Drug resistance variant calling, and population structuring was validated using 23 previously characterized samples from Cambodia and Vietnam that were genotyped using a TruSeq Custom Assay and SpotMalaria genetic report card (Rovira-Vallbona et al., 2020; Rovira-Vallbona et al., 2023).

FASTQ files were processed using an in-house bioinformatic pipeline published on GitHub (https://github.com/Ekattenberg/Plasmodium-AmpliSeq-Pipeline). Briefly, demultiplexed reads were trimmed (Trimmomatic) to remove adapter sequences and low-quality bases (Bolger et al., 2014). Trimmed reads were then aligned to the reference genome 3D7 (version 44, PlasmoDB) using Burrows-Wheeler aligner (Li and Durbin, 2009). Alignment statistics were generated with Picard’s tools and variants were called using GATK4 haplotypecaller generating a jointly-called VCF file with variants (SNPs and INDELs) detected in the targeted regions (Poplin et al., 2017; Picard Toolkit, 2019). Variants were hard filtered (QUAL>30, overall DP > 500, MQ > 50, QD > 1.0, SOR<4, ReadPosRankSum > −10.0, GT depth >5) and annotated with SnpEff (v4.3T) (Cingolani et al., 2012), resulting in 2,215 high quality genotypes (incl. all variant types, e.g., SNPs and indels). Per locus filtered depth of coverage (format field DP in the VCF) was used to calculate the median depth of all loci per sample or per amplicon. Aligned coverage was calculated as the number of bases passed filter divided by the number of bases (67,503 bp) targeted in the Pf AmpliSeq v2 Vietnam.

Samples with good quality sequences that passed our inclusion criteria (average aligned coverage >10; and missing genotype calls <37.5%) were selected for further analysis. In total 354 samples were analyzed coming from Gia Lai (231/354), Binh Thuan (37/354), Khanh Hoa (34/354), Binh Phuoc (19/354), Quang Tri (11/354), Ninh Thuan (9/354), Kon Tum (7/354), (Lam Dong (6/354) (Table 2).

We created a list of variants of interest (Supplementary Table S3) that included variants validated for resistance according to the WHO, WARNN meta-analysis and mutations in targeted genes reported in the literature as potentially associated with P. falciparum antimalarial resistance (World Health Organization, 2018; Rovira-Vallbona et al., 2023; Kattenberg et al., 2023a; Das et al., 2013). The list was supplemented with non-synonymous variants detected in the target genes that contributed to the variation in the DAPC. Haplotypes were created by combining genotypes of variants of interest.

Genetic diversity, expressed as expected heterozygosity (He), was calculated using polymorphic barcode (46 SNPs) genotypes from the VCF with the adegenet package in R (Jombart, 2008; Jombart and Ahmed, 2011). To compare the median of the genetic diversity parameters across different provinces, Kruskal-Wallis rank sum test was performed. For pairwise comparisons between groups, Wilcoxon test with Benjamini-Hochberg correction for multiple testing was used. Statistical tests were performed with the R package stats. P-values <0.05 were considered significant. PCA and DAPC with cross-validation was performed to infer population structure based on haplotypes across districts and years (Jombart et al., 2010). Associated allele loadings for the first four components in the DAPC were determined. Genetic differentiation, expressed as fixation index (F_ST), was calculated using barcode SNP genotypes using the R package hierfstat (Goudet et al., 2022). Within-host infection complexity was assessed using within-sample F-statistic (Fws) with the R package moimix using all biallelic SNPs (n = 1,276) detected by the assay (Lee and Bahlo, 2016). Fws ≥0.95 was considered a proxy for a monoclonal infection as in Auburn et al. (2012).

To measure pairwise identity-by-descent (IBD) between samples, PED and MAP file formats were generated from VCF data using VCFtools (Danecek et al., 2011). The level of IBD-sharing was calculated employing the isoRelate package in R (Henden et al., 2018). We used all biallelic SNPs (n = 1,276) identified by the Pf AmpliSeq v2 Vietnam assay, applying filtering criteria of MAF = 0.001 and SNP and individual missingness thresholds (0.6 and 0.3, respectively), resulting in 952 SNPs for subsequent IBD analysis. Furthermore, IBD thresholds were set to include a minimum of 15 SNPs per segment and a segment length of 3,000 bp, aimed at mitigating false-positive calls using an error of 0.001. Network plots (at thresholds of 90%) were created using the igraph package in R to visualize the relatedness between samples and connectivity between districts as described before (Kattenberg et al., 2022; Kattenberg et al., 2024; Csárdi et al., 2024).

The read depth was extracted for each amplicon in each individual BAM file using samtools (selecting only reads with MQ ≥ 60) (Danecek et al., 2021). For this analysis we included only samples with mean read depth between 500 and 9,000 (306/354 samples). The read depths were normalized (creating values between 0 and 1) using the normalize function of the heatmaply package in R (Galili et al., 2018). Samples were clustered based on the amplicons of the target genes (Pfpm2, Pfmdr1, Pfhrp2, or Pfhrp3) and reference amplicons with little variation in depth (Supplementary Figure S1) and where possible on the same chromosome as the target gene. Hierarchical clustering was performed using the heatmaply package in R, based on Euclidean distances between read depth profiles of samples and the complete clustering method. Clustering patterns and normalized read depth were visually inspected in the resulting heatmaps. Individual sample Pfpm2 amplifications were determined for samples belonging to two clusters that contained qPCR confirmed CNV in this gene and had clear increase in read depth compared to the reference amplicons.

The sample collections in Vietnam were approved by local ethical review boards at NIMPE and by agreement of the Ministry of Health in Vietnam. Individuals were included in this study only if they willingly signed informed consent that included an opt-in future-use clause. Secondary use of all samples was approved through the Institutional Review Board of the Institute of Tropical Medicine Antwerp (reference 1417/20).

We designed a 46-SNP Pf Vietnam Barcode with in-country resolution, able to differentiate distinct P. falciparum populations in Vietnam based on the variability of allele frequencies between provinces. This barcode was developed using P. falciparum genomes from Vietnam and Cambodia (n = 273), where neutral SNPs were selected based on their contribution to geographically distinct genetic clusters in DAPC. The minor allele frequencies of all SNPs in the barcode varied across a range of 0.11–0.49, with a median value of 0.36 (Supplementary Table S4). The Pf Vietnam Barcode was able to differentiate the samples by province in the PCA analysis, although clusters where less pronounced than with all biallelic SNPs (n = 44,553 SNPs) in the WGS data (Figure 2). With all biallelic SNPs, the strongest differentiation was observed with the cluster of samples from Dak Nong and some samples from Dak Lak (along PC1 in Figure 2C), which also separated with the barcode SNPs (along PC2 in Figure 2A). In addition, strong differentiation was seen with part of the samples from Gia Lai with both the barcode and all biallelic SNPs (Figure 2A, E, F). The samples from Khanh Hoa and Ninh Thuan clustered close together but could still be differentiated at the resolution of the SNP barcode (Figures 2A, B).

We assessed the possibility of Pfhrp2 and Pfhrp3 deletions by investigating the read depth. Both genes amplified very well in the assay, with general high read depth at the amplicons targeting these genes. We were able to detect Pfhrp2 deletions in the laboratory Dd2 control strain, but only in very few study samples (n = 5) we observed a similar lack of Pfhrp2 amplifications, indicating a deletion in this gene (Supplementary Figure S4). Similarly, Pfhrp3 amplicons amplified well in all samples, with only very few samples where all or some of the amplicons did not amplify (n = 4). In summary, deletions in Pfhrp2 and Pfhrp3 were rare in Vietnam.

Most samples (255/354, 72%) were single clone infections, based on Fws determined using all biallelic variants. There was no difference in the proportion of single clone infections by province (range 59%–90%, Χ² p-value = 0.19). Binh Thuan province had the highest proportion of multiple clone infections (40.5%), followed by Gia Lai province (29.9%).

Median genetic diversity (He) was lowest in Kon Tum (median He = 0), followed by Binh Thuan and Gia Lai province (median He = 0.054 for both, Figure 8A). Diversity was highest in Lam Dong and Quang Tri (median He = 0.28 and 0.30, respectively, Figure 8A), and significantly higher in Quang Tri compared to other provinces, except Ninh Thuan and Lam Dong (p< 0.05 pairwise Wilcoxon BH-corrected for multiple comparisons). While the burden of malaria (number of cases) in Gia Lai province is much higher than other provinces (Table 2), the genetic diversity was one of the lowest. Genetic differentiation (F_ST) was low between Gia Lai province and Binh Phuoc, and Binh Phuoc and Kon Tum (Figure 8B). Parasites from Khanh Hoa showed high genetic differentiation with other provinces. Binh Thuan is different from other Coastal provinces, and more similar to Binh Phuoc and Gia Lai populations (Figure 8B).

Population structure between provinces in Vietnam was further explored using DAPC using all biallelic SNPs for higher resolution analysis (Figure 9). DAPC is a model in which genetic variation is partitioned into a between-group and a within-group component, and yields synthetic variables that maximize the first while minimize the second. In other words, DAPC attempts to summarize the genetic differentiation between groups, while slightly overlooking within-group variation. This is similar to PCA but uses a statistical approach to maximize between-group variation. North and South-Central Coast Provinces Khanh Hoa Ninh Thuan and Quang Tri clustered away from Central Highland provinces, Gia Lai and Kon Tum, and Binh Phuoc (Southeast region), while South Central Coast Province Binh Thuan clustered close to these provinces (Figures 9A, B). The population from Khanh Hoa clustered furthest away from other populations, except from Lam Dong, consistent with the F_ST results with the barcode (Figure 9A). While Ninh Thuan Province neighbors Binh Thuan, Lam Dong and Khanh Hoa provinces, the P. falciparum parasites in this province clustered closer to isolates from more northern provinces Quang Tri and Kon Tum along PC2 (Figure 9A), while it clustered away from all provinces in PC4 (Figure 9B). Upon closer investigation of the main cluster of Central Highland provinces and Binh Phuoc and Binh Thuan, PC1 indicated a clustering pattern from northern to southern provinces, and Kon Tum separated from the other provinces along PC2 (Figure 9C). Taken together, these results indicate that the parasite populations in the provinces of Vietnam can be differentiated at distinct levels in the DAPC analysis with marked geographical clustering of isolates in Vietnam.

Clusters of highly related (>90% IDB-sharing) P. falciparum isolates were seen especially in Gia Lai and Binh Thuan province and a separate cluster in Khanh Hoa province (Figure 10A). The clonal structure of parasites in these provinces explains the observed low diversity, despite the high number of samples from these provinces. The clonal clusters in Gia Lai and Binh Thuan all belonged to the KEL1-lineage (genotyped at several loci on chromosomes 13 and 14, identifying 16 haplotypes, Supplementary Table S7) and contained the C580Y mutation in PfK13, and in some clusters also the PPQ-resistance conferring mutation N326S in Pfcrt (Figures 10B, C). Some samples from Binh Phuoc were also related to these clusters (Figure 10A). The cluster of samples from Khanh Hoa were predominantly ART- and PPQ-sensitive, while it also contained the 4 samples with the PfK13 P553L mutation, indicating this mutation emerged in the distinct local parasite population, which also shared ancestry to the KEL1-lineage, but without mutations C580Y and N326S (Figures 10B, C). In pairwise analysis of IBD between isolates, high relatedness was seen especially between isolates within a province and less so between provinces, indicating little connectivity between the P. falciparum populations between provinces (Figure 10D).

Vietnam is ambitiously targeting the elimination of malaria by 2030, a goal that necessitates rigorous monitoring of the disease to inform strategic interventions. The adoption of innovative tools like the AmpliSeq Pf Vietnam v2 assay, can significantly contribute to understanding the genetic diversity and drug resistance profiles of parasite populations in a region with a high level of multi-drug resistant P. falciparum (World Health Organization, 2022; Jacob et al., 2021). In a period of steep decline of malaria incidence in the country (2018-2020), we observed an increase in resistance to both artemisinin and its partner drug piperaquine in more provinces compared to previous studies, indicating a rapid spread. Nearly 80% of all samples analyzed (281/354) carried the C580Y mutation in the PfK13 gene, while our earlier studies conducted in similar regions before 2018 detected validated PfK13 mutations in approximately half the samples (Rovira-Vallbona et al., 2023). Other recent studies confirm our observations of high proportion of ART-R PfK13 mutations, predominantly C580Y (Verschuuren et al., 2024; Thu et al., 2023). This highlights the critical need for ongoing surveillance and adaptive strategies to effectively combat malaria.

High proportions of parasites resistant to DHA-PPQ were observed in the provinces of Binh Phuoc, Binh Thuan, Gia Lai, and Kon Tum. Observations of treatment failures in several provinces in the country prompted a change in treatment guidelines by the NMCP from DHA-PPQ followed by PA as the first-line treatment for P. falciparum malaria in these provinces in the same period: in 2019 for Dak Nong and Binh Phuoc, followed by Gia Lai, Khanh Hoa, Phu Yen, and Dak Lak Province in 2020. In addition, we show for the first time high levels of DHA-PPQ resistance in Binh Thuan Province, observed as early as 2018. Binh Thuan province has set an ambitious goal to eliminate malaria by 2026, and the decline in reported cases suggests progress, but the spread of resistant lineages from neighboring provinces may pose challenges. Fortunately, in 2023, PA was adopted as the first-line treatment for both P. falciparum and P. vivax infections in all provinces in Vietnam, building on the success of PA in prior TES conducted in 2017–2018 and its proven efficacy in other provinces (Manh et al., 2021; Quang Bui et al., 2020; Minister of Health Vietnam, 2023). In contrast, the parasite population in Quang Tri has remained sensitive to ART, PPQ, and CQ, setting it apart from other regions. While Quang Tri province is part of the Central Highlands, and one of the provinces with the highest malaria burden in the early 2000s, malaria transmission has decreased rapidly much earlier than in other provinces in the Central Highlands, and remaining pockets are in remote areas near the border with neighboring Laos that also has little ART-R in that region (Rovira-Vallbona et al., 2023; Pongvongsa et al., 2012).

The observed Pfpm2 gene amplifications and the persistence and spread of Pfcrt mutations present a complex picture of piperaquine resistance dynamics. While early studies were showing an increase in Pfpm2 gene amplifications in Vietnam (from 3% in 2000 to 11% in 2016), our study along with other recent reports indicate a stabilization and potential decline of this marker (Rovira-Vallbona et al., 2023; Verschuuren et al., 2024). However, while a recent preprint reports the decline of resistance to both DHA and PPQ in P. falciparum parasites in the GMS, this interpretation is based on the decrease of Pfpm2 gene amplifications alone and Pfcrt mutations associated with PPQ were not investigated (Verschuuren et al., 2024). With the extensive panel of markers in the Pf AmpliSeq v2 Vietnam assay, we detected a significant increase of PPQ resistance in 2019 and 2020 compared to previous years, primarily driven by the emerging Pfcrt variants, with resistance observed in over 50% of the analyzed isolates. Recent studies suggest that Pfcrt mutations may confer a more stable resistance phenotype with lower fitness costs than Pfpm2 amplifications, which could explain their selection in these parasite lineages in the face of PPQ pressure (Hastings and Donnelly, 2005; Kim et al., 2019). The shift from Pfpm2 amplifications to Pfcrt mutations observed here could therefore indicate an evolutionary adaptation, allowing parasites to survive better. In addition, the resistant parasite populations observed in Binh Phuoc, Binh Thuan, Gia Lai and Kon Tum provinces exhibited a clonal nature and low genetic diversity (median He ≤ 0.11), indicating that the highly resistant parasite lineage recently emerged in the KEL1-PLA1 genetic background. Nevertheless, the PPQ susceptibility profiles of parasites with only PPQ-R Pfcrt variants, without Pfpm2 amplifications, requires further phenotypic characterization, as this was not possible with the DBS samples collected in this study. Previous clinical studies in Vietnam included samples with both mutations, whereas in vitro studies and studies in other regions have indicated that Pfcrt mutations alone can be sufficient to confer PPQ resistance (Van Der Pluijm et al., 2019; Manh et al., 2021; Ross et al., 2018; Boonyalai et al., 2022; Florimond et al., 2024; Silva et al., 2020).

Beyond ART and PPQ resistance, this study shows high levels of multidrug resistance, with 6/8 provinces reporting very high SP resistance levels with >70% of the parasites carrying double, triple or quadruple mutations in Pfdhps and Pfdhps. The continued use of CQ for P. vivax malaria treatment, together with PPQ pressure, shapes Pfcrt haplotypes, preventing a reversion of P. falciparum parasites to CQ sensitivity as seen in Africa since the introduction of ACTs (Asare et al., 2021; Laufer et al., 2010; Balikagala et al., 2020; Kateera et al., 2016). In most provinces, besides Lam Dong, Ninh Thuan and Quang Tri, CQ-R Pfcrt CVIET and CVIDT haplotypes were reported in >90% of the samples. Furthermore, since the change to PA as first line treatment in Gia Lai province in 2020 (and TES trials in 2017–2018), we start to observe the fitness compensatory mutation C258W in Pfcrt, which reportedly reverses the PPQ resistance caused by the F145I mutation while it increases CQ resistance (Hagenah et al., 2024). This C258W mutation has been shown through competitive growth assays with pfcrt-edited parasites to reduce the fitness defect caused by F145I, observed at increasing levels in this study (48/354 samples, with 30 samples from Binh Thuan), and has been reported to emerge in laboratory cultures and Southeast Asian isolates no longer exposed to PPQ pressure (Hagenah et al., 2024). Only few samples were observed with this mutation, and future studies including samples from the years after the treatment guideline change are necessary to assess the spread and importance of this mutation. In addition, we identified several non-PfK13 mutations that may contribute to artemisinin resistance or enhanced parasite fitness of mutant parasites, as reported in Thailand and reflecting similar findings in African settings where additional resistance mechanisms have reduced ACT efficacy (Cerqueira et al., 2017). We detected E1528D in Pfubp1, associated with artemisinin resistance, in Gia Lai and Kon Tum, in parasites also carrying PfK13 C580Y. In samples with the PfK13 P553L mutation from Khanh Hoa, a unique insertion in Pfubp1 was observed alongside mutations in the Pfeps15 gene, suggesting a complex genetic landscape influencing artemisinin resistance. Further investigation into the functional significance and epidemiological relevance of these non-PfK13 markers and in relation to the evolution of PfK13-mediated ART-R is needed, as it could uncover important insights into the evolving landscape of antimalarial drug resistance.

The high prevalence of multidrug resistance, including resistance to ART, PPQ, CQ, and SP, poses a significant threat to the current malaria control and elimination strategies in Vietnam. This trend underscores the need for continuous monitoring and potential revision of treatment guidelines, as evidenced by the switch to PA for uncomplicated malaria. High levels of CQ resistance in P. falciparum from Khanh Hoa, Binh Phuoc, and Gia Lai calls for close monitoring of CQ efficacy for P. vivax, especially in view of the recent outbreaks in Lai Chau and Khanh Hoa. Malaria elimination faces significant challenges, but genetic surveillance can guide interventions like mass drug administration (MDA), targeted drug administration (TDA), or enhanced vector control. For instance, Vietnam has implemented MDA in Lai Chau and plans TDA in Khanh Hoa due to the outbreaks. Molecular evaluations before and after these interventions are crucial, as past experiences indicate MDA and TDA may not offer sustained results, however this is not currently part of the strategy. Based on the high level of resistance and its continuous evolution detected in this study, it is recommended to include genetic surveillance in these outbreak scenarios.

Pyronaridine resistance in Plasmodium has been linked to high mdr1 expression in mouse models (Kimani et al., 2014; Kimani and Shume, 2020). Although we did not investigate mdr1 expression, we found very few samples with mdr1 gene amplifications. While the resistance mechanism is not fully understood, pyronaridine resistance shares similarities with resistance to other quinolines (Dorn et al., 1998). Given the high resistance rates to PPQ and other quinolines in Vietnam, PA might not be ideal, and an alternative ACT may be needed for effective malaria management.

On the other hand, there was a notable lack of Pfhrp2/3 gene deletions, which can lead to false-negative rapid diagnostic test results, consistent with other studies in the region (Molina-de La Fuente et al., 2021; Thang et al., 2022). This contrasts with certain countries in South America and the Horn of Africa, where Pfhrp2/3 deletions are prevalent and prompted changes to RDT strategies (Molina-de La Fuente et al., 2021). The absence of such deletions in the study samples indicates that the current HRP2/3-based RDT remains suitable for the diagnosis of P. falciparum in the country.

Plasmodium falciparum populations in Vietnam are becoming more isolated through time, as indicated by the high proportion of clonal populations, with a high degree of geographical clustering by province, indicating limited connectivity between parasite populations across different provinces. These distinct, localized subpopulations may have evolved independently in response to an interplay of local environmental factors, epidemiological patterns, and population movements. Interestingly, the Central Highland provinces, and especially Gia Lai province, is the region with the highest residual burden of malaria in Vietnam, while it had the lowest diversity and very clonal population. This is likely due to the drug resistance profile of the parasites in this province, combined with strong selection pressure exerted/posed by the targeted efforts from the NCMP to reduce malaria in this province. The observed geographical isolation of these subpopulations is important for malaria elimination efforts, as it can promote the development and spread of region-specific drug resistance profiles that may require tailored intervention strategies. The NMCP efforts in Vietnam are organized however on a provincial level, and this localized approach might be one of the factors contributing to the success of the country in reducing its burden. Similar trends of geographically isolated parasite populations have been reported in Vietnam and the Greater Mekong Subregion, indicating a regional pattern of fragmented malaria transmission (Verschuuren et al., 2024; Cui et al., 2012).

The comprehensive approach of this study provides a robust framework for understanding drug resistance and its evolution. By examining a broad panel of full-length resistance genes spanning multiple drug classes and a wide range of SNPs, it provided a high-resolution perspective on the genetic diversity and resistance evolution in this country. This resulted in valuable insights that can inform local policy and decision makers to improve treatment and intervention strategies. However, the limited geographic scope and sample size highlights the need for broader surveillance in more malaria endemic provinces to capture the full spectrum of malaria resistance in Vietnam, including for example, the sites of recent outbreaks. Low incidence in some provinces made it difficult to collect large numbers of samples for statistical analysis. However, the few samples that we did include often represented >10% of the reported malaria cases in the province, making it still a meaningful assessment to inform NMCPs. In contrast, TES in these settings would not be feasible due to the low incidence, demonstrating the added value of molecular surveillance in a pre-elimination scenario. However, for genetic surveillance data to be actionable, these results need to be produced much faster and become integrated with routine epidemiological and entomological monitoring, to allow for a comprehensive understanding of malaria transmission and control. Moreover, a centralized database for sharing genetic surveillance data would facilitate coordinated responses and efficient resource allocation. Enhancing collaboration with local health authorities and communities is essential for a more representative data collection and swift outbreak response, supporting the country’s efforts towards malaria elimination.

The coordinated implementation of various approaches, including targeted drug administration, vector control measures, active surveillance and community engagement strategies, is crucial for malaria elimination efforts. However, financial constraints and changes in procurement regulations have made it challenging to obtain the necessary chemicals, antimalarial drugs, and to mobilize sufficient resources for effective implementation of these interventions and surveillance activities. This lack of coordination and resource availability can hinder the synchronization of the various control measures, thereby increasing the risk of the further spread and proliferation of drug-resistant malaria strains across the region. Addressing these systemic challenges through increased and sustained funding, streamlined procurement processes, and a comprehensive, multi-faceted strategy is essential to effectively combat the escalating threat of drug-resistant malaria and achieve the goal of malaria elimination in Vietnam.